3

I was wondering of there is a QGIS equivalent of the following Turf.js nearest point function, which can be called from a Python Script.

https://www.npmjs.com/package/@turf/nearest-point

The JavaScript function takes a reference point and a FeatureCollection of Features with Point geometries and returns the point from the FeatureCollection closest to the reference.

I am trying to write a Python equivalent of the following JavaScript block of code.

nearestPoint = turf.nearestPoint(currentFeature, entire_dataset)

As it happens, in QGIS, there is not exactly the same function but there is a Distance to nearest hub (points) tool in the processing toolbox that will give the distance to all of the points and then you can choose the one with the shortest distance. We can also use the select by expression using the Closest Point function.

However, these QGIS features need to be invoked manually, and I was hoping that these functions (or a similar one) could somehow be called using a Python script from within QGIS. Maybe someone with knowledge of QGIS internal functions can advice on how these functions can be called.

1

1 Answer 1

4

You can use the QgsSpatialIndex class which has a nearestNeighbor() method.

E.g. to iterate over a point layer, finding the nearest point and distance for each feature, you could do something like:

project = QgsProject.instance()

lyr = project.mapLayersByName('Your Layer Name')[0]

# Create a QgsSpatialIndex instance & load layer features
sp_idx = QgsSpatialIndex(lyr.getFeatures())

# Iterate layer features
for ft in lyr.getFeatures():
    print(f'Current ID: {ft.id()}')
    current_geom = ft.geometry()
    # Remove current feature so nn search doesn't just return itself
    sp_idx.deleteFeature(ft)
    # Find nearest point to current geom (returns feature id)
    # (default is 1 but you can add additional arguments for number of neighbors/max search distance)
    nn_id = sp_idx.nearestNeighbor(current_geom)
    print(f'Nearest neighbor ID: {nn_id}')
    nn = lyr.getFeature(nn_id[0])
    nn_geom = nn.geometry()
    # We can also find the distance to the current point
    # This will return cartesian distance in the layer's CRS units
    dist = current_geom.distance(nn_geom)
    print(f'Distance to nearest neighbor: {dist}')
    # Insert the current feature back into the spatial index
    sp_idx.insertFeature(ft)

As commented in the above example, the distance calculation is cartesian, performed in the layer's crs. To calculate ellipsoidal distance, you can use the QgsDistanceArea class.

project = QgsProject.instance()

lyr = project.mapLayersByName('Your Layer Name')[0]

# Create a QgsSpatialIndex instance & load layer features
sp_idx = QgsSpatialIndex(lyr.getFeatures())

# Create a QgsDistanceArea instance, set CRS & ellipsoid
da = QgsDistanceArea()
da.setSourceCrs(lyr.crs(), project.transformContext())
da.setEllipsoid(lyr.crs().ellipsoidAcronym())

# Iterate layer features
for ft in lyr.getFeatures():
    print(f'Current ID: {ft.id()}')
    current_geom = ft.geometry()
    # Remove current feature so nn search doesn't just return itself
    sp_idx.deleteFeature(ft)
    nn_id = sp_idx.nearestNeighbor(current_geom)
    print(f'Nearest neighbor ID: {nn_id}')
    nn = lyr.getFeature(nn_id[0])
    nn_geom = nn.geometry()
    # Use QgsDistanceArea to calculate ellipsoidal distance to nearest point
    # (Distance returned is in meters)
    dist = da.measureLine(current_geom.asPoint(), nn_geom.asPoint())
    print(f'Distance to nearest neighbor: {round(dist, 2)}m')
    # Insert the current feature back into the spatial index
    sp_idx.insertFeature(ft)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.