Functional blocks

momepy/functional/_elements.py

@@ -10,6 +10,7 @@
 from libpysal.cg import voronoi_frames
 from libpysal.graph import Graph
 from packaging.version import Version
+from pandas import Series
 
 GPD_GE_013 = Version(gpd.__version__) >= Version("0.13.0")
 GPD_GE_10 = Version(gpd.__version__) >= Version("1.0dev")
@@ -20,6 +21,7 @@
     "enclosed_tessellation",
     "verify_tessellation",
     "get_nearest_street",
+    "generate_blocks",
     "buffered_limit",
 ]
 
@@ -354,6 +356,86 @@ def get_nearest_street(
     return ids
 
 
+def generate_blocks(
+    tessellation: GeoDataFrame, edges: GeoDataFrame, buildings: GeoDataFrame
+) -> tuple[Series, Series, Series]:
+    """
+    Generate blocks based on buildings, tessellation, and street network.
+    Dissolves tessellation cells based on street-network based polygons.
+    Links resulting ID to ``buildings`` and ``tessellation`` and returns
+    ``blocks``, ``buildings_ds`` and ``tessellation`` ids.
+
+    Parameters
+    ----------
+    tessellation : GeoDataFrame
+        A GeoDataFrame containing morphological tessellation.
+    edges : GeoDataFrame
+        A GeoDataFrame containing a street network.
+    buildings : GeoDataFrame
+        A GeoDataFrame containing buildings.
+
+    Returns
+    -------
+    blocks : GeoDataFrame
+        A GeoDataFrame containing generated blocks.
+    buildings_ids : Series
+        A Series derived from buildings with block ID.
+    tessellation_ids : Series
+        A Series derived from morphological tessellation with block ID.
+
+    Examples
+    --------
+    >>> blocks, buildings_id, tessellation_id = mm.generate_blocks(tessellation_df,
+    ... streets_df, buildings_df)
+    >>> blocks.head()
+            geometry
+    0	    POLYGON ((1603560.078648818 6464202.366899694,...
+    1	    POLYGON ((1603457.225976106 6464299.454696888,...
+    2	    POLYGON ((1603056.595487018 6464093.903488506,...
+    3	    POLYGON ((1603260.943782872 6464141.327631323,...
+    4	    POLYGON ((1603183.399594798 6463966.109982309,...
+    """
+
+    id_name: str = "bID"
+
+    # slice the tessellations by the street network
+    cut = gpd.overlay(
+        tessellation,
+        gpd.GeoDataFrame(geometry=edges.buffer(0.001)),
+        how="difference",
+    )
+    cut = cut.explode(ignore_index=True)
+    # touching tessellations form a block
+    weights = Graph.build_contiguity(cut, rook=False)
+    cut["component"] = weights.component_labels
+
+    # generate block geometries
+    buildings_c = buildings.copy()
+    buildings_c.geometry = buildings_c.representative_point()  # make points
+    centroids_temp_id = gpd.sjoin(
+        buildings_c,
+        cut[[cut.geometry.name, "component"]],
+        how="left",
+        predicate="within",
+    )
+    cells_copy = tessellation[[tessellation.geometry.name]].merge(
+        centroids_temp_id[["component"]], right_index=True, left_index=True, how="left"
+    )
+    blocks = cells_copy.dissolve(by="component").explode(ignore_index=True)
+
+    # assign block ids to buildings and tessellations
+    centroids_w_bl_id2 = gpd.sjoin(buildings_c, blocks, how="left", predicate="within")
+    # return blocks, centroids_w_bl_id2
+    buildings_id = centroids_w_bl_id2["index_right"]
+    buildings_id.name = id_name
+    cells_m = tessellation.merge(
+        buildings_id, left_index=True, right_index=True, how="left"
+    )
+    tessellation_id = cells_m[id_name]
+
+    return blocks, buildings_id, tessellation_id
+
+
 def buffered_limit(
     gdf: GeoDataFrame | GeoSeries,
     buffer: float | str = 100,

momepy/functional/tests/test_elements.py

@@ -5,7 +5,7 @@
 import pytest
 from geopandas.testing import assert_geodataframe_equal
 from packaging.version import Version
-from pandas.testing import assert_index_equal
+from pandas.testing import assert_index_equal, assert_series_equal
 from shapely import affinity
 from shapely.geometry import MultiPoint, Polygon, box
 
@@ -256,3 +256,57 @@ def test_buffered_limit_error(self):
             ValueError, match="`buffer` must be either 'adaptive' or a number."
         ):
             mm.buffered_limit(self.df_buildings, "invalid")
+
+    def test_blocks(self):
+        blocks, buildings_id, tessellation_id = mm.generate_blocks(
+            self.df_tessellation, self.df_streets, self.df_buildings
+        )
+        assert not tessellation_id.isna().any()
+        assert not buildings_id.isna().any()
+        assert len(blocks) == 8
+
+    def test_blocks_inner(self):
+        streets = self.df_streets.copy()
+        streets.loc[35, "geometry"] = (
+            self.df_buildings.geometry.iloc[141]
+            .representative_point()
+            .buffer(20)
+            .exterior
+        )
+        blocks, buildings_id, tessellation_id = mm.generate_blocks(
+            self.df_tessellation, streets, self.df_buildings
+        )
+        assert not tessellation_id.isna().any()
+        assert not buildings_id.isna().any()
+        assert len(blocks) == 9
+        if GPD_GE_013:
+            assert len(blocks.sindex.query(blocks.geometry, "overlaps")[0]) == 0
+        else:
+            assert len(blocks.sindex.query_bulk(blocks.geometry, "overlaps")[0]) == 0
+
+
+class TestElementsEquivalence:
+    def setup_method(self):
+        test_file_path = mm.datasets.get_path("bubenec")
+        self.df_buildings = gpd.read_file(test_file_path, layer="buildings")
+        self.df_tessellation = gpd.read_file(test_file_path, layer="tessellation")
+        self.df_streets = gpd.read_file(test_file_path, layer="streets")
+        self.limit = mm.buffered_limit(self.df_buildings, 50)
+        self.enclosures = mm.enclosures(
+            self.df_streets,
+            gpd.GeoSeries([self.limit.exterior], crs=self.df_streets.crs),
+        )
+
+    def test_blocks(self):
+        blocks, buildings_id, tessellation_id = mm.generate_blocks(
+            self.df_tessellation, self.df_streets, self.df_buildings
+        )
+        res = mm.Blocks(
+            self.df_tessellation, self.df_streets, self.df_buildings, "bID", "uID"
+        )
+
+        assert_geodataframe_equal(
+            blocks.geometry.to_frame(), res.blocks.geometry.to_frame()
+        )
+        assert_series_equal(buildings_id, res.buildings_id)
+        assert_series_equal(tessellation_id, res.tessellation_id)