Nearest neighbor between a point layer and a line layer

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?

EDIT

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? – WolfOdrade Jul 27 '10 at 19:22
@wolfOdrade - Overall, it doesn't matter. – dassouki Jul 27 '10 at 19:33

7 Answers

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 '13 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 '13 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 '13 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 '13 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 '13 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
``````
-

What kind of use case do you have?

Depending on your interests, 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 – dassouki Jul 27 '10 at 19:00
Ah, so your goal is to automate creation of this "dummy road link"? – underdark Jul 27 '10 at 19:42
indeed :) or dummy link(s) – dassouki Jul 27 '10 at 19:52

Another PostGIS answer

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.

HTH Nicklas

-

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 '11 at 5:27

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)
``````
-

Our RW Net 4 product has highly optimized algorithms for such situations, due to a very efficient spatial index, created for this in particular. Ideal if you have to run the calculations many times (even billions!), for GPS coordinates for instance. It returns distance, coordinates of projection, ID of line, side of line, distance from start of line etc.

-