From 452298dc26426cbf282963a99b624114eeff2983 Mon Sep 17 00:00:00 2001
From: Martin Fleischmann <martin@martinfleischmann.net>
Date: Thu, 29 Feb 2024 18:44:41 +0100
Subject: [PATCH] ENH: add orientation and shared_walls functional versions
 (#553)

---
 momepy/__init__.py                           |  1 +
 momepy/functional/_distribution.py           | 91 ++++++++++++++++++++
 momepy/functional/tests/test_distribution.py | 41 +++++++++
 3 files changed, 133 insertions(+)
 create mode 100644 momepy/functional/_distribution.py
 create mode 100644 momepy/functional/tests/test_distribution.py

diff --git a/momepy/__init__.py b/momepy/__init__.py
index f9a41de5..60f51162 100644
--- a/momepy/__init__.py
+++ b/momepy/__init__.py
@@ -8,6 +8,7 @@
 from .diversity import *
 from .elements import *
 from .functional._dimension import *
+from .functional._distribution import *
 from .functional._shape import *
 from .graph import *
 from .intensity import *
diff --git a/momepy/functional/_distribution.py b/momepy/functional/_distribution.py
new file mode 100644
index 00000000..55732de5
--- /dev/null
+++ b/momepy/functional/_distribution.py
@@ -0,0 +1,91 @@
+import geopandas as gpd
+import numpy as np
+import shapely
+from geopandas import GeoDataFrame, GeoSeries
+from packaging.version import Version
+from pandas import Series
+
+__all__ = ["orientation", "shared_walls"]
+
+GPD_GE_013 = Version(gpd.__version__) >= Version("0.13.0")
+
+
+def orientation(geometry: GeoDataFrame | GeoSeries) -> Series:
+    """Calculate the orientation of objects.
+
+    The 'orientation' is defined as the deviation of orientation of the bounding
+    rectangle from cardinal directions. As such it is within a range 0 - 45. The
+    orientation of LineStrings is represented by the orientation of the line connecting
+    the first and the last point of the segment.
+
+    Adapted from :cite:`schirmer2015`.
+
+    Parameters
+    ----------
+    geometry : GeoDataFrame | GeoSeries
+        A GeoDataFrame or GeoSeries containing polygons to analyse.
+
+    Returns
+    -------
+    Series
+    """
+    geom_types = geometry.geom_type
+    poly_mask = geom_types.str.contains("Polygon")
+    line_mask = geom_types.str.contains("Line")
+
+    result = np.full(len(geometry), np.nan, dtype=float)
+
+    if poly_mask.any():
+        bboxes = shapely.minimum_rotated_rectangle(geometry.geometry.loc[poly_mask])
+        coords = shapely.get_coordinates(bboxes)
+        pt0 = coords[::5]
+        pt1 = coords[1::5]
+        angle = np.arctan2(pt1[:, 0] - pt0[:, 0], pt1[:, 1] - pt0[:, 1])
+        result[poly_mask] = np.degrees(angle)
+
+    if line_mask.any():
+        first = shapely.get_point(geometry.geometry, 0)
+        last = shapely.get_point(geometry.geometry, -1)
+        pt0 = shapely.get_coordinates(first)
+        pt1 = shapely.get_coordinates(last)
+        angle = np.arctan2(pt1[:, 0] - pt0[:, 0], pt1[:, 1] - pt0[:, 1])
+        result[line_mask] = np.degrees(angle)
+
+    return Series(
+        np.abs((result + 45) % 90 - 45),
+        index=geometry.index,
+        dtype=float,
+        name="orientation",
+    )
+
+
+def shared_walls(geometry: GeoDataFrame | GeoSeries) -> Series:
+    """Calculate the length of shared walls of adjacent elements (typically buildings).
+
+    Note that data needs to be topologically correct. Overlapping polygons will lead to
+    incorrect results.
+
+    Adapted from :cite:`hamaina2012a`.
+
+    Parameters
+    ----------
+    geometry : GeoDataFrame | GeoSeries
+        A GeoDataFrame or GeoSeries containing polygons to analyse.
+
+    Returns
+    -------
+    Series
+    """
+    if GPD_GE_013:
+        inp, res = geometry.sindex.query(geometry.geometry, predicate="touches")
+    else:
+        inp, res = geometry.sindex.query_bulk(geometry.geometry, predicate="touches")
+    left = geometry.geometry.take(inp).reset_index(drop=True)
+    right = geometry.geometry.take(res).reset_index(drop=True)
+    intersections = left.intersection(right).length
+    walls = intersections.groupby(inp).sum()
+    walls.index = geometry.index.take(walls.index)
+
+    results = Series(0.0, index=geometry.index, name="shared_walls")
+    results.loc[walls.index] = walls
+    return results
diff --git a/momepy/functional/tests/test_distribution.py b/momepy/functional/tests/test_distribution.py
new file mode 100644
index 00000000..50657fa6
--- /dev/null
+++ b/momepy/functional/tests/test_distribution.py
@@ -0,0 +1,41 @@
+import geopandas as gpd
+
+import momepy as mm
+
+from .test_shape import assert_result
+
+
+class TestDistribution:
+    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_streets = gpd.read_file(test_file_path, layer="streets")
+
+    def test_orientation(self):
+        expected = {
+            "mean": 20.983859394267952,
+            "sum": 3021.6757527745854,
+            "min": 7.968673890244247,
+            "max": 42.329365250279125,
+        }
+        r = mm.orientation(self.df_buildings)
+        assert_result(r, expected, self.df_buildings)
+
+        expected = {
+            "mean": 21.176405050561755,
+            "sum": 741.1741767696615,
+            "min": 0.834911325974133,
+            "max": 44.83357900046826,
+        }
+        r = mm.orientation(self.df_streets)
+        assert_result(r, expected, self.df_streets)
+
+    def test_shared_walls(self):
+        expected = {
+            "mean": 36.87618331446485,
+            "sum": 5310.17039728293,
+            "min": 0,
+            "max": 106.20917523555639,
+        }
+        r = mm.shared_walls(self.df_buildings)
+        assert_result(r, expected, self.df_buildings)