python - Nearest neighbor between point layer and 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?

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

Answer

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...