ENH: mean_interbuilding_distance and building_adjacency (#556)

* mean_interbuilding_distance

* building_adjacency

* refactor it all

* cleanup

* bump min libpysal to 4.10

* assing self weights only when needed

* clean a path that will never happen

* properly bump libpysal pin

---------

Co-authored-by: James Gaboardi <[email protected]>

ci/envs/310-latest.yaml

@@ -5,7 +5,7 @@ dependencies:
   - python=3.10
   - geopandas
   - inequality
-  - libpysal>=4.8.0
+  - libpysal>=4.10.0
   - mapclassify
   - networkx
   - packaging

ci/envs/310-oldest.yaml

@@ -5,7 +5,7 @@ dependencies:
   - python=3.10
   - geopandas=0.12
   - inequality
-  - libpysal=4.8.0
+  - libpysal=4.10.0
   - mapclassify
   - networkx=2.7
   - numpy=1.22

ci/envs/311-latest.yaml

@@ -5,7 +5,7 @@ dependencies:
   - python=3.11
   - geopandas
   - inequality
-  - libpysal>=4.8.0
+  - libpysal>=4.10.0
   - mapclassify
   - networkx
   - packaging

ci/envs/312-dev.yaml

@@ -6,7 +6,7 @@ dependencies:
   - dask
   - geopandas
   - inequality
-  - libpysal>=4.8.0
+  - libpysal>=4.10.0
   - networkx
   - osmnx
   - packaging

ci/envs/312-latest.yaml

@@ -5,7 +5,7 @@ dependencies:
   - python=3.12
   - geopandas
   - inequality
-  - libpysal>=4.8.0
+  - libpysal>=4.10.0
   - mapclassify
   - networkx
   - osmnx

docs/install.rst

@@ -109,7 +109,7 @@ Dependencies
 Required dependencies:
 
 - `geopandas`_ (>= 0.12.0)
-- `libpysal`_ (>= 4.8.0)
+- `libpysal`_ (>= 4.10.0)
 - `networkx`_
 - `tqdm`_
 

environment.yml

@@ -6,7 +6,7 @@ dependencies:
   - geopandas
   - inequality
   - jupyter
-  - libpysal>=4.8.0
+  - libpysal>=4.10.0
   - mapclassify
   - matplotlib
   - momepy

momepy/functional/_distribution.py

@@ -1,12 +1,21 @@
 import geopandas as gpd
 import numpy as np
+import pandas as pd
 import shapely
 from geopandas import GeoDataFrame, GeoSeries
 from libpysal.graph import Graph
 from packaging.version import Version
 from pandas import Series
+from scipy import sparse
 
-__all__ = ["orientation", "shared_walls", "alignment", "neighbor_distance"]
+__all__ = [
+    "orientation",
+    "shared_walls",
+    "alignment",
+    "neighbor_distance",
+    "mean_interbuilding_distance",
+    "building_adjacency",
+]
 
 GPD_GE_013 = Version(gpd.__version__) >= Version("0.13.0")
 
@@ -152,3 +161,116 @@ def neighbor_distance(geometry: GeoDataFrame | GeoSeries, graph: Graph) -> Serie
     )
     mean_distance.loc[graph.isolates] = np.nan
     return mean_distance
+
+
+def mean_interbuilding_distance(
+    geometry: GeoDataFrame | GeoSeries,
+    adjacency_graph: Graph,
+    neighborhood_graph: Graph,
+) -> Series:
+    """Calculate the mean distance between adjacent geometries within a set neighborhood
+
+    For each building, this function takes a neighborhood based on the neighbors within
+    a ``neighborhood_graph`` and calculates the mean distance between adjacent buildings
+    within this neighborhood where adjacency is captured by ``adjacency_graph``.
+
+    Notes
+    -----
+    The index of ``geometry`` must match the index along which both of the graphs are
+    built.
+
+    Parameters
+    ----------
+    geometry : GeoDataFrame | GeoSeries
+        A GeoDataFrame or GeoSeries containing geometries to analyse.
+    adjacency_graph : libpysal.graph.Graph
+        Graph representing the adjacency of geometries. Typically, this is a contiguity
+        graph derived from tessellation cells linked to buildings.
+    neighborhood_graph : libpysal.graph.Graph
+        Graph representing the extent around each geometry within which to calculate
+        the mean interbuilding distance. This can be a distance based graph, KNN graph,
+        higher order contiguity, etc.
+
+    Returns
+    -------
+    Series
+    """
+    distance = pd.Series(
+        shapely.distance(
+            geometry.geometry.loc[
+                adjacency_graph._adjacency.index.get_level_values(0)
+            ].values,
+            geometry.geometry.loc[
+                adjacency_graph._adjacency.index.get_level_values(1)
+            ].values,
+        ),
+        index=adjacency_graph._adjacency.index,
+        name="distance",
+    )
+
+    distance_matrix = (
+        distance.astype("Sparse[float]").sparse.to_coo(sort_labels=True)[0].tocsr()
+    )
+    neighborhood_matrix = sparse.coo_matrix(neighborhood_graph.sparse).tocsr()
+
+    mean_distances = np.zeros(distance_matrix.shape[0], dtype=float)
+
+    for i in range(distance_matrix.shape[0]):
+        neighborhood_indices = np.append(neighborhood_matrix[i].indices, i)
+        sub_matrix = distance_matrix[neighborhood_indices][:, neighborhood_indices]
+        mean_distances[i] = sub_matrix.sum() / sub_matrix.nnz
+
+    return Series(
+        mean_distances, index=geometry.index, name="mean_interbuilding_distance"
+    )
+
+
+def building_adjacency(
+    contiguity_graph: Graph,
+    neighborhood_graph: Graph,
+) -> Series:
+    """Calculate the level of building adjacency.
+
+    Building adjacency reflects how much buildings tend to join together into larger
+    structures. It is calculated as a ratio of joined built-up structures captured by
+    ``contiguity_graph`` and buildings within the neighborhood defined in
+    ``neighborhood_graph``.
+
+    Adapted from :cite:`vanderhaegen2017`.
+
+    Notes
+    -----
+    Both graphs must be built on the same index.
+
+    Parameters
+    ----------
+    contiguity_graph : libpysal.graph.Graph
+        Graph representing contiguity between geometries, typically a rook contiguity
+        graph derived from buildings.
+    neighborhood_graph : libpysal.graph.Graph
+        Graph representing the extent around each geometry within which to calculate
+        the level of building adjacency. This can be a distance based graph, KNN graph,
+        higher order contiguity, etc.
+
+    Returns
+    -------
+    Series
+    """
+    components = contiguity_graph.component_labels
+
+    # check if self-weights are present, otherwise assign them to treat self as part of
+    # the neighborhood
+    has_self_weights = (
+        neighborhood_graph._adjacency.index.get_level_values("focal")
+        == neighborhood_graph._adjacency.index.get_level_values("neighbor")
+    ).sum() == neighborhood_graph.n
+    if not has_self_weights:
+        neighborhood_graph = neighborhood_graph.assign_self_weight()
+
+    grouper = components.loc[
+        neighborhood_graph._adjacency.index.get_level_values(1)
+    ].groupby(neighborhood_graph._adjacency.index.get_level_values(0))
+    result = grouper.agg("nunique") / grouper.agg("count")
+    result.name = "building_adjacency"
+    result.index.name = None
+    return result

momepy/functional/tests/test_distribution.py

@@ -13,6 +13,8 @@ def setup_method(self):
         self.df_buildings = gpd.read_file(test_file_path, layer="buildings")
         self.df_streets = gpd.read_file(test_file_path, layer="streets")
         self.graph = Graph.build_knn(self.df_buildings.centroid, k=5)
+        self.contiguity = Graph.build_contiguity(self.df_buildings)
+        self.neighborhood_graph = self.graph.higher_order(3, lower_order=True)
 
     def test_orientation(self):
         expected = {
@@ -64,15 +66,45 @@ def test_neighbor_distance(self):
         r = mm.neighbor_distance(self.df_buildings, self.graph)
         assert_result(r, expected, self.df_buildings)
 
+    def test_mean_interbuilding_distance(self):
+        expected = {
+            "mean": 13.018190603684694,
+            "sum": 1874.6194469305958,
+            "min": 6.623582625492466,
+            "max": 22.513464171665948,
+        }
+        r = mm.mean_interbuilding_distance(
+            self.df_buildings, self.graph, self.neighborhood_graph
+        )
+        assert_result(r, expected, self.df_buildings)
+
+    def test_building_adjacency(self):
+        expected = {
+            "mean": 0.3784722222222222,
+            "sum": 54.5,
+            "min": 0.16666666666666666,
+            "max": 0.8333333333333334,
+        }
+        r = mm.building_adjacency(self.contiguity, self.graph)
+        assert_result(r, expected, self.df_buildings, exact=False)
+
 
 class TestEquality:
     def setup_method(self):
         test_file_path = mm.datasets.get_path("bubenec")
         self.df_buildings = gpd.read_file(test_file_path, layer="buildings").set_index(
             "uID"
         )
+        self.df_tessellation = gpd.read_file(
+            test_file_path, layer="tessellation"
+        ).set_index("uID")
         self.graph = Graph.build_knn(self.df_buildings.centroid, k=5)
         self.df_buildings["orientation"] = mm.orientation(self.df_buildings)
+        self.contiguity = Graph.build_contiguity(self.df_buildings)
+        self.tessellation_contiguity = Graph.build_contiguity(self.df_tessellation)
+        self.neighborhood_graph = self.tessellation_contiguity.higher_order(
+            3, lower_order=True
+        )
 
     def test_alignment(self):
         new = mm.alignment(self.df_buildings["orientation"], self.graph)
@@ -91,3 +123,22 @@ def test_neighbor_distance(self):
             self.df_buildings.reset_index(), self.graph.to_W(), "uID", verbose=False
         ).series
         assert_series_equal(new, old, check_names=False, check_index=False)
+
+    def test_mean_interbuilding_distance(self):
+        new = mm.mean_interbuilding_distance(
+            self.df_buildings, self.tessellation_contiguity, self.neighborhood_graph
+        )
+        old = mm.MeanInterbuildingDistance(
+            self.df_buildings.reset_index(),
+            self.tessellation_contiguity.to_W(),
+            "uID",
+            verbose=False,
+        ).series
+        assert_series_equal(new, old, check_names=False, check_index=False)
+
+    def test_building_adjacency(self):
+        new = mm.building_adjacency(self.contiguity, self.graph)
+        old = mm.BuildingAdjacency(
+            self.df_buildings.reset_index(), self.graph.to_W(), "uID", verbose=False
+        ).series
+        assert_series_equal(new, old, check_names=False, check_index=False)

momepy/functional/tests/test_shape.py

@@ -10,12 +10,12 @@
 GPD_013 = Version(gpd.__version__) >= Version("0.13")
 
 
-def assert_result(result, expected, geometry):
+def assert_result(result, expected, geometry, **kwargs):
     """Check the expected values and types of the result."""
     for key, value in expected.items():
         assert getattr(result, key)() == pytest.approx(value)
     assert isinstance(result, pd.Series)
-    assert_index_equal(result.index, geometry.index)
+    assert_index_equal(result.index, geometry.index, **kwargs)
 
 
 class TestShape:

pyproject.toml

@@ -35,7 +35,7 @@ classifiers = [
 requires-python = ">=3.10"
 dependencies = [
     "geopandas>=0.12.0",
-    "libpysal>=4.8.0",
+    "libpysal>=4.10.0",
     "networkx>=2.7",
     "packaging",
     "pandas!=1.5.0",