diff --git a/momepy/functional/_intensity.py b/momepy/functional/_intensity.py index 0b4fe6bb..09a6cf4c 100644 --- a/momepy/functional/_intensity.py +++ b/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) diff --git a/momepy/functional/tests/test_intensity.py b/momepy/functional/tests/test_intensity.py index 18e25122..2e305cd1 100644 --- a/momepy/functional/tests/test_intensity.py +++ b/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 + )