Nearest neighbor between a point layer and a line layer

nn.line.point.py

#!/usr/bin/env python
#Author: scw
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)