Find shortest path from one line to other in Shapely

Question

Given two lines line_a and line_b, how do I find the pair of points that represents the smaller route from line_a to line_b?

Answer 1

Unfortunately, there is no shapely operation to do this. I was thinking about this problem when was asked to extend the solution in the question find-coordinate-of-closest-point-on-polygon-shapely , to handle two polygons. The answer depends on geometry theorem that is intuitively true (although a formal proof requires some calculus and it is a little long to be written here without LaTex). The theorem tells:

The minimum distance between two line strings that do not cross each other (being a line string a concatenated sequence of segments), is always achieved in one of the edge points of the line strings.

With this in mind, the problem is reduced to compute the minimum distance between each edge of a LineString to the other LineString.

Another way of seeing this problem is to reduce it to compute the minimum distance between each pair of segments. And to notice that the distance between two segments that do not intersect, is achieved between two end points, or between an endpoint and the projection of that end point over the other segment.

The code has been optimized a little to avoid redundant computation, but perhaps you can get a more elegant version. Or if you are familiar with numpy, you can probably get a shorter version that use numpy vector distance and dot product.

Beware, If you are handling several thousands of points in the polygons, this routine should be optimized to avoid the computation of all distances between edges. Perhaps, while computing edge to edge distance, you could discard edge by introducing some clever filtering.

from shapely.geometry import LineString, Point
import math 

def get_min_distance_pair_points(l1, l2):
    """Returns the minimum distance between two shapely LineStrings.

    It also returns two points of the lines that have the minimum distance.
    It assumes the lines do not intersect.

    l2 interior point case:
    >>> l1=LineString([(0,0), (1,1), (1,0)])
    >>> l2=LineString([(0,1), (1,1.5)])
    >>> get_min_distance_pair_points(l1, l2)
    ((1.0, 1.0), (0.8, 1.4), 0.4472135954999578)

    change l2 slope to see the point changes accordingly:
    >>> l2=LineString([(0,1), (1,2)])
    >>> get_min_distance_pair_points(l1, l2)
    ((1.0, 1.0), (0.5, 1.5), 0.7071067811865476)

    l1 interior point case:
    >>> l2=LineString([(0.3,.1), (0.6,.1)])
    >>> get_min_distance_pair_points(l1, l2)
    ((0.2, 0.2), (0.3, 0.1), 0.1414213562373095)

    Both edges case:
    >>> l2=LineString([(5,0), (6,3)])
    >>> get_min_distance_pair_points(l1, l2)
    ((1.0, 0.0), (5.0, 0.0), 4.0)

    Parallels case:
    >>> l2=LineString([(0,5), (5,0)])
    >>> get_min_distance_pair_points(l1, l2)
    ((1.0, 1.0), (2.5, 2.5), 2.1213203435596424)

    Catch intersection with the assertion:
    >>> l2=LineString([(0,1), (1,0.8)])
    >>> get_min_distance_pair_points(l1, l2)
    Traceback (most recent call last):
      ...
      assert( not l1.intersects(l2))
    AssertionError

    """ 

    def distance(a, b):
        return math.sqrt( (a[0]-b[0])**2 + (a[1]-b[1])**2 ) 

    def get_proj_distance(apoint, segment):
        '''
        Checks if the ortogonal projection of the point is inside the segment.

        If True, it returns the projected point and the distance, otherwise 
        returns None.
        '''
        a = ( float(apoint[0]), float(apoint[1]) )
        b, c = segment
        b = ( float(b[0]), float(b[1]) )
        c = ( float(c[0]), float(c[1]) )
        # t = <a-b, c-b>/|c-b|**2
        # because p(a) = t*(c-b)+b is the ortogonal projection of vector a 
        # over the rectline that includes the points b and c. 
        t = (a[0]-b[0])*(c[0]-b[0]) + (a[1]-b[1])*(c[1]-b[1])
        t = t / ( (c[0]-b[0])**2 + (c[1]-b[1])**2 )
        # Only if t 0 <= t <= 1 the projection is in the interior of 
        # segment b-c, and it is the point that minimize the distance 
        # (by pitagoras theorem).
        if 0 < t < 1:
            pcoords = (t*(c[0]-b[0])+b[0], t*(c[1]-b[1])+b[1])
            dmin = distance(a, pcoords)
            return a, pcoords, dmin
        elif t <= 0:
            return a, b, distance(a, b)
        elif 1 <= t:
            return a, c, distance(a, c)

    def get_min(items1, items2, distance_func, revert=False):
        "Minimum of all distances (with points) achieved using distance_func."
        a_min, b_min, d_min = None, None, None
        for p in items1:
            for s in items2:
                a, b, d = distance_func(p, s)
                if d_min == None or d < d_min:
                    a_min, b_min, d_min = a, b, d 
        if revert:
            return b_min, a_min, d_min
        return a_min, b_min, d_min

    assert( not l1.intersects(l2))

    l1p = list(l1.coords)
    l2p = list(l2.coords)
    l1s = zip(l1p, l1p[1:])
    l2s = zip(l2p, l2p[1:])

    edge1_min = get_min(l1p, l2s, get_proj_distance)
    edge2_min = get_min(l2p, l1s, get_proj_distance, revert=True)

    if edge1_min[2] <= edge2_min[2]:
        return edge1_min
    else:
        return edge2_min

if __name__ == "__main__":
    import doctest
    doctest.testmod()