functional node_density implementation (#576)

* functional node_density implementation

* revert back to geodataframe functional implementation

* docstring change

* code review changes

* Update momepy/functional/_intensity.py

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

---------

Co-authored-by: Martin Fleischmann <martin@martinfleischmann.net>
Co-authored-by: James Gaboardi <jgaboardi@gmail.com>

momepy/functional/_intensity.py

@@ -1,9 +1,11 @@
+import numpy as np
+import pandas as pd
 import shapely
 from geopandas import GeoDataFrame, GeoSeries
 from libpysal.graph import Graph
 from pandas import Series
 
-__all__ = ["courtyards"]
+__all__ = ["courtyards", "node_density"]
 
 
 def courtyards(geometry: GeoDataFrame | GeoSeries, graph: Graph) -> Series:
@@ -48,3 +50,64 @@ def _calculate_courtyards(group):
     )
 
     return result
+
+
+def node_density(
+    nodes: GeoDataFrame, edges: GeoDataFrame, graph: Graph, weighted: bool = False
+) -> Series:
+    """Calculate the density of a node's neighbours (for all nodes)
+    on the street network defined in ``graph``.
+
+    Calculated as the number of neighbouring
+    nodes / cumulative length of street network within neighbours.
+    ``node_start``,  ``node_end``, is standard output of
+    :py:func:`momepy.nx_to_gdf` and is compulsory for ``edges`` to have
+    these columns.
+    If ``weighted``, a ``degree`` column is also required in ``nodes``.
+
+    Adapted from :cite:`dibble2017`.
+
+    Parameters
+    ----------
+    nodes : GeoDataFrame
+        A GeoDataFrame containing nodes of a street network.
+    edges : GeoDataFrame
+        A GeoDataFrame containing edges of a street network.
+    graph :  libpysal.graph.Graph
+        A spatial weights matrix capturing relationship between nodes.
+    weighted : bool (default False)
+        If ``True``, density will take into account node degree as ``k-1``.
+
+    Returns
+    -------
+    Series
+        A Series containing resulting values.
+
+    Examples
+    --------
+    >>> nodes['density'] = mm.node_density(nodes, edges, graph)
+    """
+
+    required_cols = ["node_start", "node_end"]
+    for col in required_cols:
+        if col not in edges.columns:
+            raise ValueError(f"Column {col} is needed in the edges GeoDataframe.")
+
+    if weighted and ("degree" not in nodes.columns):
+        raise ValueError("Column degree is needed in nodes GeoDataframe.")
+
+    def _calc_nodedensity(group, edges):
+        """Helper function to calculate group values."""
+        neighbours = group.index.values
+        locs = np.in1d(edges["node_start"], neighbours) & np.in1d(
+            edges["node_end"], neighbours
+        )
+        lengths = edges.loc[locs].geometry.length.sum()
+        return group.sum() / lengths if lengths else 0
+
+    if weighted:
+        summation_values = nodes["degree"] - 1
+    else:
+        summation_values = pd.Series(np.ones(nodes.shape[0]), index=nodes.index)
+
+    return graph.apply(summation_values, _calc_nodedensity, edges=edges)

momepy/functional/tests/test_intensity.py

@@ -1,5 +1,6 @@
 import geopandas as gpd
 import numpy as np
+import pytest
 from libpysal.graph import Graph
 from pandas.testing import assert_series_equal
 
@@ -38,6 +39,57 @@ def test_courtyards(self):
         expected = {"mean": 0.6805555555555556, "sum": 98, "min": 0, "max": 1}
         assert_result(courtyards, expected, self.df_buildings)
 
+    def test_node_density(self):
+        nx = mm.gdf_to_nx(self.df_streets, integer_labels=True)
+        nx = mm.node_degree(nx)
+        nodes, edges, w = mm.nx_to_gdf(nx, spatial_weights=True)
+        g = Graph.from_W(w).higher_order(k=3, lower_order=True).assign_self_weight()
+
+        density = mm.node_density(nodes, edges, g)
+        expected_density = {
+            "count": 29,
+            "mean": 0.005534125924228438,
+            "max": 0.010177844322387136,
+            "min": 0.00427032489140038,
+        }
+        assert_result(density, expected_density, nodes, check_names=False)
+
+        weighted = mm.node_density(nodes, edges, g, weighted=True)
+        expected_weighted = {
+            "count": 29,
+            "mean": 0.010090861332429164,
+            "max": 0.020355688644774272,
+            "min": 0.0077472994887720905,
+        }
+        assert_result(weighted, expected_weighted, nodes, check_names=False)
+
+        island = mm.node_density(nodes, edges, Graph.from_W(w).assign_self_weight())
+        expected_island = {
+            "count": 29,
+            "mean": 0.01026753724860306,
+            "max": 0.029319191032027746,
+            "min": 0.004808273240207287,
+        }
+        assert_result(island, expected_island, nodes, check_names=False)
+
+        with pytest.raises(
+            ValueError,
+            match=("Column node_start is needed in the edges GeoDataframe."),
+        ):
+            mm.node_density(nodes, nodes, g)
+
+        with pytest.raises(
+            ValueError,
+            match=("Column node_end is needed in the edges GeoDataframe."),
+        ):
+            mm.node_density(nodes, edges["node_start"].to_frame(), g)
+
+        with pytest.raises(
+            ValueError,
+            match=("Column degree is needed in nodes GeoDataframe."),
+        ):
+            mm.node_density(edges, edges, g, weighted=True)
+
 
 class TestIntensityEquality:
     def setup_method(self):
@@ -70,3 +122,32 @@ def test_courtyards(self):
         assert_series_equal(
             new_courtyards, old_courtyards, check_names=False, check_dtype=False
         )
+
+    def test_node_density(self):
+        nx = mm.gdf_to_nx(self.df_streets, integer_labels=True)
+        nx = mm.node_degree(nx)
+        nodes, edges, w = mm.nx_to_gdf(nx, spatial_weights=True)
+        sw = mm.sw_high(k=3, weights=w)
+        g = Graph.from_W(w).higher_order(k=3, lower_order=True).assign_self_weight()
+
+        density_old = mm.NodeDensity(nodes, edges, sw).series
+        density_new = mm.node_density(nodes, edges, g)
+        assert_series_equal(
+            density_old, density_new, check_names=False, check_dtype=False
+        )
+
+        weighted_old = mm.NodeDensity(
+            nodes, edges, sw, weighted=True, node_degree="degree"
+        ).series
+        weighted_new = mm.node_density(nodes, edges, g, weighted=True)
+        assert_series_equal(
+            weighted_old, weighted_new, check_names=False, check_dtype=False
+        )
+
+        islands_old = mm.NodeDensity(nodes, edges, w).series
+        islands_new = mm.node_density(
+            nodes, edges, Graph.from_W(w).assign_self_weight()
+        )
+        assert_series_equal(
+            islands_old, islands_new, check_names=False, check_dtype=False
+        )