Finden Sie den kürzesten Weg von einer Geometrie zu anderen auf shapely

Question

Wie zwei Linien line_a und line_b, wie finde ich das Paar der Punkte, die die kleinere Route von line_a zu line_b darstellt?

Answer 1

Leider gibt es keine formschöne Operation, um dies zu tun. Ich dachte über dieses Problem nach, als ich gebeten wurde, die Lösung in der Frage find-coordinate-of-closest-point-on-polygon-shapely zu erweitern, um zwei Polygone zu behandeln. Die Antwort hängt vom Geometriesatz ab, der intuitiv wahr ist (obwohl ein formaler Beweis einen Kalkül erfordert und es ist ein bisschen zu lang, hier ohne LaTex geschrieben zu werden). Das Theorem sagt:

Der Mindestabstand zwischen zwei Zeilenzeichenfolgen, die (im Falle einer Polygonzugs eine verkettete Sequenz von Segmenten) nicht gegen kreuzen, wird immer in einem der Randpunkte der Linie erreicht, Saiten.

dies mit im Sinn, ist das Problem der Mindestabstand zwischen den einzelnen Kante eines Linestring zu dem anderen Linestring zu berechnen, verringert.

Eine weitere Möglichkeit, dieses Problem zu sehen, besteht darin, den Abstand zwischen jedem Segmentpaar auf zu reduzieren. Und zu bemerken, dass der Abstand zwischen zwei Segmenten, die sich nicht schneiden, zwischen zwei Endpunkten, oder zwischen einem Endpunkt und der Projektion dieses Endpunktes über das andere Segment erreicht wird.

Der Code wurde ein wenig optimiert, um redundante Berechnungen zu vermeiden, aber vielleicht können Sie eine elegantere Version erhalten. Oder wenn Sie mit numpy vertraut sind, können Sie wahrscheinlich eine kürzere Version erhalten, die numpy Vektorabstand und Punktprodukt verwendet.

Vorsicht, wenn Sie mehrere tausend Punkte in den Polygonen behandeln, sollte diese Routine optimiert werden, um die Berechnung aller Abstände zwischen Kanten zu vermeiden. Vielleicht könnten Sie, während Sie die Entfernung von Kante zu Kante berechnen, die Kante ablegen, indem Sie eine clevere Filterung einführen.

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