From b0235548aa716242228713331100921a32be4838 Mon Sep 17 00:00:00 2001
From: Krasen Samardzhiev <krasensam@gmail.com>
Date: Sat, 25 May 2024 23:30:33 +0200
Subject: [PATCH 1/9] generate tesselations for individual buildings

---
 momepy/functional/_elements.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py
index 97cfe585..76d067de 100644
--- a/momepy/functional/_elements.py
+++ b/momepy/functional/_elements.py
@@ -241,7 +241,7 @@ def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id):
             > (shapely.area(blg.geometry.array) * threshold)
         ]
 
-    if len(blg) > 1:
+    if len(blg) >= 1:
         tess = voronoi_frames(
             blg,
             clip=poly,

From 965895ede987ee59c23dd957d9200e3e979c7f11 Mon Sep 17 00:00:00 2001
From: Krasen Samardzhiev <krasensam@gmail.com>
Date: Wed, 12 Jun 2024 15:00:51 +0200
Subject: [PATCH 2/9] simplify code

---
 momepy/functional/_elements.py | 27 +++++++++++++++++++++++++--
 1 file changed, 25 insertions(+), 2 deletions(-)

diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py
index 43b19dae..13088503 100644
--- a/momepy/functional/_elements.py
+++ b/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_205 = Version(shapely.__version__) >= Version("2.0.5")
 
 __all__ = [
     "morphological_tessellation",
@@ -101,6 +102,7 @@ def enclosed_tessellation(
     shrink: float = 0.4,
     segment: float = 0.5,
     threshold: float = 0.05,
+    simplify: bool = False,
     n_jobs: int = -1,
 ) -> GeoDataFrame:
     """Generate enclosed tessellation
@@ -148,6 +150,9 @@ def enclosed_tessellation(
         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 False.
     n_jobs : int, optional
         The number of jobs to run in parallel. -1 means using all available cores.
         By default -1
@@ -190,6 +195,18 @@ def enclosed_tessellation(
     single = unique[counts == 1]
     altered = unique[counts > 0]
 
+    ##simplify enclosures now that the spatial queries are
+    if simplify:
+        if not SHPLY_GE_205:
+            raise ImportError(
+                "Coverage simplification requires shapely 2.0.5 or higher."
+            )
+        from shapely import coverage_simplify
+
+        to_simplify = np.union1d(single, splits)
+        simplified = coverage_simplify(enclosures.iloc[to_simplify], 1e-1)
+        enclosures.iloc[to_simplify, 0] = simplified
+
     # prepare input for parallel processing
     tuples = [
         (
@@ -202,7 +219,8 @@ def enclosed_tessellation(
 
     # generate tessellation in parallel
     new = Parallel(n_jobs=n_jobs)(
-        delayed(_tess)(*t, threshold, shrink, segment, index_name) for t in tuples
+        delayed(_tess)(*t, threshold, shrink, segment, index_name, to_simplify)
+        for t in tuples
     )
 
     new_df = pd.concat(new, axis=0)
@@ -234,7 +252,7 @@ def enclosed_tessellation(
     return pd.concat([new_df, singles.drop(columns="position"), clean_blocks])
 
 
-def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id):
+def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, to_simplify):
     """Generate tessellation for a single enclosure. Helper for enclosed_tessellation"""
     # check if threshold is set and filter buildings based on the threshold
     if threshold:
@@ -252,6 +270,11 @@ def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id):
             return_input=False,
             as_gdf=True,
         )
+        if to_simplify:
+            from shapely import coverage_simplify
+
+            tess.geometry = coverage_simplify(tess.geometry, 1e-1)
+
         tess[enclosure_id] = ix
         return tess
 

From f92acd37d5bcde2e6745d4dccf79a5634884092e Mon Sep 17 00:00:00 2001
From: Krasen Samardzhiev <krasensam@gmail.com>
Date: Wed, 12 Jun 2024 15:04:55 +0200
Subject: [PATCH 3/9] simplify code

---
 momepy/functional/_elements.py | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py
index 13088503..1f746122 100644
--- a/momepy/functional/_elements.py
+++ b/momepy/functional/_elements.py
@@ -195,7 +195,6 @@ def enclosed_tessellation(
     single = unique[counts == 1]
     altered = unique[counts > 0]
 
-    ##simplify enclosures now that the spatial queries are
     if simplify:
         if not SHPLY_GE_205:
             raise ImportError(
@@ -203,9 +202,9 @@ def enclosed_tessellation(
             )
         from shapely import coverage_simplify
 
-        to_simplify = np.union1d(single, splits)
-        simplified = coverage_simplify(enclosures.iloc[to_simplify], 1e-1)
-        enclosures.iloc[to_simplify, 0] = simplified
+        ##simplify singles here and let the threads simplify the actual tesselations
+        simplified = coverage_simplify(enclosures.iloc[single], 1e-1)
+        enclosures.iloc[single, 0] = simplified
 
     # prepare input for parallel processing
     tuples = [
@@ -219,7 +218,7 @@ def enclosed_tessellation(
 
     # generate tessellation in parallel
     new = Parallel(n_jobs=n_jobs)(
-        delayed(_tess)(*t, threshold, shrink, segment, index_name, to_simplify)
+        delayed(_tess)(*t, threshold, shrink, segment, index_name, simplify)
         for t in tuples
     )
 
@@ -260,6 +259,8 @@ def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, to_simplify):
             shapely.area(shapely.intersection(blg.geometry.array, poly))
             > (shapely.area(blg.geometry.array) * threshold)
         ]
+    if to_simplify:
+        from shapely import coverage_simplify
 
     if len(blg) >= 1:
         tess = voronoi_frames(
@@ -271,13 +272,13 @@ def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, to_simplify):
             as_gdf=True,
         )
         if to_simplify:
-            from shapely import coverage_simplify
-
             tess.geometry = coverage_simplify(tess.geometry, 1e-1)
-
         tess[enclosure_id] = ix
         return tess
 
+    if to_simplify:
+        poly = coverage_simplify(poly, 1e-1)
+
     return GeoDataFrame(
         {enclosure_id: ix},
         geometry=[poly],

From a383cc09120281dac6d79c3a93d9ad96e8dd95ac Mon Sep 17 00:00:00 2001
From: Krasen Samardzhiev <krasensam@gmail.com>
Date: Wed, 12 Jun 2024 15:41:27 +0200
Subject: [PATCH 4/9] PR changes

---
 momepy/functional/_elements.py | 27 +++++++++------------------
 1 file changed, 9 insertions(+), 18 deletions(-)

diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py
index 1f746122..41dc6d94 100644
--- a/momepy/functional/_elements.py
+++ b/momepy/functional/_elements.py
@@ -15,7 +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_205 = Version(shapely.__version__) >= Version("2.0.5")
+SHPLY_GE_210 = Version(shapely.__version__) >= Version("2.1")
 
 __all__ = [
     "morphological_tessellation",
@@ -171,6 +171,9 @@ def enclosed_tessellation(
     momepy.verify_tessellation
     """
 
+    if simplify and not SHPLY_GE_210:
+        raise ImportError("Coverage simplification requires shapely 2.0.5 or higher.")
+
     # convert to GeoDataFrame and add position (we will need it later)
     enclosures = enclosures.geometry.to_frame()
     enclosures["position"] = range(len(enclosures))
@@ -195,17 +198,6 @@ def enclosed_tessellation(
     single = unique[counts == 1]
     altered = unique[counts > 0]
 
-    if simplify:
-        if not SHPLY_GE_205:
-            raise ImportError(
-                "Coverage simplification requires shapely 2.0.5 or higher."
-            )
-        from shapely import coverage_simplify
-
-        ##simplify singles here and let the threads simplify the actual tesselations
-        simplified = coverage_simplify(enclosures.iloc[single], 1e-1)
-        enclosures.iloc[single, 0] = simplified
-
     # prepare input for parallel processing
     tuples = [
         (
@@ -259,8 +251,6 @@ def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, to_simplify):
             shapely.area(shapely.intersection(blg.geometry.array, poly))
             > (shapely.area(blg.geometry.array) * threshold)
         ]
-    if to_simplify:
-        from shapely import coverage_simplify
 
     if len(blg) >= 1:
         tess = voronoi_frames(
@@ -272,13 +262,14 @@ def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, to_simplify):
             as_gdf=True,
         )
         if to_simplify:
-            tess.geometry = coverage_simplify(tess.geometry, 1e-1)
+            from shapely import coverage_simplify
+
+            tess.geometry = coverage_simplify(
+                tess.geometry, tolerance=1e-1, simplify_boundary=False
+            )
         tess[enclosure_id] = ix
         return tess
 
-    if to_simplify:
-        poly = coverage_simplify(poly, 1e-1)
-
     return GeoDataFrame(
         {enclosure_id: ix},
         geometry=[poly],

From 8199333bd29841ce95748e67858d52867fd789d7 Mon Sep 17 00:00:00 2001
From: Krasen Samardzhiev <krasensam@gmail.com>
Date: Wed, 12 Jun 2024 15:51:12 +0200
Subject: [PATCH 5/9] changed assignment

---
 momepy/functional/_elements.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py
index 41dc6d94..44b50276 100644
--- a/momepy/functional/_elements.py
+++ b/momepy/functional/_elements.py
@@ -252,7 +252,7 @@ def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, to_simplify):
             > (shapely.area(blg.geometry.array) * threshold)
         ]
 
-    if len(blg) >= 1:
+    if len(blg) > 1:
         tess = voronoi_frames(
             blg,
             clip=poly,
@@ -270,10 +270,12 @@ def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, to_simplify):
         tess[enclosure_id] = ix
         return tess
 
+    assigned_ix = ix if len(blg) == 1 else -1
+
     return GeoDataFrame(
         {enclosure_id: ix},
         geometry=[poly],
-        index=[-1],
+        index=[assigned_ix],
         crs=blg.crs,
     )
 

From 653a8c6f646d8f698f8d38b5617b4f8143c0ae87 Mon Sep 17 00:00:00 2001
From: Krasen Samardzhiev <krasensam@gmail.com>
Date: Wed, 12 Jun 2024 15:51:23 +0200
Subject: [PATCH 6/9] changed assignment

---
 momepy/functional/_elements.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py
index 44b50276..0909eca5 100644
--- a/momepy/functional/_elements.py
+++ b/momepy/functional/_elements.py
@@ -270,6 +270,7 @@ def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, to_simplify):
         tess[enclosure_id] = ix
         return tess
 
+    ## in case a single building is left in blg
     assigned_ix = ix if len(blg) == 1 else -1
 
     return GeoDataFrame(

From 7afc8e57960d5dfa025e3c5c884ab2d8475a16c4 Mon Sep 17 00:00:00 2001
From: Krasen Samardzhiev <krasensam@gmail.com>
Date: Thu, 13 Jun 2024 17:34:11 +0200
Subject: [PATCH 7/9] simplification & single building catch

---
 momepy/functional/_elements.py           | 12 +++--
 momepy/functional/tests/test_elements.py | 68 ++++++++++++++++++++++++
 2 files changed, 75 insertions(+), 5 deletions(-)

diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py
index 0909eca5..48dd4028 100644
--- a/momepy/functional/_elements.py
+++ b/momepy/functional/_elements.py
@@ -15,7 +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")
+SHPLY_GE_210 = Version(shapely.__version__) >= Version("2.1.0dev")
 
 __all__ = [
     "morphological_tessellation",
@@ -262,16 +262,18 @@ def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, to_simplify):
             as_gdf=True,
         )
         if to_simplify:
-            from shapely import coverage_simplify
-
-            tess.geometry = coverage_simplify(
+            simpl_collection = shapely.coverage_simplify(
                 tess.geometry, tolerance=1e-1, simplify_boundary=False
             )
+            tess.geometry = gpd.GeoSeries(simpl_collection.geoms).values
         tess[enclosure_id] = ix
         return tess
 
     ## in case a single building is left in blg
-    assigned_ix = ix if len(blg) == 1 else -1
+    if len(blg) == 1:
+        assigned_ix = blg.index[0]
+    else:
+        assigned_ix = -1
 
     return GeoDataFrame(
         {enclosure_id: ix},
diff --git a/momepy/functional/tests/test_elements.py b/momepy/functional/tests/test_elements.py
index c63d6002..1109cf8a 100644
--- a/momepy/functional/tests/test_elements.py
+++ b/momepy/functional/tests/test_elements.py
@@ -284,6 +284,74 @@ def test_blocks_inner(self):
         else:
             assert len(blocks.sindex.query_bulk(blocks.geometry, "overlaps")[0]) == 0
 
+    def test_simplified_tesselations(self):
+        n_workers = -1
+        tessellations = mm.enclosed_tessellation(
+            self.df_buildings, self.enclosures.geometry, 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],
+        )
+
+        import shapely
+
+        ## 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
+            ## this needs to be redone
+            assert np.isclose(simpl_group.difference(enclosure).area, 0)
+
+            # simplified is not different to original tess
+            ## this needs to be redone
+            assert np.isclose(simpl_group.difference(orig_group).area, 0)
+
+    def test_tess_single_building_edge_case(self):
+        tessellations = mm.enclosed_tessellation(
+            self.df_buildings, self.enclosures.geometry, n_jobs=-1
+        )
+        orig_grouper = tessellations.groupby("enclosure_index")
+        idxs = ~self.df_buildings.index.isin(orig_grouper.get_group(8).index)
+        idxs[1] = True
+        idxs[21] = False
+        idxs[23] = False
+
+        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)
+
+        ##assert that buildings 1 and 22 intersect the same enclosure
+        inp, res = self.enclosures.sindex.query(
+            new_blg.geometry, predicate="intersects"
+        )
+        assert np.isclose(new_blg.iloc[inp[res == 8]].index.values, [1, 22]).all()
+
+        ### assert that there is a tessellation for building 1
+        assert 1 in new_tess.index
+
 
 class TestElementsEquivalence:
     def setup_method(self):

From c5a2ee984aeb2325d80b34f7eaa723b3663c8e51 Mon Sep 17 00:00:00 2001
From: Krasen Samardzhiev <krasensam@gmail.com>
Date: Thu, 13 Jun 2024 17:41:05 +0200
Subject: [PATCH 8/9] test versionning

---
 momepy/functional/tests/test_elements.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/momepy/functional/tests/test_elements.py b/momepy/functional/tests/test_elements.py
index 1109cf8a..7c0f93e3 100644
--- a/momepy/functional/tests/test_elements.py
+++ b/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.0dev")
 
 
 class TestElements:
@@ -284,6 +286,7 @@ 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(
@@ -297,9 +300,6 @@ def test_simplified_tesselations(self):
             tessellations[tessellations.index < 0],
             simplified_tessellations[simplified_tessellations.index < 0],
         )
-
-        import shapely
-
         ## simplification should result in less total points
         orig_points = shapely.get_coordinates(
             tessellations[tessellations.index >= 0].geometry

From 1891b01e91c324704814b2649d300cf81c96771c Mon Sep 17 00:00:00 2001
From: Krasen Samardzhiev <krasensam@gmail.com>
Date: Thu, 13 Jun 2024 20:04:45 +0200
Subject: [PATCH 9/9] test versionning

---
 momepy/functional/_elements.py           | 6 +++---
 momepy/functional/tests/test_elements.py | 4 ++--
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py
index 48dd4028..997fa33a 100644
--- a/momepy/functional/_elements.py
+++ b/momepy/functional/_elements.py
@@ -15,7 +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.0dev")
+SHPLY_GE_250 = Version(shapely.__version__) >= Version("2.5.0dev")
 
 __all__ = [
     "morphological_tessellation",
@@ -171,8 +171,8 @@ def enclosed_tessellation(
     momepy.verify_tessellation
     """
 
-    if simplify and not SHPLY_GE_210:
-        raise ImportError("Coverage simplification requires shapely 2.0.5 or higher.")
+    if simplify and not SHPLY_GE_250:
+        raise ImportError("Coverage simplification requires shapely 2.5 or higher.")
 
     # convert to GeoDataFrame and add position (we will need it later)
     enclosures = enclosures.geometry.to_frame()
diff --git a/momepy/functional/tests/test_elements.py b/momepy/functional/tests/test_elements.py
index 7c0f93e3..e330d62a 100644
--- a/momepy/functional/tests/test_elements.py
+++ b/momepy/functional/tests/test_elements.py
@@ -14,7 +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.0dev")
+SHPLY_GE_250 = Version(shapely.__version__) >= Version("2.5.0dev")
 
 
 class TestElements:
@@ -286,7 +286,7 @@ 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")
+    @pytest.mark.skipif(not SHPLY_GE_250, reason="coverage_simplify required")
     def test_simplified_tesselations(self):
         n_workers = -1
         tessellations = mm.enclosed_tessellation(