# Nearest neighbor between point layer and line layer? [closed]

I've asked this question several times on stackoverflow and irc between #qgis and #postgis and I also tried to code it or implement it my self in postgis with no real answer.

Using programming (most preferably python), I'd like to draw a line from a point layer, to its projection on the nearest line of a line or polygon layer.

As of now most of my data is in ESRI's shape and postgis formats; however, I'd rather stay away from a postgis solution as I'm predominantly a shp + qgis user.

An ideal solution would be to implement GDAL/OGR with python or similar libraries

• Using the GDAL/OGR libraries where should I start? would it be possible to give a solution plan?
• Can I use NetworkX to do the nearest neighbor analysis?
• Is this actually possible?

If it's easier, the points could connect to the segment end point instead of a projected point

• can the line be restricted to being orthagonal to the line segment? Jul 27, 2010 at 19:22
• @wolfOdrade - Overall, it doesn't matter. Jul 27, 2010 at 19:33

This question turned out to be a bit trickier than I thought to get right. There are many implementations of the shortest distance itself, such as the Shapely provided distance (from GEOS). Few of the solutions provide the intersection point itself, however, but only the distance.

My first attempt buffered the point by the distance between the point and polygon, and looked for intersections, but rounding errors prevent this from giving an exact answer.

Here's a complete solution using Shapely, based on these equations:

``````#!/usr/bin/env python
from shapely.geometry import Point, Polygon
from math import sqrt
from sys import maxint

# define our polygon of interest, and the point we'd like to test
# for the nearest location
polygon = Polygon(((0, 0), (0, 1), (1, 1), (1, 0), (0, 0)))
point = Point(0.5, 1.5)

# pairs iterator:
# http://stackoverflow.com/questions/1257413/1257446#1257446
def pairs(lst):
i = iter(lst)
first = prev = i.next()
for item in i:
yield prev, item
prev = item
yield item, first

# these methods rewritten from the C version of Paul Bourke's
# geometry computations:
# http://local.wasp.uwa.edu.au/~pbourke/geometry/pointline/
def magnitude(p1, p2):
vect_x = p2.x - p1.x
vect_y = p2.y - p1.y
return sqrt(vect_x**2 + vect_y**2)

def intersect_point_to_line(point, line_start, line_end):
line_magnitude =  magnitude(line_end, line_start)
u = ((point.x - line_start.x) * (line_end.x - line_start.x) +
(point.y - line_start.y) * (line_end.y - line_start.y)) \
/ (line_magnitude ** 2)

# closest point does not fall within the line segment,
# take the shorter distance to an endpoint
if u < 0.00001 or u > 1:
ix = magnitude(point, line_start)
iy = magnitude(point, line_end)
if ix > iy:
return line_end
else:
return line_start
else:
ix = line_start.x + u * (line_end.x - line_start.x)
iy = line_start.y + u * (line_end.y - line_start.y)
return Point([ix, iy])

nearest_point = None
min_dist = maxint

for seg_start, seg_end in pairs(list(polygon.exterior.coords)[:-1]):
line_start = Point(seg_start)
line_end = Point(seg_end)

intersection_point = intersect_point_to_line(point, line_start, line_end)
cur_dist =  magnitude(point, intersection_point)

if cur_dist < min_dist:
min_dist = cur_dist
nearest_point = intersection_point

print "Closest point found at: %s, with a distance of %.2f units." % \
(nearest_point, min_dist)
``````

For posterity, it looks like this ArcView extension handles this problem quite nicely, too bad its on a dead platform written in a dead language...

• I wonder if there is a way to index polygon points to avoid explicit enumeration...
– mlt
May 15, 2013 at 0:04
• @mlt not sure exactly what you're thinking of, but there are some approaches that can help depending on the geometry. Could do some basic ray-casting to determine relevant nearest segments, if performance was an issue. In that vein, moving this into C or Pyrex would improve things.
– scw
May 15, 2013 at 11:24
• I mean with `pairs` it is algorithmically O(n) or something. @eprand solution perhaps can be modified to use KNN however I managed to live without PostGIS so far...
– mlt
May 15, 2013 at 17:26
• I can't edit my previous comment any longer:( Perhaps solution of Nicklas Avén with ST_Closestpoint & ST_Shortestline are the fastest if PostGIS is an option.
– mlt
May 15, 2013 at 17:40
• Right, you could use a KNN algorithm in Python directly. I don't believe that ST_Shortestline uses KNN, it just iterates as well based on my reading of postgis.refractions.net/documentation/postgis-doxygen/d1/dbf/…
– scw
Jun 4, 2013 at 5:17

A PostGIS answer (for multilinestring, if linestring, remove st_geometryn function)

``````select t2.gid as point_gid, t1.gid as line_gid,
st_makeline(t2.geom,st_line_interpolate_point(st_geometryn(t1.geom,1),st_line_locate_point(st_geometryn(t1.geom,1),t2.geom))) as geom
from your_line_layer t1, your_point_layer t2,
(
select gid as point_gid,
(select gid
from your_line_layer
order by st_distance(your_line_layer.geom, your_point_layer.geom)
limit 1 ) as line_gid
from your_point_layer
) as t3
where t1.gid = t3.line_gid
and t2.gid = t3.point_gid
``````

This is a bit old, but I was searching for solutions to this problem today (point --> line). The simplest solution I've come across for this related problem is:

``````>>> from shapely.geometry import Point, LineString
>>> line = LineString([(0, 0), (1, 1), (2, 2)])
>>> point = Point(0.3, 0.7)
>>> point
POINT (0.3000000000000000 0.7000000000000000)
>>> line.interpolate(line.project(point))
POINT (0.5000000000000000 0.5000000000000000)
``````

If I understand you right the functionality you are asking for is built in PostGIS.

To get a point projected on a line you can use ST_Closestpoint (on PostGIS 1.5)

Some hints about how to use it you can read here: http://blog.jordogskog.no/2010/02/07/how-to-use-the-new-distance-related-functions-in-postgis-part1/

It is usable also to find the closest point on a polygon to another polygon for instance.

If you want the line between the two closest points on both geometries you can use ST_Shortestline. ST_Closestpoint is the first point in ST_Shortestline

The length of ST_Shortestline between two geometries is the same as ST_Distance between the geometries.

See the comment below concerning how my answer should not be considered a reliable solution... I'll leave this original post here just so others can examine the problem.

If I understand the question, this general procedure should work.

To find the shortest path between a point (as defined by x,y or x,y,z) and a polyine (as defined by a connecting set of x,y or x,y,z) within Euclidean space:

1) From a given user-defined point (I'll call it pt0), find the nearest vertex of polyline (pt1). OGRinfo can poll the vertices of a polyline, and then distance calculations can be made via standard methods. For example, iterate over a distance calc like: distance_in_radians=2*math.asin(math.sqrt(math.pow((math.sin((pt0_radians-ptx_radians)/2)),2) + math.cos(pt0_radians)*math.cos(ptx_radians)*math.pow((math.sin((pt0_radians-ptx_radians)/2)),2)))

2) Store the associated minimum distance value (d1) and (pt1)

3) look at the two segments stemming away from pt1 (in the ogrinfo linestring, these will be the prior and subsequent vertices). Record these vertices (n2 and n3).

4) create y = mx + b formula for each segment

5) Relate your point (pt0) to the perpendicular for each of those two formulas

6) Calculate distance and intersections (d2 and d3; pt2, pt3)

7) Compare the three distances (d1, d2, d3) for the shortest. Your pt0 to the associated node (pt1, pt2, or pt3) is the shortest link.

That's a stream of consciousness answer -- hopefully, my mental picture of the problem and solution fits.

• This won't work in general. E.g. point = (1,1), line = ((0,2),(0,3),(3,0),(2,0)). If you sketch that, you can see the "closest" vertices on the line are not adjacent to the segment that passes closest to the point... I think the only way to handle this is check every segment (possibly using bounding boxes to avoid optimize it a bit). HTH.
– Tom
Mar 22, 2011 at 5:27

Here is a python script for QGIS >2.0 made from the hints and solutions given above. It works fine for a reasonable amount of points and lines. But I did not try it with huge amount of objects.

Of course it had to be copy in idle or whatever other less "pythonic solution" and save it as "closest.point.py".

In QGIS tool box go for script,tools, add a script, an choose it.

``````##Vector=group
##CLosest_Point_V2=name
##Couche_de_Points=vector
##Couche_de_Lignes=vector

"""
This script intent to provide a count as for the SQL Funciton CLosestPoint
Ce script vise a recrÃ©er dans QGIS la Focntion SQL : CLosest Point
It rely on the solutions provided in "Nearest neighbor between a point layer and a line layer"
http://gis.stackexchange.com/questions/396/nearest-pojected-point-from-a-point-                               layer-on-a-line-or-polygon-outer-ring-layer
V2 du  8 aout 2016
[email protected]
"""
from qgis.core import *
from qgis.gui import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
import os
import sys
import unicodedata
from osgeo import ogr
from math import sqrt
from sys import maxint

from processing import *

def magnitude(p1, p2):
if p1==p2: return 1
else:
vect_x = p2.x() - p1.x()
vect_y = p2.y() - p1.y()
return sqrt(vect_x**2 + vect_y**2)

def intersect_point_to_line(point, line_start, line_end):
line_magnitude =  magnitude(line_end, line_start)
u = ((point.x()-line_start.x())*(line_end.x()-line_start.x())+(point.y()-line_start.y())*(line_end.y()-line_start.y()))/(line_magnitude**2)
# closest point does not fall within the line segment,
# take the shorter distance to an endpoint
if u < 0.0001 or u > 1:
ix = magnitude(point, line_start)
iy = magnitude(point, line_end)
if ix > iy:
return line_end
else:
return line_start
else:
ix = line_start.x() + u * (line_end.x() - line_start.x())
iy = line_start.y() + u * (line_end.y() - line_start.y())
return QgsPoint(ix, iy)

layerP = processing.getObject(Couche_de_Points)
providerP = layerP.dataProvider()
fieldsP = providerP.fields()
inFeatP = QgsFeature()

layerL = processing.getObject(Couche_de_Lignes)
providerL = layerL.dataProvider()
fieldsL = providerL.fields()
inFeatL = QgsFeature()

counterP = counterL= nElement=0

for featP in layerP.selectedFeatures():
counterP+=1
if counterP==0:
QMessageBox.information(None,"information:","Choose at least one point from point layer_"+ str(layerP.name()))

indexLine=QgsSpatialIndex()
for featL in layerL.selectedFeatures():
indexLine.insertFeature(featL)
counterL+=1
if counterL==0:
QMessageBox.information(None,"information:","Choose at least one line from point layer_"+ str(layerL.name()))
#QMessageBox.information(None,"DEBUGindex:",str(indexBerge))
ClosestP=QgsVectorLayer("Point", "Projected_ Points_From_"+ str(layerP.name()), "memory")
prClosestP = ClosestP.dataProvider()

for f in fieldsP:
znameField= f.name()
Type= str(f.typeName())
if Type == 'Integer': prClosestP.addAttributes([ QgsField( znameField, QVariant.Int)])
if Type == 'Real': prClosestP.addAttributes([ QgsField( znameField, QVariant.Double)])
if Type == 'String': prClosestP.addAttributes([ QgsField( znameField, QVariant.String)])
else : prClosestP.addAttributes([ QgsField( znameField, QVariant.String)])
QgsField("XDep", QVariant.Double),
QgsField("YDep", QVariant.Double),
QgsField("XProj", QVariant.Double),
QgsField("YProj", QVariant.Double),
QgsField("Xmed", QVariant.Double),
QgsField("Ymed", QVariant.Double)])
featsP = processing.features(layerP)
nFeat = len(featsP)
"""
for inFeatP in featsP:
progress.setPercentage(int(100 * nElement / nFeatL))
nElement += 1
# pour avoir l'attribut d'un objet/feat ....
attributs = inFeatP.attributes()
"""

for inFeatP in layerP.selectedFeatures():
progress.setPercentage(int(100 * nElement / counterL))
nElement += 1
attributs=inFeatP.attributes()
geomP=inFeatP.geometry()
nearest_point = None
minVal=0.0
counterSelec=1
first= True
nearestsfids=indexLine.nearestNeighbor(geomP.asPoint(),counterSelec)
#layer.getFeatures( QgsFeatureRequest().setFilterFid( fid ) )
request = QgsFeatureRequest().setFilterFids( nearestsfids )
#list = [ feat for feat in CoucheL.getFeatures( request ) ]
# QMessageBox.information(None,"DEBUGnearestIndex:",str(list))
NBNodes=0
Dist=DistT=minValT=Distance=0.0
for featL in  layerL.getFeatures(request):
geomL=featL.geometry()
firstM=True
geomL2=geomL.asPolyline()
NBNodes=len(geomL2)
for i in range(1,NBNodes):
lineStart,lineEnd=geomL2[i-1],geomL2[i]
ProjPoint=intersect_point_to_line(geomP.asPoint(),QgsPoint(lineStart),QgsPoint(lineEnd))
Distance=magnitude(geomP.asPoint(),ProjPoint)
toto=''
toto=toto+ 'lineStart :'+ str(lineStart)+ '  lineEnd : '+ str(lineEnd)+ '\n'+ '\n'
toto=toto+ 'ProjPoint '+ str(ProjPoint)+ '\n'+ '\n'
toto=toto+ 'Distance '+ str(Distance)
#QMessageBox.information(None,"DEBUG", toto)
if firstM:
minValT,nearest_pointT,firstM = Distance,ProjPoint,False
else:
if Distance < minValT:
minValT=Distance
nearest_pointT=ProjPoint
#at the end of the loop save the nearest point for a line object
#min_dist=magnitude(ObjetPoint,PProjMin)
#QMessageBox.information(None,"DEBUG", " Dist min: "+ str(minValT))
if first:
minVal,nearest_point,first = minValT,nearest_pointT,False
else:
if minValT < minVal:
minVal=minValT
nearest_point=nearest_pointT
#at loop end give the nearest Projected points on Line nearest Line
PProjMin=nearest_point
Geom= QgsGeometry().fromPoint(PProjMin)
min_dist=minVal
PX=geomP.asPoint().x()
PY=geomP.asPoint().y()
Xmed=(PX+PProjMin.x())/2
Ymed=(PY+PProjMin.y())/2
newfeat = QgsFeature()
newfeat.setGeometry(Geom)
Values=[]
#Values.extend(attributs)
fields=layerP.pendingFields()
Values=[attributs[i] for i in range(len(fields))]
Values.append(min_dist)
Values.append(PX)
Values.append(PY)
Values.append(PProjMin.x())
Values.append(PProjMin.y())
Values.append(Xmed)
Values.append(Ymed)
newfeat.setAttributes(Values)
ClosestP.startEditing()
#prClosestP.updateExtents()
ClosestP.commitChanges()
iface.mapCanvas().refresh()
``````

!!! WARNING !!! Pay attention that some "weird"/ wrong projected points may be produce because of this line command:

``````nearestsfids=indexLine.nearestNeighbor(geomP.asPoint(),counterSelec)
``````

The `counterSelec` value in it set how many nearestNeighbor are returned. In fact each point should be projected at the shortest distance possible to each line object; and the minimum distance found would give the correct line and projected point as the nearest neighbours we seek. In order to reduce the looping time, the nearest Neighbour command is used. Choosing a `counterSelec` value reduce to 1 will return the "first" object met (it's bounding box more exactly) and it may not be the right one. Different line size objects may oblige to choose may be 3 or 5, or even more nearest objects in order to determine the shortest distance. The higher the value, the longer it takes. With hundreds of points and lines it start to get very slow with 3 or 5 nearest neighbour, with thousands it may bug with such values.

Depending on your interests and use case, it might be useful to look into "map matching algorithms". For example, there is a RoadMatcher project on OSM wiki: http://wiki.openstreetmap.org/wiki/Roadmatcher.

• It's for travel demand and forecasting. Usually we divide areas into traffic analysis zones (polygons) and we establish the centroid of the polygon as the "dummy" originator of all traffic in that zone. We then draw x or y "dummy road link" lines from that point to the nearest roads and distribute the traffic equally from that zone onto those dummy links and onto the actual road layer Jul 27, 2010 at 19:00
• Ah, so your goal is to automate creation of this "dummy road link"? Jul 27, 2010 at 19:42
• indeed :) or dummy link(s) Jul 27, 2010 at 19:52