make cut_polygon_with_line work() correctly

src/c3nav/mapdata/utils/geometry.py

@@ -1,14 +1,15 @@
+import math
 from collections import deque, namedtuple
 from itertools import chain
-from typing import List
+from typing import List, Union
 
 import matplotlib.pyplot as plt
 from django.core import checks
 from matplotlib.patches import PathPatch
 from matplotlib.path import Path
 from shapely import speedups
-from shapely.geometry import LinearRing, LineString, Point, Polygon
+from shapely.geometry import GeometryCollection, LinearRing, LineString, MultiLineString, MultiPolygon, Point, Polygon
-from shapely.ops import unary_union
+from shapely.geometry.polygon import orient
 
 if speedups.available:
     speedups.enable()
@@ -41,7 +42,7 @@ def clean_geometry(geometry):
     return geometry
 
 
-def assert_multipolygon(geometry):
+def assert_multipolygon(geometry: Union[Polygon, MultiPolygon, GeometryCollection]) -> List[Polygon]:
     """
     given a Polygon or a MultiPolygon, return a list of Polygons
     :param geometry: a Polygon or a MultiPolygon
@@ -54,7 +55,7 @@ def assert_multipolygon(geometry):
     return [geom for geom in geometry.geoms if isinstance(geom, Polygon)]
 
 
-def assert_multilinestring(geometry):
+def assert_multilinestring(geometry: Union[LineString, MultiLineString, GeometryCollection]) -> List[LineString]:
     """
     given a LineString or MultiLineString, return a list of LineStrings
     :param geometry: a LineString or a MultiLineString
@@ -115,23 +116,29 @@ def get_rings(geometry):
 cutpoint = namedtuple('cutpoint', ('point', 'polygon', 'ring'))
 
 
-def cut_line_with_point(line, point):
+def cut_line_with_point(line: LineString, point: Point):
     distance = line.project(point)
+    if distance <= 0 or distance >= line.length:
+        return line,
     pointlist = [(point.x, point.y)]
+    subdistance = 0
+    last = None
     for i, p in enumerate(line.coords):
-        subdistance = line.project(Point(p))
+        if last is not None:
+            subdistance += ((p[0]-last[0])**2 + (p[1]-last[1])**2)**0.5
+        last = p
         if subdistance >= distance:
             return (LineString(line.coords[:i] + pointlist),
                     LineString(pointlist + line.coords[i+(1 if subdistance == distance else 0):]))
 
 
-def cut_polygon_with_line(polygon, line):
+def cut_polygon_with_line(polygon: Polygon, line: LineString):
-    polygons: List[List[LinearRing]] = [[polygon.exterior, *polygon.interiors]
+    polygons = (orient(polygon) for polygon in assert_multipolygon(polygon))
-                                        for polygon in assert_multipolygon(polygon)]
+    polygons: List[List[LinearRing]] = [[polygon.exterior, *polygon.interiors] for polygon in polygons]
 
     points = deque()
-    for i, polygon in polygons:
+    for i, polygon in enumerate(polygons):
-        for j, ring in polygon:
+        for j, ring in enumerate(polygon):
             intersection = ring.intersection(line)
             for item in getattr(intersection, 'geoms', (intersection, )):
                 if isinstance(item, Point):
@@ -149,15 +156,19 @@ def cut_polygon_with_line(polygon, line):
         current = points.popleft()
         if current.polygon == last.polygon:
             polygon = polygons[current.polygon]
+            segment = cut_line_with_point(cut_line_with_point(line, last.point)[-1], current.point)[0]
             if current.ring != last.ring:
+                # connect rings
                 ring1 = cut_line_with_point(polygon[last.ring], last.point)
                 ring2 = cut_line_with_point(polygon[current.ring], current.point)
-                new_ring = LinearRing(ring1[:-1] + ring1[0:] + ring2[:-1] + ring2[0:] + ring1[:1])
+                new_ring = LinearRing(ring1[1].coords[:-1] + ring1[0].coords[:-1] + segment.coords[:-1] +
+                                      ring2[1].coords[:-1] + ring2[0].coords[:-1] + segment.coords[::-1])
                 if current.ring == 0 or last.ring == 0:
+                    new_i = 0
                     polygon[0] = new_ring
                     interior = current.ring if last.ring == 0 else last.ring
                     polygon[interior] = None
-                    mapping = {interior: 0}
+                    mapping = {interior: new_i}
                 else:
                     new_i = len(polygon)
                     mapping = {last.ring: new_i, current.ring: new_i}
@@ -168,19 +179,41 @@ def cut_polygon_with_line(polygon, line):
                 points = deque((cutpoint(item.point, item.polygon, mapping[item.ring])
                                 if (item.polygon == current.polygon and item.ring in mapping) else item)
                                for item in points)
+                current = cutpoint(current.point, current.polygon, new_i)
+
             elif current.ring == 0:
+                # cut polygon
                 new_i = len(polygons)
                 exterior = cut_line_with_point(polygon[0], current.point)
-                exterior = cut_line_with_point(LinearRing(exterior[:-1] + exterior[0:]), last.point)
+                exterior = cut_line_with_point(LinearRing(exterior[1].coords[:-1] +
-                exterior1 = LinearRing(exterior[0][0:] + exterior[0][:1])
+                                                          exterior[0].coords[0:]), last.point)
-                exterior2 = LinearRing(exterior[1][0:] + exterior[1][:1])
+
+                point_forwards = exterior[1].coords[1]
+                point_backwards = exterior[0].coords[-2]
+                angle_forwards = math.atan2(point_forwards[0] - last.point.x, point_forwards[1] - last.point.y)
+                angle_backwards = math.atan2(point_backwards[0] - last.point.x, point_backwards[1] - last.point.y)
+                angle_segment = math.atan2(current.point.x - last.point.x, current.point.y - last.point.y)
+
+                while angle_forwards <= angle_backwards:
+                    angle_forwards += 2*math.pi
+                if angle_segment < angle_backwards:
+                    while angle_segment < angle_backwards:
+                        angle_segment += 2*math.pi
+                else:
+                    while angle_segment > angle_forwards:
+                        angle_segment -= 2*math.pi
+
+                # cut polygon only if segment is inside polygon
+                if angle_backwards < angle_segment < angle_forwards:
+                    exterior1 = LinearRing(exterior[0].coords[:-1] + segment.coords[0:])
+                    exterior2 = LinearRing(exterior[1].coords[:-1] + segment.coords[::-1])
                     geom = Polygon(exterior1)
                     polygon[0] = exterior1
                     new_polygon = [exterior2]
                     polygons.append(new_polygon)
                     mapping = {}
                     for i, interior in enumerate(polygon[1:]):
-                    if not geom.contains(interior):
+                        if interior is not None and not geom.contains(interior):
                             mapping[i] = len(new_polygon)
                             new_polygon.append(interior)
 
@@ -191,5 +224,7 @@ def cut_polygon_with_line(polygon, line):
                                     if (item.polygon == current.polygon and item.ring == 0 and
                                         not exterior1.contains(item.point)) else item)
                                    for item in points)
-
+                    current = cutpoint(current.point, new_i, 0)
-    return unary_union(tuple(Polygon(polygon[0], polygon[1:]) for polygon in polygons))
+        last = current
+    return tuple(Polygon(polygon[0], tuple(ring for ring in polygon[1:] if ring is not None))
+                 for polygon in polygons)