在 Shapely 中找到从一条线到另一条线的最短路径

Question

给定两条线line_a 和line_b，我如何找到代表从line_a 到line_b 的较小路线的点对？

Answer 1

不幸的是，没有合适的操作可以做到这一点。 当被要求扩展解决方案时，我正在考虑这个问题 问题find-coordinate-of-closest-point-on-polygon-shapely，处理两个多边形。答案取决于几何定理 这在直觉上是正确的（尽管正式证明需要一些微积分 没有 LaTex 在这里写有点长）。 定理告诉我们：

两个不相交的线串之间的最小距离 其他（作为线串连接的段序列）， 总是在线串的边缘点之一实现。

考虑到这一点，问题被简化为计算最小距离 一个 LineString 的每个边缘到另一个 LineString 之间。

查看此问题的另一种方法是减少它以计算最小值 每对线段之间的距离。并注意到距离 在两个不相交的线段之间，在两个端点之间实现， 或在一个端点和该端点在另一个端点上的投影之间 分割。

代码做了一些优化以避免冗余计算，但是 也许你可以得到一个更优雅的版本。或者如果你熟悉numpy，你可能会得到一个使用numpy向量距离和点积的更短的版本。

注意，如果您要处理多边形中的数千个点，则应优化此例程以避免计算边缘之间的所有距离。也许，在计算边缘到边缘的距离时，您可以通过引入一些巧妙的过滤来丢弃边缘。

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