Simplify edges in enclosed_tessellation

momepy/functional/_elements.py

@@ -15,6 +15,7 @@
 GPD_GE_013 = Version(gpd.__version__) >= Version("0.13.0")
 GPD_GE_10 = Version(gpd.__version__) >= Version("1.0dev")
 LPS_GE_411 = Version(libpysal.__version__) >= Version("4.11.dev")
+SHPLY_GE_210 = Version(shapely.__version__) >= Version("2.1.0")
 
 __all__ = [
     "morphological_tessellation",
@@ -133,6 +134,7 @@
     shrink: float = 0.4,
     segment: float = 0.5,
     threshold: float = 0.05,
+    simplify: bool = True,
     n_jobs: int = -1,
     **kwargs,
 ) -> GeoDataFrame:
@@ -181,6 +183,9 @@
         inlude it in the tessellation of that enclosure. Resolves sliver geometry
         issues. If None, the check is skipped and all intersecting buildings are
         considered. By default 0.05
+    simplify: bool, optional
+        Whether to attempt to simplify the resulting tesselation boundaries with
+        ``shapely.coverage_simplify``. By default True.
     n_jobs : int, optional
         The number of jobs to run in parallel. -1 means using all available cores.
         By default -1
@@ -226,13 +231,19 @@
 
     >>> momepy.enclosed_tessellation(buildings, enclosures).head()
                                                   geometry  enclosure_index
-    0    POLYGON ((1603572.779 6464354.58, 1603572.505 ...                0
-    113  POLYGON ((1603543.601 6464322.376, 1603543.463...                0
-    114  POLYGON ((1603525.157 6464283.592, 1603524.725...                0
-    125  POLYGON ((1603601.446 6464256.455, 1603600.982...                0
-    126  POLYGON ((1603528.593 6464221.033, 1603527.796...                0
+    0    POLYGON ((1603546.697 6464383.596, 1603585.64 ...                0
+    113  POLYGON ((1603517.131 6464349.296, 1603546.697...                0
+    114  POLYGON ((1603517.87 6464285.864, 1603515.152 ...                0
+    125  POLYGON ((1603586.269 6464256.691, 1603581.813...                0
+    126  POLYGON ((1603499.92 6464243.917, 1603493.299 ...                0
     """
 
+    if simplify and not SHPLY_GE_210:
+        raise ImportError(
+            "`simplify=True` requires shapely 2.1 or higher. "
+            "Update shapely or set `simplify` to False."
+        )
+
     if isinstance(geometry.index, MultiIndex):
         raise ValueError(
             "MultiIndex is not supported in `momepy.enclosed_tessellation`."
@@ -274,7 +285,7 @@
 
     # generate tessellation in parallel
     new = Parallel(n_jobs=n_jobs)(
-        delayed(_tess)(*t, threshold, shrink, segment, index_name, kwargs)
+        delayed(_tess)(*t, threshold, shrink, segment, index_name, simplify, kwargs)
         for t in tuples
     )
 
@@ -307,7 +318,7 @@
     return pd.concat([new_df, singles.drop(columns="position"), clean_blocks])
 
 
-def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, kwargs):
+def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, to_simplify, kwargs):
     """Generate tessellation for a single enclosure. Helper for enclosed_tessellation"""
     # check if threshold is set and filter buildings based on the threshold
     if threshold:
@@ -326,6 +337,11 @@
             as_gdf=True,
             **kwargs,
         )
+        if to_simplify:
+            simpl_collection = shapely.coverage_simplify(
+                tess.geometry, tolerance=segment / 2, simplify_boundary=False
+            )
+            tess.geometry = gpd.GeoSeries(simpl_collection).values
         tess[enclosure_id] = ix
         return tess
 

momepy/functional/tests/test_elements.py

@@ -3,6 +3,7 @@
 import numpy as np
 import pandas as pd
 import pytest
+import shapely
 from geopandas.testing import assert_geodataframe_equal
 from packaging.version import Version
 from pandas.testing import assert_index_equal, assert_series_equal
@@ -13,6 +14,7 @@
 
 GPD_GE_013 = Version(gpd.__version__) >= Version("0.13.0")
 LPS_GE_411 = Version(libpysal.__version__) >= Version("4.11.dev")
+SHPLY_GE_210 = Version(shapely.__version__) >= Version("2.1.0")
 
 
 class TestElements:
@@ -94,8 +96,7 @@ def test_morphological_tessellation_errors(self):
 
     def test_enclosed_tessellation(self):
         tessellation = mm.enclosed_tessellation(
-            self.df_buildings,
-            self.enclosures.geometry,
+            self.df_buildings, self.enclosures.geometry, simplify=False
         )
         assert (tessellation.geom_type == "Polygon").all()
         assert tessellation.crs == self.df_buildings.crs
@@ -105,6 +106,7 @@ def test_enclosed_tessellation(self):
         sparser = mm.enclosed_tessellation(
             self.df_buildings,
             self.enclosures.geometry,
+            simplify=False,
             segment=2,
         )
         if GPD_GE_013:
@@ -114,7 +116,11 @@ def test_enclosed_tessellation(self):
             )
 
         no_threshold_check = mm.enclosed_tessellation(
-            self.df_buildings, self.enclosures.geometry, threshold=None, n_jobs=1
+            self.df_buildings,
+            self.enclosures.geometry,
+            simplify=False,
+            threshold=None,
+            n_jobs=1,
         )
 
         assert_geodataframe_equal(tessellation, no_threshold_check)
@@ -135,7 +141,11 @@ def test_enclosed_tessellation(self):
         )
 
         threshold_elimination = mm.enclosed_tessellation(
-            buildings, self.enclosures.geometry, threshold=0.99, n_jobs=1
+            buildings,
+            self.enclosures.geometry,
+            simplify=False,
+            threshold=0.99,
+            n_jobs=1,
         )
         assert not threshold_elimination.index.duplicated().any()
         assert_index_equal(threshold_elimination.index, tessellation.index)
@@ -148,6 +158,7 @@ def test_enclosed_tessellation(self):
         tessellation_df = mm.enclosed_tessellation(
             self.df_buildings,
             self.enclosures,
+            simplify=False,
         )
         assert_geodataframe_equal(tessellation, tessellation_df)
 
@@ -156,6 +167,7 @@ def test_enclosed_tessellation(self):
         tessellation_custom_index = mm.enclosed_tessellation(
             self.df_buildings,
             custom_index,
+            simplify=False,
         )
         assert (tessellation_custom_index.geom_type == "Polygon").all()
         assert tessellation_custom_index.crs == self.df_buildings.crs
@@ -330,6 +342,50 @@ def test_blocks_inner(self):
         else:
             assert len(blocks.sindex.query_bulk(blocks.geometry, "overlaps")[0]) == 0
 
+    @pytest.mark.skipif(not SHPLY_GE_210, reason="coverage_simplify required")
+    def test_simplified_tesselations(self):
+        n_workers = -1
+        tessellations = mm.enclosed_tessellation(
+            self.df_buildings,
+            self.enclosures.geometry,
+            simplify=False,
+            n_jobs=n_workers,
+        )
+        simplified_tessellations = mm.enclosed_tessellation(
+            self.df_buildings, self.enclosures.geometry, simplify=True, n_jobs=n_workers
+        )
+        ## empty enclosures should be unmodified
+        assert_geodataframe_equal(
+            tessellations[tessellations.index < 0],
+            simplified_tessellations[simplified_tessellations.index < 0],
+        )
+        ## simplification should result in less total points
+        orig_points = shapely.get_coordinates(
+            tessellations[tessellations.index >= 0].geometry
+        ).shape
+        simpl_points = shapely.get_coordinates(
+            simplified_tessellations[simplified_tessellations.index >= 0].geometry
+        ).shape
+        assert orig_points > simpl_points
+
+        ## simplification should not modify the external borders of tesselation cells
+        orig_grouper = tessellations.groupby("enclosure_index")
+        simpl_grouper = simplified_tessellations.groupby("enclosure_index")
+        for idx in np.union1d(
+            tessellations["enclosure_index"].unique(),
+            simplified_tessellations["enclosure_index"].unique(),
+        ):
+            orig_group = orig_grouper.get_group(idx).dissolve().boundary
+            enclosure = self.enclosures.loc[[idx]].dissolve().boundary
+
+            simpl_group = simpl_grouper.get_group(idx).dissolve().boundary
+
+            ## simplified is not different to enclosure
+            assert np.isclose(simpl_group.difference(enclosure).area, 0)
+
+            # simplified is not different to original tess
+            assert np.isclose(simpl_group.difference(orig_group).area, 0)
+
     def test_multi_index(self):
         buildings = self.df_buildings.set_index(["uID", "uID"])
         with pytest.raises(
@@ -341,7 +397,7 @@ def test_multi_index(self):
             ValueError,
             match="MultiIndex is not supported in `momepy.enclosed_tessellation`.",
         ):
-            mm.enclosed_tessellation(buildings, self.enclosures)
+            mm.enclosed_tessellation(buildings, self.enclosures, simplify=False)
         with pytest.raises(
             ValueError,
             match="MultiIndex is not supported in `momepy.verify_tessellation`.",
@@ -363,7 +419,7 @@ def test_multi_index(self):
 
     def test_tess_single_building_edge_case(self):
         tessellations = mm.enclosed_tessellation(
-            self.df_buildings, self.enclosures.geometry, n_jobs=-1
+            self.df_buildings, self.enclosures.geometry, simplify=False, n_jobs=-1
         )
         orig_grouper = tessellations.groupby("enclosure_index")
         idxs = ~self.df_buildings.index.isin(orig_grouper.get_group(8).index)
@@ -373,7 +429,9 @@ def test_tess_single_building_edge_case(self):
 
         new_blg = self.df_buildings[idxs]
         new_blg.loc[22, "geometry"] = new_blg.loc[22, "geometry"].buffer(20)
-        new_tess = mm.enclosed_tessellation(new_blg, self.enclosures.geometry, n_jobs=1)
+        new_tess = mm.enclosed_tessellation(
+            new_blg, self.enclosures.geometry, simplify=False, n_jobs=1
+        )
 
         # assert that buildings 1 and 22 intersect the same enclosure
         inp, res = self.enclosures.sindex.query(