ENH: Add a function (geometry.split) to split polygons with polygons

CHANGES.md

@@ -4,6 +4,7 @@
 
 - **ENH: Add a function (`arcpy.feature_layer_to_path`) to retrieve the path from an ArcGIS feature layer**
 - **ENH: Add a class `ArcPyLogger` to better handle logs in Arcgis Pro tools**
+- **ENH: Add a function (`geometry.split`) to split polygons with polygons**
 - FIX: Fix `rasters.sieve` with integer data
 
 ## 1.34.2 (2024-01-23)

CI/SCRIPTS/test_geometry.py

@@ -18,7 +18,7 @@
 
 from CI.SCRIPTS.script_utils import geometry_path, s3_env, vectors_path
 from sertit import ci, geometry, vectors
-from sertit.geometry import fill_polygon_holes, get_wider_exterior
+from sertit.geometry import fill_polygon_holes, get_wider_exterior, split
 
 ci.reduce_verbosity()
 
@@ -88,3 +88,25 @@ def test_fill_polygon_holes():
     ci.assert_geom_equal(
         fill_polygon_holes(water, threshold=1000), vectors.read(water_1000_path)
     )
+
+
+@s3_env
+def test_split():
+    """Test fill_polygon_holes"""
+    water_path = geometry_path().joinpath("water.geojson")
+    water = vectors.read(water_path)
+
+    # No MultiLineStrings
+    footprint_path = geometry_path().joinpath("footprint_split.geojson")
+    water_split_path = geometry_path().joinpath("water_split.geojson")
+    ci.assert_geom_equal(
+        split(water, vectors.read(footprint_path)), vectors.read(water_split_path)
+    )
+
+    # With MultiLineStrings
+    footprint_raw_path = geometry_path().joinpath("footprint_raw.geojson")
+    water_split_raw_path = geometry_path().joinpath("water_split_raw.geojson")
+    ci.assert_geom_equal(
+        split(water, vectors.read(footprint_raw_path)),
+        vectors.read(water_split_raw_path),
+    )

sertit/geometry.py

@@ -22,12 +22,14 @@
 import logging
 
 import numpy as np
+from shapely.errors import GeometryTypeError
 from tqdm import tqdm
 
 from sertit.types import AnyPolygonType
 
 try:
     import geopandas as gpd
+    from shapely import ops
     from shapely.geometry import Polygon, box
 except ModuleNotFoundError as ex:
     raise ModuleNotFoundError(
@@ -252,3 +254,31 @@ def _fill_polygon(polygon: Polygon, threshold: float = None):
 
     # Write back to file
     return gpd_results
+
+
+def split(polygons: gpd.GeoDataFrame, splitter: gpd.GeoDataFrame):
+    """
+    Split polygons with polygons
+
+    Args:
+        polygons (gpd.GeoDataFrame): Polygons to split
+        splitter (gpd.GeoDataFrame): Splitter to split the polygons
+
+    Returns:
+        gpd.GeoDataFrame: Splitted GeoDataFrame
+    """
+    out = polygons.geometry
+    for _, split in splitter.iterrows():
+        # Compute the boundary of the splitter polygon (to get a LineString)
+        boundary = split.geometry.boundary
+
+        # Explode to prevent FeatureCollections
+        try:
+            # LineStrings
+            out = out.map(lambda geom: ops.split(geom, boundary).geoms).explode()
+        except GeometryTypeError:
+            # MultiLineStrings
+            for line in boundary.geoms:
+                out = out.map(lambda geom: ops.split(geom, line).geoms).explode()
+
+    return gpd.GeoDataFrame(geometry=out.explode(), crs=polygons.crs)