9

Given a point and few lines, how should I go about finding the nearest line to the point?

enter image description here

I am aware of Finding nearest line to point using ArcGIS Desktop (ArcObjects/ArcPy)?, but that is in ArcGIS and uses its functions.


The "duplicate" question doesn't answer this question properly, it's just a pointer to MMQGIS.

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
PeterBorook
  • 639
  • 6
  • 11
  • 1
    After you've found the nearest line, what do you want to do with it? – R.K. Apr 26 '13 at 02:03
  • @R.K. Then I want to find the perpendicular distance from the point to that line – PeterBorook Apr 26 '13 at 02:21
  • 2
    There has been some solution in QGIS mailing list a few years ago:http://lists.osgeo.org/pipermail/qgis-user/2010-May/008745.html. Carson provide "meta-python-code": – posiki Apr 26 '13 at 05:21
  • The NNJoin QGIS Plugin can help you with the job. For each feature of the input layer it adds all the attributes of the nearest feature in the join layer and also adds a distance attribute with the distance to this feature. So if you only want to join the attributes of the nearest line to each point, the NNJoin result is all you need. Otherwise, the line can be identified through its attributes. – Håvard Tveite Mar 25 '17 at 11:04

4 Answers4

15

You can use QgsSpatialIndex class for finding nearest objects. First you will need to create a new object of this class. Then add the required features to the index. Then you should be able to use QgsSpatialIndex.nearestNeighbor (QgsPoint point, int neighbors) method to retrieve the nearest ones.

Here is sample that I managed in python console.

lineLayer = iface.activeLayer()
provider = lineLayer.dataProvider()

spIndex = QgsSpatialIndex() #create spatial index object

feat = QgsFeature() fit = provider.getFeatures() #gets all features in layer

insert features to index

while fit.nextFeature(feat): spIndex.insertFeature(feat)

pt = QgsPoint(-0.00062201,0.00001746)

QgsSpatialIndex.nearestNeighbor (QgsPoint point, int neighbors)

nearestIds = spIndex.nearestNeighbor(pt,1) # we need only one neighbour print nearestIds

Edit:

To get the actual QgsFeature object from the python list, you can do this,

featureId = nearestIds[0]
fit2 = lineLayer.getFeatures(QgsFeatureRequest().setFilterFid(featureId))
ftr = QgsFeature()
fit2.nextFeature(ftr)
# ftr now contains the QgsFeature object for the id
Matt
  • 16,843
  • 3
  • 21
  • 52
vinayan
  • 7,282
  • 3
  • 37
  • 75
  • Thanks for introducing me to QgsSpatialIndex! it all looks wonderful; just let me ask a further question: How can I get the actual feature in "lineLayer" from this "nearestIds"?( i.e. get featureid from indexid?) – PeterBorook Apr 26 '13 at 06:50
  • @PeterBorook - i updated the answer..the code is based on version 1.9. If you need to get it working for older version, please refer http://www.qgis.org/pyqgis-cookbook/ – vinayan Apr 26 '13 at 07:19
  • 5
    As far as I know, QgsSpatialIndex.nearestNeighbor finds the nearest neighbours among the index geometries. The index consists of approximations of the real geometries (bounding box for lines and polygons), so the nearest neighbours found in the index are not necessarily the real nearest neighbour. – Håvard Tveite Sep 23 '14 at 06:59
4

Almost 10 years later it is probably worth adding an update to @vinayan's answer:

By storing the features geometries inside the QgsSpatialIndex() using QgsSpatialIndex.FlagStoreFeatureGeometries flag, you can get the real closest features, not just the approximate ones by their bounding box. This adresses the comment by Håvard Tveite. Also, nowadays it is much more efficient to bulkload the index at once; this is much faster than adding features one by one after the creation of the index in a loop. Additionally you can get the nearest geometries directly from the index, which is much faster than getting the whole feature via a QgsFeatureRequest(). If you still want to get the whole feature, you can now also use directly .getFeature(id), which is shorter. Apart from that, the syntax is now for PyQGIS 3. The code:

layer = iface.activeLayer()

Bulkloading the index is much faster than adding features one by one

Storing the geometry in the index allows finding the exact nearest neighbor, also for not single-point-features

index = QgsSpatialIndex(layer.getFeatures(), flags=QgsSpatialIndex.FlagStoreFeatureGeometries)

QgsSpatialIndex.nearestNeighbor(QgsPointXY(), int neighbors)

pt = QgsPointXY(-0.00062201,0.00001746) nearestIds = index.nearestNeighbor(pt,1) # we need only one neighbour

Or, not stated in the docs, if you store the geometry in the index as done here, you can also: QgsSpatialIndex.nearestNeighbor(QgsGeometry(), int neighbors)

geom = QgsGeometry.fromPointXY(QgsPointXY(-0.00062201,0.00001746)) nearestIds = index.nearestNeighbor(geom,1)

note that in reality, if both layers are the same, the nearest feature will always be itself

in that case you should request 2 nearest neighbors and remove the one with the same id as your current source feature

print(nearestIds)

for nearestId in nearestIds: # get the geometry of nearest feature efficiently directly from the index without a slow feature request nearestGeom = index.geometry(nearestId) print(nearestGeom) # or get the whole feature via a short getFeature() request nearestFeat = layer.getFeature(nearestId) print(nearestFeat)

Glorfindel
  • 1,096
  • 2
  • 9
  • 14
MrXsquared
  • 34,292
  • 21
  • 67
  • 117
  • Thanks! Without this flag, only lines bounding boxes are considered. I'm sure it should be faster, but gives the wrong results sometimes – Andrii Liekariev Nov 17 '23 at 21:41
2

Meanwhile, there is an easy way to do that using QGIS expressions with overlay_nearest() function:

overlay_nearest ('line_layer', $geometry)[0]
Babel
  • 71,072
  • 14
  • 78
  • 208
2

Since QGIS 3.8 there is also a native processing algorithm: Join Attributes By Nearest

enter image description here

MrXsquared
  • 34,292
  • 21
  • 67
  • 117