ENH: street_alignment and get_nearest_street (#566)

* ENH: street_alignment and get_nearest_street

* Apply suggestions from code review

Co-authored-by: James Gaboardi <jgaboardi@gmail.com>

* fix

---------

Co-authored-by: James Gaboardi <jgaboardi@gmail.com>

momepy/functional/_distribution.py

@@ -16,6 +16,7 @@
     "mean_interbuilding_distance",
     "building_adjacency",
     "neighbors",
+    "street_alignment",
 ]
 
 GPD_GE_013 = Version(gpd.__version__) >= Version("0.13.0")
@@ -279,7 +280,7 @@ def building_adjacency(
 
 def neighbors(
     geometry: GeoDataFrame | GeoSeries, graph: Graph, weighted=False
-) -> pd.Series:
+) -> Series:
     """Calculate the number of neighbours captured by ``graph``.
 
     If ``weighted=True``, the number of neighbours will be divided by the perimeter of
@@ -304,7 +305,7 @@ def neighbors(
 
     Returns
     -------
-    pd.Series
+    Series
     """
     if weighted:
         r = graph.cardinalities / geometry.length
@@ -313,3 +314,29 @@ def neighbors(
 
     r.name = "neighbors"
     return r
+
+
+def street_alignment(
+    building_orientation: Series,
+    street_orientation: Series,
+    street_index: Series,
+) -> Series:
+    """Calulate the deviation of the building orientation from the street orientation.
+
+    Parameters
+    ----------
+    building_orientation : Series
+        Series with the orientation of buildings. Can be measured using
+        :func:`orientation`.
+    street_orientation : Series
+        Series with the orientation of streets. Can be measured using
+        :func:`orientation`.
+    street_index : Series
+        Series with the index of the street to which the building belongs. Can be
+        retrieved using :func:`momepy.get_nearest_street`.
+
+    Returns
+    -------
+    Series
+    """
+    return (building_orientation - street_orientation.loc[street_index].values).abs()

momepy/functional/_elements.py

@@ -15,6 +15,7 @@
     "morphological_tessellation",
     "enclosed_tessellation",
     "verify_tessellation",
+    "get_nearest_street",
 ]
 
 
@@ -295,3 +296,44 @@ def verify_tessellation(tesselation, geometry):
             stacklevel=2,
         )
     return collapsed, multipolygons
+
+
+def get_nearest_street(
+    buildings: gpd.GeoSeries | gpd.GeoDataFrame,
+    streets: gpd.GeoSeries | gpd.GeoDataFrame,
+    max_distance: float | None = None,
+) -> np.ndarray:
+    """Identify the nearest street for each building.
+
+    Parameters
+    ----------
+    buildings : gpd.GeoSeries | gpd.GeoDataFrame
+        GeoSeries or GeoDataFrame of buildings
+    streets : gpd.GeoSeries | gpd.GeoDataFrame
+        GeoSeries or GeoDataFrame of streets
+    max_distance : float | None, optional
+        Maximum distance within which to query for nearest street. Must be
+        greater than 0. By default None, indicating no distance limit. Note that it is
+        advised to set a limit to avoid long processing times.
+
+    Notes
+    -----
+    In case of multiple streets within the same distance, only one is returned.
+
+    Returns
+    -------
+    np.ndarray
+        array containing the index of the nearest street for each building
+    """
+    blg_idx, str_idx = streets.sindex.nearest(
+        buildings.geometry, return_all=False, max_distance=max_distance
+    )
+
+    if streets.index.dtype == "object":
+        ids = np.empty(len(buildings), dtype=object)
+    else:
+        ids = np.empty(len(buildings), dtype=np.float32)
+        ids[:] = np.nan
+
+    ids[blg_idx] = streets.index[str_idx]
+    return ids

momepy/functional/tests/test_distribution.py

@@ -109,6 +109,19 @@ def test_neighbors(self):
         r = mm.neighbors(self.df_tessellation, self.tess_contiguity, weighted=True)
         assert_result(r, expected, self.df_buildings, exact=False, check_names=False)
 
+    def test_street_alignment(self):
+        building_orientation = mm.orientation(self.df_buildings)
+        street_orientation = mm.orientation(self.df_streets)
+        street_index = mm.get_nearest_street(self.df_buildings, self.df_streets)
+        expected = {
+            "mean": 2.024707906317863,
+            "sum": 291.5579385097722,
+            "min": 0.0061379200252815,
+            "max": 20.357934749623894,
+        }
+        r = mm.street_alignment(building_orientation, street_orientation, street_index)
+        assert_result(r, expected, self.df_buildings)
+
 
 class TestEquality:
     def setup_method(self):
@@ -119,8 +132,10 @@ def setup_method(self):
         self.df_tessellation = gpd.read_file(
             test_file_path, layer="tessellation"
         ).set_index("uID")
+        self.df_streets = gpd.read_file(test_file_path, layer="streets")
         self.graph = Graph.build_knn(self.df_buildings.centroid, k=5)
         self.df_buildings["orientation"] = mm.orientation(self.df_buildings)
+        self.df_streets["orientation"] = mm.orientation(self.df_streets)
         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(
@@ -188,3 +203,20 @@ def test_neighbors(self):
             verbose=False,
         ).series
         assert_series_equal(new, old, check_names=False, check_index=False)
+
+    def test_street_alignment(self):
+        street_index = mm.get_nearest_street(self.df_buildings, self.df_streets)
+        self.df_buildings["nID"] = street_index
+        new = mm.street_alignment(
+            self.df_buildings["orientation"],
+            self.df_streets["orientation"],
+            street_index,
+        )
+        old = mm.StreetAlignment(
+            self.df_buildings.reset_index(),
+            self.df_streets.reset_index(),
+            "orientation",
+            left_network_id="nID",
+            right_network_id="index",
+        ).series
+        assert_series_equal(new, old, check_names=False, check_index=False)

momepy/functional/tests/test_elements.py

@@ -148,3 +148,29 @@ def test_verify_tessellation(self):
         assert_index_equal(
             multi, pd.Index([1, 46, 57, 62, 103, 105, 129, 130, 134, 136, 137])
         )
+
+    def test_get_nearest_street(self):
+        streets = self.df_streets.copy()
+        nearest = mm.get_nearest_street(self.df_buildings, streets)
+        assert len(nearest) == len(self.df_buildings)
+        expected = np.array(
+            [0, 1, 2, 5, 6, 8, 10, 11, 12, 14, 16, 19, 21, 24, 25, 26, 28, 32, 33, 34]
+        )
+        expected_counts = np.array(
+            [9, 1, 12, 5, 7, 15, 1, 3, 4, 1, 3, 9, 9, 6, 5, 5, 15, 6, 10, 18]
+        )
+        unique, counts = np.unique(nearest, return_counts=True)
+        np.testing.assert_array_equal(unique, expected)
+        np.testing.assert_array_equal(counts, expected_counts)
+
+        # induce missing
+        nearest = mm.get_nearest_street(self.df_buildings, streets, 10)
+        expected = np.array([2.0, 34.0, np.nan])
+        expected_counts = np.array([3, 4, 137])
+        unique, counts = np.unique(nearest, return_counts=True)
+        np.testing.assert_array_equal(unique, expected)
+        np.testing.assert_array_equal(counts, expected_counts)
+
+        streets.index = streets.index.astype(str)
+        nearest = mm.get_nearest_street(self.df_buildings, streets, 10)
+        assert (nearest == None).sum() == 137  # noqa: E711