diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py index f783cb17..88df5d67 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_210 = Version(shapely.__version__) >= Version("2.1.0") __all__ = [ "morphological_tessellation", @@ -133,6 +134,7 @@ def enclosed_tessellation( 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 @@ 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 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 @@ def enclosed_tessellation( >>> 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 @@ def enclosed_tessellation( # 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 @@ def enclosed_tessellation( 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 @@ def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, kwargs): 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 diff --git a/momepy/functional/tests/test_elements.py b/momepy/functional/tests/test_elements.py index 9a5c1425..00f432b7 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.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(