CIF-221 add recreational space per capita indicator

city_metrix/layers/layer.py

@@ -127,7 +127,7 @@ class LayerGroupBy:
     def __init__(self, aggregate, zones, layer=None, masks=[]):
         self.aggregate = aggregate
         self.masks = masks
-        self.zones = zones.reset_index()
+        self.zones = zones.reset_index(drop=True)
         self.layer = layer
 
     def mean(self):
@@ -136,6 +136,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])
@@ -315,6 +318,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/layers/open_street_map.py

@@ -7,6 +7,13 @@
 
 
 class OpenStreetMapClass(Enum):
+    # ALL includes all 29 primary features https://wiki.openstreetmap.org/wiki/Map_features
+    ALL = {'aerialway': True, 'aeroway': True, 'amenity': True, 'barrier': True, 'boundary': True, 
+           'building': True, 'craft': True, 'emergency': True, 'geological': True, 'healthcare': True,
+           'highway': True, 'historic': True, 'landuse': True, 'leisure': True, 'man_made': True,
+           'military': True, 'natural': True, 'office': True, 'place': True, 'power': True,
+           'public_transport': True, 'railway': True, 'route': True, 'shop': True, 'sport': True,
+           'telecom': True, 'tourism': True, 'water': True, 'waterway': True}
     OPEN_SPACE = {'leisure': ['park', 'nature_reserve', 'common', 'playground', 'pitch', 'track'],
                   'boundary': ['protected_area', 'national_park']}
     OPEN_SPACE_HEAT = {'leisure': ['park', 'nature_reserve', 'common', 'playground', 'pitch', 'garden', 'golf_course', 'dog_park', 'recreation_ground', 'disc_golf_course'],
@@ -23,7 +30,7 @@ class OpenStreetMapClass(Enum):
 							'building': ['office', 'commercial', 'industrial', 'retail', 'supermarket'],
 							'shop': True}
     SCHOOLS = {'building': ['school',],
-				'amenity': ['school', 'kindergarten']}
+			   'amenity': ['school', 'kindergarten']}
     HIGHER_EDUCATION = {'amenity': ['college', 'university'],
 						'building': ['college', 'university']}
     TRANSIT_STOP = {'amenity':['ferry_terminal'],
@@ -35,19 +42,19 @@ class OpenStreetMapClass(Enum):
 
 
 class OpenStreetMap(Layer):
-    def __init__(self, osm_class=None, **kwargs):
+    def __init__(self, osm_class=OpenStreetMapClass.ALL, **kwargs):
         super().__init__(**kwargs)
         self.osm_class = osm_class
 
     def get_data(self, bbox):
-        north, south, east, west = bbox[3], bbox[1], bbox[0], bbox[2]
+        left, top, right, bottom = bbox
         # Set the OSMnx configuration to disable caching
         ox.settings.use_cache = False
         try:
-            osm_feature = ox.features_from_bbox(bbox=(north, south, east, west), tags=self.osm_class.value)
+            osm_feature = ox.features_from_bbox(bbox=(left, bottom, right, top), tags=self.osm_class.value)
         # When no feature in bbox, return an empty gdf
         except ox._errors.InsufficientResponseError as e:
-            osm_feature = gpd.GeoDataFrame(pd.DataFrame(columns=['osmid', 'geometry']+list(self.osm_class.value.keys())), geometry='geometry')
+            osm_feature = gpd.GeoDataFrame(pd.DataFrame(columns=['id', 'geometry']+list(self.osm_class.value.keys())), geometry='geometry')
             osm_feature.crs = "EPSG:4326"
 
         # Filter by geo_type
@@ -62,7 +69,7 @@ def get_data(self, bbox):
             osm_feature = osm_feature[osm_feature.geom_type.isin(['Polygon', 'MultiPolygon'])]
 
         # keep only columns desired to reduce file size
-        keep_col = ['osmid', 'geometry']
+        keep_col = ['id', 'geometry']
         for key in self.osm_class.value:
             if key in osm_feature.columns:
                 keep_col.append(key)

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

environment.yml

@@ -5,7 +5,8 @@ dependencies:
   - python=3.10
   - earthengine-api=0.1.411
   - geocube=0.4.2
-  - geopandas=0.14.4
+  - geopandas=1.0.1
+  - gdal=3.10.0
   - xarray=2024.7.0
   - rioxarray=0.15.0
   - odc-stac=0.3.8
@@ -14,7 +15,7 @@ dependencies:
   - xarray-spatial=0.3.7
   - xee=0.0.15
   - utm=0.7.0
-  - osmnx=1.9.0
+  - osmnx=2.0.0
   - dask[complete]=2023.11.0
   - matplotlib=3.8.2
   - jupyterlab=4.0.10

setup.py

@@ -17,8 +17,9 @@
         "xee",
         "rioxarray",
         "utm",
-        "osmnx",
+        "osmnx>=2.0.0",
         "geopandas",
+        "xarray",
         "s3fs",
         "dask>=2023.11.0",
         "boto3",

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