Merge branch 'main' of https://github.com/wri/cities-cif into add-ind…

…icator-percentpop-euclidean-proximity-to-openspace

.github/workflows/dev_ci_cd_conda.yml

@@ -33,4 +33,4 @@ jobs:
         GOOGLE_APPLICATION_USER: ${{ secrets.GOOGLE_APPLICATION_USER }}
         GOOGLE_APPLICATION_CREDENTIALS: ${{ secrets.GOOGLE_APPLICATION_CREDENTIALS }}
       run: |
-          pytest tests
+          pytest tests

README.md

@@ -43,6 +43,7 @@ To run the module,
 
   1. You need access to Google Earth Engine
   2. Install <https://cloud.google.com/sdk/docs/install>
+  3. If you want to use the ERA5 layer, you need to install the [Climate Data Store (CDS) Application Program Interface (API)](https://cds.climate.copernicus.eu/api-how-to)
 
 ### Interactive development
 

city_metrix/layers/__init__.py

@@ -19,4 +19,5 @@
 from .alos_dsm import AlosDSM
 from .overture_buildings import OvertureBuildings
 from .nasa_dem import NasaDEM
+from .era_5_hottest_day import Era5HottestDay
 from .impervious_surface import ImperviousSurface

city_metrix/layers/era_5_hottest_day.py

@@ -0,0 +1,115 @@
+import ee
+from timezonefinder import TimezoneFinder
+from pytz import timezone
+from datetime import datetime
+import pytz
+import cdsapi
+import os
+import xarray as xr
+
+from .layer import Layer
+
+class Era5HottestDay(Layer):
+    def __init__(self, start_date="2023-01-01", end_date="2024-01-01", **kwargs):
+        super().__init__(**kwargs)
+        self.start_date = start_date
+        self.end_date = end_date
+
+    def get_data(self, bbox):
+        dataset = ee.ImageCollection("ECMWF/ERA5_LAND/HOURLY")
+
+        # Function to find the city mean temperature of each hour
+        def hourly_mean_temperature(image):
+            hourly_mean = image.select('temperature_2m').reduceRegion(
+                reducer=ee.Reducer.mean(),
+                geometry=ee.Geometry.BBox(*bbox),
+                scale=11132,
+                bestEffort=True
+            ).values().get(0)
+
+            return image.set('hourly_mean_temperature', hourly_mean)
+
+        era5 = ee.ImageCollection(dataset
+                                  .filterBounds(ee.Geometry.BBox(*bbox))
+                                  .filterDate(self.start_date, self.end_date)
+                                  .select('temperature_2m')
+                                  )
+
+        era5_hourly_mean = era5.map(hourly_mean_temperature)
+
+        # Sort the collection based on the highest temperature and get the first image
+        highest_temperature_day = era5_hourly_mean.sort('hourly_mean_temperature', False).first()
+        highest_temperature_day = highest_temperature_day.get('system:index').getInfo()
+
+        # system:index in format 20230101T00
+        year = highest_temperature_day[0:4]
+        month = highest_temperature_day[4:6]
+        day = highest_temperature_day[6:8]
+        time = highest_temperature_day[-2:]
+
+        min_lon, min_lat, max_lon, max_lat = bbox
+        center_lon = (min_lon + max_lon) / 2
+        center_lat = (min_lat + max_lat) / 2
+
+        # Initialize TimezoneFinder
+        tf = TimezoneFinder()
+        # Find the timezone of the center point
+        tz_name = tf.timezone_at(lng=center_lon, lat=center_lat)
+        # Get the timezone object
+        local_tz = timezone(tz_name)
+        # Define the UTC time
+        utc_time = datetime.strptime(f'{year}-{month}-{day} {time}:00:00', "%Y-%m-%d %H:%M:%S")
+
+        # Convert UTC time to local time
+        local_time = utc_time.replace(tzinfo=pytz.utc).astimezone(local_tz)
+        local_date = local_time.date()
+
+        utc_times = []
+        for i in range(0, 24):
+            local_time_hourly = local_tz.localize(datetime(local_date.year, local_date.month, local_date.day, i, 0))
+            utc_time_hourly = local_time_hourly.astimezone(pytz.utc)
+            utc_times.append(utc_time_hourly)
+
+        utc_dates = list(set([dt.date() for dt in utc_times]))
+
+        dataarray_list = []
+        c = cdsapi.Client()
+        for i in range(len(utc_dates)):
+            c.retrieve(
+                'reanalysis-era5-single-levels',
+                {
+                    'product_type': 'reanalysis',
+                    'variable': [
+                        '10m_u_component_of_wind', '10m_v_component_of_wind', '2m_dewpoint_temperature',
+                        '2m_temperature', 'clear_sky_direct_solar_radiation_at_surface', 'mean_surface_direct_short_wave_radiation_flux_clear_sky',
+                        'mean_surface_downward_long_wave_radiation_flux_clear_sky', 'sea_surface_temperature', 'total_precipitation',
+                    ],
+                    'year': utc_dates[i].year,
+                    'month': utc_dates[i].month,
+                    'day': utc_dates[i].day,
+                    'time': ['00:00', '01:00', '02:00', '03:00', '04:00', '05:00', '06:00', '07:00', '08:00', '09:00',
+                             '10:00', '11:00', '12:00', '13:00', '14:00', '15:00', '16:00', '17:00', '18:00', '19:00',
+                             '20:00', '21:00', '22:00', '23:00'],
+                    'area': [max_lat, min_lon, min_lat, max_lon],
+                    'data_format': 'netcdf',
+                    'download_format': 'unarchived'
+                },
+                f'download_{i}.nc')
+
+            dataarray = xr.open_dataset(f'download_{i}.nc')
+
+            # Subset times for the day
+            times = [valid_time.astype('datetime64[s]').astype(datetime).replace(tzinfo=pytz.UTC) for valid_time in dataarray['valid_time'].values]
+            indices = [i for i, value in enumerate(times) if value in utc_times]
+            subset_dataarray = dataarray.isel(valid_time=indices)
+
+            dataarray_list.append(subset_dataarray)
+
+            # Remove local file
+            os.remove(f'download_{i}.nc')
+
+        data = xr.concat(dataarray_list, dim='valid_time')
+        # xarray.Dataset to xarray.DataArray
+        data = data.to_array()
+
+        return data

city_metrix/layers/layer.py

@@ -6,10 +6,11 @@
 
 import ee
 import boto3
+import math
 from dask.diagnostics import ProgressBar
 from ee import ImageCollection
 from geocube.api.core import make_geocube
-from shapely.geometry import box
+from shapely.geometry import box, polygon
 from xrspatial import zonal_stats
 import geopandas as gpd
 import xarray as xr
@@ -54,35 +55,74 @@ def groupby(self, zones, layer=None):
         """
         return LayerGroupBy(self.aggregate, zones, layer, self.masks)
 
-    def write(self, bbox, output_path, tile_degrees=None, **kwargs):
+    def write(self, bbox, output_path, tile_degrees=None, buffer_meters=None, **kwargs):
         """
         Write the layer to a path. Does not apply masks.
 
         :param bbox: (min x, min y, max x, max y)
         :param output_path: local or s3 path to output to
         :param tile_degrees: optional param to tile the results into multiple files with a VRT.
             Degrees to tile by. `output_path` should be a folder path to store the tiles.
+        : param buffer_meters: tile buffer distance in meters
         :return:
         """
 
-        if tile_degrees is not None:
-            tiles = create_fishnet_grid(*bbox, tile_degrees)
+        if tile_degrees is None:
+            clipped_data = self.aggregate.get_data(bbox)
+            write_layer(output_path, clipped_data)
+        else:
+            tile_grid, unbuffered_tile_grid = _get_tile_boundaries(bbox, tile_degrees, buffer_meters)
 
+            # write raster data to files
             if not os.path.exists(output_path):
                 os.makedirs(output_path)
 
-            file_names = []
-            for tile in tiles["geometry"]:
-                data = self.aggregate.get_data(tile.bounds)
+            for tile in tile_grid:
+                tile_name = tile['tile_name']
+                tile_geom = tile['geometry']
 
-                file_name = f"{output_path}/{uuid4()}.tif"
-                file_names.append(file_name)
+                file_path = os.path.join(output_path, tile_name)
+                clipped_data = self.aggregate.get_data(tile_geom.bounds)
+                write_layer(file_path, clipped_data)
 
-                write_layer(file_name, data)
-        else:
-            data = self.aggregate.get_data(bbox)
-            write_layer(output_path, data)
+            # write tile grid to geojson file
+            _write_tile_grid(tile_grid, output_path, 'tile_grid.geojson')
+
+            # if tiles were buffered, also write unbuffered tile grid to geojson file
+            if unbuffered_tile_grid:
+                _write_tile_grid(unbuffered_tile_grid, output_path, 'tile_grid_unbuffered.geojson')
+
+
+def _get_tile_boundaries(bbox, tile_degrees, buffer_meters):
+    has_buffer = True if buffer_meters is not None and buffer_meters != 0 else False
+    if has_buffer:
+        lon_degree_offset, lat_degree_offset = meters_to_offset_degrees(bbox, buffer_meters)
+        tiles = create_fishnet_grid(*bbox, tile_degrees, lon_degree_offset, lat_degree_offset)
+        unbuffered_tiles = create_fishnet_grid(*bbox, tile_degrees)
+    else:
+        tiles = create_fishnet_grid(*bbox, tile_degrees)
+        unbuffered_tiles = None
+
+    tile_grid = []
+    unbuffered_tile_grid = []
+    for index in range(0, len(tiles)):
+        tile_serial_id = index + 1
+        tile_suffix = str(tile_serial_id).zfill(3)
+        tile_name = f'tile_{tile_suffix}.tif'
 
+        tile_geom = tiles.iloc[index]['geometry']
+        tile_grid.append({"tile_name": tile_name, "geometry": tile_geom})
+
+        if has_buffer:
+            unbuffered_tile_geom = unbuffered_tiles.iloc[index]['geometry']
+            unbuffered_tile_grid.append({"tile_name": tile_name, "geometry": unbuffered_tile_geom})
+
+    return tile_grid, unbuffered_tile_grid
+
+def _write_tile_grid(tile_grid, output_path, target_file_name):
+    tile_grid = gpd.GeoDataFrame(tile_grid, crs='EPSG:4326')
+    tile_grid_file_path = str(os.path.join(output_path, target_file_name))
+    tile_grid.to_file(tile_grid_file_path)
 
 class LayerGroupBy:
     def __init__(self, aggregate, zones, layer=None, masks=[]):
@@ -222,20 +262,24 @@ def get_utm_zone_epsg(bbox) -> str:
     return f"EPSG:{epsg}"
 
 
-def create_fishnet_grid(min_x, min_y, max_x, max_y, cell_size):
+def create_fishnet_grid(min_x, min_y, max_x, max_y, cell_size, lon_degree_buffer=0, lat_degree_buffer=0):
     x, y = (min_x, min_y)
     geom_array = []
 
     # Polygon Size
     while y < max_y:
         while x < max_x:
+            cell_min_x = x - lon_degree_buffer
+            cell_min_y = y - lat_degree_buffer
+            cell_max_x = x + cell_size + lon_degree_buffer
+            cell_max_y = y + cell_size + lat_degree_buffer
             geom = geometry.Polygon(
                 [
-                    (x, y),
-                    (x, y + cell_size),
-                    (x + cell_size, y + cell_size),
-                    (x + cell_size, y),
-                    (x, y),
+                    (cell_min_x, cell_min_y),
+                    (cell_min_x, cell_max_y),
+                    (cell_max_x, cell_max_y),
+                    (cell_max_x, cell_min_y),
+                    (cell_min_x, cell_min_y),
                 ]
             )
             geom_array.append(geom)
@@ -327,3 +371,15 @@ def write_dataarray(path, data):
         os.remove(tmp_path)
     else:
         data.rio.to_raster(raster_path=path, driver="COG")
+
+def meters_to_offset_degrees(bbox, offset_meters):
+    min_lat = bbox[1]
+    max_lat = bbox[3]
+    center_lat = (min_lat + max_lat) / 2
+
+    meters_per_degree_of_lat = 111111
+    earth_circumference_meters = 40075000
+    lon_degree_offset = offset_meters/ (earth_circumference_meters * math.cos(center_lat) / 360)
+    lat_degree_offset = offset_meters / meters_per_degree_of_lat
+
+    return abs(lon_degree_offset), abs(lat_degree_offset)

city_metrix/layers/world_pop.py

@@ -5,25 +5,45 @@
 
 from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
+
 class WorldPop(Layer):
     """
     Attributes:
+        agesex_classes: list of age-sex classes to retrieve (see https://airtable.com/appDWCVIQlVnLLaW2/tblYpXsxxuaOk3PaZ/viwExxAgTQKZnRfWU/recFjH7WngjltFMGi?blocks=hide)
+        year: year used for data retrieval
         spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
     """
 
-    def __init__(self, spatial_resolution=100, **kwargs):
+    def __init__(self, agesex_classes=[], year=2020, spatial_resolution=100, **kwargs):
         super().__init__(**kwargs)
+        # agesex_classes options:
+        # M_0, M_1, M_5, M_10, M_15, M_20, M_25, M_30, M_35, M_40, M_45, M_50, M_55, M_60, M_65, M_70, M_75, M_80
+        # F_0, F_1, F_5, F_10, F_15, F_20, F_25, F_30, F_35, F_40, F_45, F_50, F_55, F_60, F_65, F_70, F_75, F_80
+        self.agesex_classes = agesex_classes
+        self.year = year
         self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
-        # load population
-        dataset = ee.ImageCollection('WorldPop/GP/100m/pop')
-        world_pop = ee.ImageCollection(dataset
-                     .filterBounds(ee.Geometry.BBox(*bbox))
-                     .filter(ee.Filter.inList('year', [2020]))
-                     .select('population')
-                     .mean()
-                     )
-
-        data = get_image_collection(world_pop, bbox, self.spatial_resolution, "world pop")
-        return data.population
+        if not self.agesex_classes:
+            # total population
+            dataset = ee.ImageCollection('WorldPop/GP/100m/pop')
+            world_pop = ee.ImageCollection(dataset
+                                           .filterBounds(ee.Geometry.BBox(*bbox))
+                                           .filter(ee.Filter.inList('year', [self.year]))
+                                           .select('population')
+                                           .mean()
+                                           )
+
+            data = get_image_collection(world_pop, bbox, self.spatial_resolution, "world pop").population
+
+        else:
+            # sum population for selected age-sex groups
+            dataset = ee.ImageCollection('WorldPop/GP/100m/pop_age_sex')
+            world_pop = dataset.filterBounds(ee.Geometry.BBox(*bbox))\
+                .filter(ee.Filter.inList('year', [self.year]))\
+                .select(self.agesex_classes)\
+                .mean()
+            world_pop = ee.ImageCollection(world_pop.reduce(ee.Reducer.sum()).rename('sum_age_sex_group'))
+            data = get_image_collection(world_pop, bbox, self.spatial_resolution, "world pop age sex").sum_age_sex_group
+
+        return data

city_metrix/metrics/__init__.py

@@ -1,6 +1,7 @@
 from .built_land_without_tree_cover import built_land_without_tree_cover
 from .built_land_with_low_surface_reflectivity import built_land_with_low_surface_reflectivity
 from .built_land_with_high_land_surface_temperature import built_land_with_high_land_surface_temperature
+from .era_5_met_preprocessing import era_5_met_preprocessing
 from .mean_tree_cover import mean_tree_cover
 from .natural_areas import natural_areas
 from .pop_open_space import pop_open_space