add recreational space per capita indicator

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):

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

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

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