Grid.split() is now cutline order independent

threedigrid_builder/grid/grid.py

@@ -7,7 +7,7 @@
 from osgeo import gdal, ogr, osr
 from pyproj import CRS
 from pyproj.exceptions import CRSError
-from shapely.ops import linemerge, split
+from shapely.ops import linemerge, polygonize
 
 import threedigrid_builder
 from threedigrid_builder.base import (
@@ -1006,32 +1006,27 @@ def apply_cutlines(self, cutlines, dem_path):
     def split(
         polygon: shapely.Polygon, cutlines: List[shapely.LineString]
     ) -> List[shapely.Polygon]:
-        """Sequentially applies cutlines to the polygon"""
-        fragment_collection = [
-            polygon
-        ]  # List of all fragments of the cell (each cutline can cause new fragments), initially the cell itself
-
-        for cutline_linestring in cutlines:
-            new_fragments = []  # List of newly cut fragments (might include original)
-            # Iterate over the fragments
-            for fragment in fragment_collection:
-                intermediate_fragments = split(fragment, cutline_linestring)
-                # According to docs: If the splitter does not split the geometry, a collection with a single geometry equal to the input geometry is returned.
-                # Note that the original input geometry does not need to be returned, we'll explicitely add that one again.
-                if len(intermediate_fragments.geoms) == 1:
-                    assert intermediate_fragments.geoms[0].equals(fragment)
-                    new_fragments.append(fragment)
-                elif len(intermediate_fragments.geoms) > 1:
-                    for intermediate_fragment in intermediate_fragments.geoms:
-                        assert intermediate_fragment.geom_type == "Polygon"
-                        new_fragments.append(intermediate_fragment)
-                else:
-                    raise RuntimeError("No resulting geometry after cutting")
+        """Integral cutting, derived from shapely.ops.cut"""
+
+        union = shapely.unary_union([polygon.boundary] + cutlines)
+
+        # Some polygonized geometries may be holes, we do not want them
+        # # that's why we test if the original polygon (poly) contains
+        # an inner point of polygonized geometry (pg)
+        result = [
+            pg
+            for pg in polygonize(union)
+            if polygon.contains(pg.representative_point())
+        ]
 
-            # Set new fragment collection and move to the next cutline
-            fragment_collection = new_fragments
+        # In case there is no cut, the resulting polygon is an equal shape, but not
+        # numerical the same (different order, additional points). In this case we
+        # explicitely return the original polygon
+        if len(result) == 1:
+            assert result[0].equals(polygon)
+            return [polygon]
 
-        return fragment_collection
+        return result
 
     def finalize(self):
         """Finalize the Grid, computing and setting derived attributes"""

threedigrid_builder/tests/test_grid.py

@@ -389,13 +389,13 @@ def test_split_partial_cut(cell_polygon):
 
 
 def test_split_order(cell_polygon):
-    # This test is to indicate that the order matters
+    # This test is to indicate that the cutline order does not matter
     cutlines_2 = [
-        shapely.LineString([(5, 0), (5, 6)]),
+        shapely.LineString([(5, 0), (5, 5)]),
         shapely.LineString([(0, 0), (10, 10)]),
     ]
     fragments = Grid.split(cell_polygon, cutlines_2)
-    assert len(fragments) == 2
+    assert len(fragments) == 3
 
     cutlines_3 = [
         shapely.LineString([(0, 0), (10, 10)]),
@@ -404,11 +404,6 @@ def test_split_order(cell_polygon):
     fragments = Grid.split(cell_polygon, cutlines_3)
     assert len(fragments) == 3
 
-    # Combining the lines in a single linestring solves the problem
-    cutlines_3 = [shapely.LineString([(0, 0), (10, 10), (5, 5), (5, 0)])]
-    fragments = Grid.split(cell_polygon, cutlines_3)
-    assert len(fragments) == 3
-
 
 def test_split_double_cut(cell_polygon):
     cutlines = [
@@ -417,18 +412,41 @@ def test_split_double_cut(cell_polygon):
     ]
     fragments = Grid.split(cell_polygon, cutlines)
     assert len(fragments) == 4
-    assert shapely.Polygon(((10, 10), (10, 0), (5, 5), (10, 10))).equals(
-        fragments[0]
-    )  # right triangle
-    assert shapely.Polygon(((10, 0), (0, 0), (5, 5), (10, 0))).equals(
-        fragments[1]
-    )  # bottom triangle
-    assert shapely.Polygon(((0, 0), (0, 10), (5, 5), (0, 0))).equals(
-        fragments[2]
-    )  # left triangle
-    assert shapely.Polygon(((0, 10), (10, 10), (5, 5), (0, 10))).equals(
-        fragments[3]
-    )  # top triangle
+    assert shapely.Polygon(((10, 10), (10, 0), (5, 5), (10, 10))).equals(fragments[0])
+    assert shapely.Polygon(((10, 0), (0, 0), (5, 5), (10, 0))).equals(fragments[3])
+    assert shapely.Polygon(((0, 0), (0, 10), (5, 5), (0, 0))).equals(fragments[2])
+    assert shapely.Polygon(((0, 10), (10, 10), (5, 5), (0, 10))).equals(fragments[1])
+
+
+def test_split_double_cut_outer_intersection(cell_polygon):
+    # these cutlines meet outside the clone cell, potentially yielding an
+    # additional cell in polygonize()
+
+    cutlines = [
+        shapely.LineString([(0, 0), (5, 15)]),
+        shapely.LineString([(5, 15), (10, 0)]),
+    ]
+    fragments = Grid.split(cell_polygon, cutlines)
+    assert len(fragments) == 3
+    assert shapely.Polygon(
+        ((6.666666666666667, 10), (10, 10), (10, 0), (6.666666666666667, 10))
+    ).equals(
+        fragments[0]  # right triangle
+    )
+    assert shapely.Polygon(
+        (
+            (3.3333333333333335, 10),
+            (6.666666666666667, 10),
+            (10, 0),
+            (0, 0),
+            (3.3333333333333335, 10),
+        )
+    ).equals(
+        fragments[1]  # trapezium
+    )
+    assert shapely.Polygon(((0, 0), (0, 10), (3.3333333333333335, 10), (0, 0))).equals(
+        fragments[2]  # left triangle
+    )
 
 
 def test_split_double_cut_miss(cell_polygon):
@@ -440,18 +458,10 @@ def test_split_double_cut_miss(cell_polygon):
     ]
     fragments = Grid.split(cell_polygon, cutlines)
     assert len(fragments) == 4
-    assert shapely.Polygon(((10, 10), (10, 0), (5, 5), (10, 10))).equals(
-        fragments[0]
-    )  # right triangle
-    assert shapely.Polygon(((10, 0), (0, 0), (5, 5), (10, 0))).equals(
-        fragments[1]
-    )  # bottom triangle
-    assert shapely.Polygon(((0, 0), (0, 10), (5, 5), (0, 0))).equals(
-        fragments[2]
-    )  # left triangle
-    assert shapely.Polygon(((0, 10), (10, 10), (5, 5), (0, 10))).equals(
-        fragments[3]
-    )  # top triangle
+    assert shapely.Polygon(((10, 10), (10, 0), (5, 5), (10, 10))).equals(fragments[0])
+    assert shapely.Polygon(((10, 0), (0, 0), (5, 5), (10, 0))).equals(fragments[3])
+    assert shapely.Polygon(((0, 0), (0, 10), (5, 5), (0, 0))).equals(fragments[2])
+    assert shapely.Polygon(((0, 10), (10, 10), (5, 5), (0, 10))).equals(fragments[1])
 
 
 def test_split_triple_cut_miss(cell_polygon):
@@ -463,18 +473,17 @@ def test_split_triple_cut_miss(cell_polygon):
     ]
     fragments = Grid.split(cell_polygon, cutlines)
     assert len(fragments) == 5
-    assert shapely.Polygon(((10, 10), (10, 0), (5, 5), (10, 10))).equals(
-        fragments[0]
-    )  # right triangle
-    assert shapely.Polygon(((10, 0), (5, 0), (5, 5), (10, 0))).equals(
-        fragments[1]
-    )  # right half of bottom triangle
-    assert shapely.Polygon(((5, 0), (0, 0), (5, 5), (5, 0))).equals(
-        fragments[2]
-    )  # left half of bottom triangle
-    assert shapely.Polygon(((0, 0), (0, 10), (5, 5), (0, 0))).equals(
-        fragments[3]
-    )  # left triangle
-    assert shapely.Polygon(((0, 10), (10, 10), (5, 5), (0, 10))).equals(
-        fragments[4]
-    )  # top triangle
+    assert shapely.Polygon(((10, 10), (10, 0), (5, 5), (10, 10))).equals(fragments[0])
+    assert shapely.Polygon(((10, 0), (5, 0), (5, 5), (10, 0))).equals(fragments[4])
+    assert shapely.Polygon(((5, 0), (0, 0), (5, 5), (5, 0))).equals(fragments[3])
+    assert shapely.Polygon(((0, 0), (0, 10), (5, 5), (0, 0))).equals(fragments[2])
+    assert shapely.Polygon(((0, 10), (10, 10), (5, 5), (0, 10))).equals(fragments[1])
+
+    cutlines.reverse()
+    fragments = Grid.split(cell_polygon, cutlines)
+    assert len(fragments) == 5
+    assert shapely.Polygon(((10, 10), (10, 0), (5, 5), (10, 10))).equals(fragments[0])
+    assert shapely.Polygon(((10, 0), (5, 0), (5, 5), (10, 0))).equals(fragments[4])
+    assert shapely.Polygon(((5, 0), (0, 0), (5, 5), (5, 0))).equals(fragments[3])
+    assert shapely.Polygon(((0, 0), (0, 10), (5, 5), (0, 0))).equals(fragments[2])
+    assert shapely.Polygon(((0, 10), (10, 10), (5, 5), (0, 10))).equals(fragments[1])