From 42bc1bbb16d4be23e7fb05cf3b082b650521e28c Mon Sep 17 00:00:00 2001
From: weiqi-tori <zhouwq@umich.edu>
Date: Tue, 26 Nov 2024 16:47:04 +0800
Subject: [PATCH] add recreational space per capita indicator

---
 city_metrix/layers/layer.py                      |  5 +++++
 city_metrix/metrics/__init__.py                  |  1 +
 .../metrics/recreational_space_per_capita.py     | 16 ++++++++++++++++
 tests/test_metrics.py                            |  7 +++++++
 4 files changed, 29 insertions(+)
 create mode 100644 city_metrix/metrics/recreational_space_per_capita.py

diff --git a/city_metrix/layers/layer.py b/city_metrix/layers/layer.py
index 83fa9eb..cad46bf 100644
--- a/city_metrix/layers/layer.py
+++ b/city_metrix/layers/layer.py
@@ -135,6 +135,9 @@ def mean(self):
     def count(self):
         return self._zonal_stats("count")
 
+    def sum(self):
+        return self._zonal_stats("sum")
+
     def _zonal_stats(self, stats_func):
         if box(*self.zones.total_bounds).area <= MAX_TILE_SIZE_DEGREES**2:
             stats = self._zonal_stats_tile(self.zones, [stats_func])
@@ -314,6 +317,8 @@ def _aggregate_stats(df, stats_func):
     elif stats_func == "mean":
         # mean must weight by number of pixels used for each tile
         return (df["mean"] * df["count"]).sum() / df["count"].sum()
+    elif stats_func == "sum":
+        return df["sum"].sum()
 
 
 def get_stats_funcs(stats_func):
diff --git a/city_metrix/metrics/__init__.py b/city_metrix/metrics/__init__.py
index d95cfa5..6efe2ae 100644
--- a/city_metrix/metrics/__init__.py
+++ b/city_metrix/metrics/__init__.py
@@ -5,3 +5,4 @@
 from .urban_open_space import urban_open_space
 from .natural_areas import natural_areas
 from .era_5_met_preprocessing import era_5_met_preprocessing
+from .recreational_space_per_capita import recreational_space_per_capita
diff --git a/city_metrix/metrics/recreational_space_per_capita.py b/city_metrix/metrics/recreational_space_per_capita.py
new file mode 100644
index 0000000..2e066d7
--- /dev/null
+++ b/city_metrix/metrics/recreational_space_per_capita.py
@@ -0,0 +1,16 @@
+from geopandas import GeoDataFrame, GeoSeries
+
+from city_metrix.layers import WorldPop, OpenStreetMap, OpenStreetMapClass
+
+
+def recreational_space_per_capita(zones: GeoDataFrame, spatial_resolution=100) -> GeoSeries:
+    world_pop = WorldPop(spatial_resolution=spatial_resolution)
+    open_space = OpenStreetMap(osm_class=OpenStreetMapClass.OPEN_SPACE)
+
+    # per 1000 people
+    world_pop_sum = world_pop.groupby(zones).sum() / 1000
+    # convert square meter to hectare
+    open_space_counts = open_space.mask(world_pop).groupby(zones).count()
+    open_space_area = open_space_counts.fillna(0) * spatial_resolution ** 2 / 10000
+
+    return open_space_area / world_pop_sum
diff --git a/tests/test_metrics.py b/tests/test_metrics.py
index 8fd42cc..a68b8ea 100644
--- a/tests/test_metrics.py
+++ b/tests/test_metrics.py
@@ -44,6 +44,13 @@ def test_natural_areas():
     assert expected_zone_size == actual_indicator_size
 
 
+def test_recreational_space_per_capita():
+    indicator = recreational_space_per_capita(ZONES)
+    expected_zone_size = ZONES.geometry.size
+    actual_indicator_size = indicator.size
+    assert expected_zone_size == actual_indicator_size
+
+
 def test_urban_open_space():
     indicator = urban_open_space(ZONES)
     expected_zone_size = ZONES.geometry.size