Functional density

momepy/functional/_intensity.py

@@ -5,7 +5,7 @@
 from libpysal.graph import Graph
 from pandas import Series
 
-__all__ = ["courtyards", "node_density"]
+__all__ = ["courtyards", "node_density", "density"]
 
 
 def courtyards(geometry: GeoDataFrame | GeoSeries, graph: Graph) -> Series:
@@ -111,3 +111,56 @@
         summation_values = pd.Series(np.ones(nodes.shape[0]), index=nodes.index)
 
     return graph.apply(summation_values, _calc_nodedensity, edges=edges)
+
+
+def density(
+    values: Series | np.ndarray, areas: Series | np.ndarray, graph: Graph
+) -> Series:
+    """Calculate the gross density.
+
+    .. math::
+        \\frac{\\sum \\text {values}}{\\sum \\text {areas}}
+
+    Adapted from :cite:`dibble2017`.
+
+    Parameters
+    ----------
+    values : pd.Series | np.ndarray
+        The character values for density calculations.
+        The index is used to arrange the final results.
+    areas : np.array | pd.Series
+        The area values for the density calculations,
+        an ``np.ndarray``, or ``pd.Series``.
+    graph : libpysal.graph.Graph
+        A spatial weights matrix for the geodataframe,
+        it is used to denote adjacent elements.
+
+    Returns
+    -------
+    DataFrame
+
+
+    Examples
+    --------
+    >>> tessellation_df['floor_area_dens'] = mm.density(tessellation_df['floor_area'],
+    ...                                                 tessellation_df['area'],
+    ...                                                 graph)
+    """
+
+    if isinstance(values, np.ndarray):
+        values = pd.DataFrame(values)
+    elif isinstance(values, pd.Series):
+        values = values.to_frame()
+
+    if isinstance(areas, np.ndarray):
+        areas = pd.Series(values)
+
+    stats = graph.apply(
+        pd.concat((values, areas.rename("area")), axis=1),
+        lambda x: (x.loc[:, x.columns != "area"].sum() / x["area"].sum()),
+    )
+    result = pd.DataFrame(
+        np.full(values.shape, np.nan), index=values.index, columns=values.columns
+    )
+    result[values.columns] = stats[values.columns]
+    return result

momepy/functional/tests/test_intensity.py

@@ -39,6 +39,23 @@ def test_courtyards(self):
         expected = {"mean": 0.6805555555555556, "sum": 98, "min": 0, "max": 1}
         assert_result(courtyards, expected, self.df_buildings)
 
+    def test_density(self):
+        graph = (
+            Graph.build_contiguity(self.df_tessellation, rook=False)
+            .higher_order(k=3, lower_order=True)
+            .assign_self_weight()
+        )
+        dens_new = mm.density(
+            self.df_buildings["fl_area"], self.df_tessellation.geometry.area, graph
+        )
+        dens_expected = {
+            "count": 144,
+            "mean": 1.6615871155383324,
+            "max": 2.450536855278486,
+            "min": 0.9746481727569978,
+        }
+        assert_result(dens_new["fl_area"], dens_expected, self.df_tessellation)
+
     def test_node_density(self):
         nx = mm.gdf_to_nx(self.df_streets, integer_labels=True)
         nx = mm.node_degree(nx)
@@ -123,6 +140,29 @@ def test_courtyards(self):
             new_courtyards, old_courtyards, check_names=False, check_dtype=False
         )
 
+    def test_density(self):
+        sw = mm.sw_high(k=3, gdf=self.df_tessellation, ids="uID")
+        graph = (
+            Graph.build_contiguity(self.df_tessellation, rook=False)
+            .higher_order(k=3, lower_order=True)
+            .assign_self_weight()
+        )
+        dens_new = mm.density(
+            self.df_buildings["fl_area"], self.df_tessellation.geometry.area, graph
+        )
+
+        dens_old = mm.Density(
+            self.df_tessellation,
+            self.df_buildings["fl_area"],
+            sw,
+            "uID",
+            self.df_tessellation.area,
+        ).series
+
+        assert_series_equal(
+            dens_new["fl_area"], dens_old, 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)