add ERA 5 highest temperature layer

README.md

@@ -41,6 +41,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

@@ -18,3 +18,4 @@
 from .alos_dsm import AlosDSM
 from .overture_buildings import OvertureBuildings
 from .nasa_dem import NasaDEM
+from .era_5_hottest_day import Era5HottestDay

city_metrix/layers/era_5_hottest_day.py

@@ -0,0 +1,112 @@
+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],
+                    'format': 'netcdf',
+                },
+                f'download_{i}.nc')
+
+            dataarray = xr.open_dataset(f'download_{i}.nc')
+
+            # Subset times for the day
+            times = [time.astype('datetime64[s]').astype(datetime).replace(tzinfo=pytz.UTC) for time in dataarray['time'].values]
+            indices = [i for i, value in enumerate(times) if value in utc_times]
+            subset_dataarray = dataarray.isel(time=indices)
+
+            dataarray_list.append(subset_dataarray)
+
+            # Remove local file
+            os.remove(f'download_{i}.nc')
+
+        data = xr.concat(dataarray_list, dim='time')
+
+        return data

city_metrix/metrics/era_5_met_preprocessing.py

@@ -0,0 +1,69 @@
+
+from datetime import datetime
+import pandas as pd
+import numpy as np
+from geopandas import GeoDataFrame, GeoSeries
+
+from city_metrix.layers import Era5HottestDay
+
+
+def era_5_met_preprocessing(zones: GeoDataFrame) -> GeoSeries:
+    """
+    Get ERA 5 data for the hottest day
+    :param zones: GeoDataFrame with geometries to collect zonal stats on
+    :return: Pandas Dataframe of data
+    """
+    era_5_data = Era5HottestDay().get_data(zones.total_bounds)
+
+    t2m_var = era_5_data['t2m'].values
+    u10_var = era_5_data['u10'].values
+    v10_var = era_5_data['v10'].values
+    sst_var = era_5_data['sst'].values
+    cdir_var = era_5_data['cdir'].values
+    sw_var = era_5_data['msdrswrfcs'].values
+    lw_var = era_5_data['msdwlwrfcs'].values
+    d2m_var = era_5_data['d2m'].values
+    time_var = era_5_data['time'].values
+    lat_var = era_5_data['latitude'].values
+    lon_var = era_5_data['longitude'].values
+
+    # temps go from K to C; global rad (cdir) goes from /hour to /second; wind speed from vectors (pythagorean)
+    # rh calculated from temp and dew point; vpd calculated from tepm and rh
+    times = [time.astype('datetime64[s]').astype(datetime) for time in time_var]
+    t2m_vals = (t2m_var[:]-273.15)
+    d2m_vals = (d2m_var[:]-273.15)
+    rh_vals = (100*(np.exp((17.625*d2m_vals)/(243.04+d2m_vals))/np.exp((17.625*t2m_vals)/(243.04+t2m_vals))))
+    grad_vals = (cdir_var[:]/3600)
+    dir_vals = (sw_var[:])
+    dif_vals = (lw_var[:])
+    wtemp_vals = (sst_var[:]-273.15)
+    wind_vals = (np.sqrt(((np.square(u10_var[:]))+(np.square(v10_var[:])))))
+    # calc vapor pressure deficit in hPa for future utci conversion. first, get svp in pascals and then get vpd
+    svp_vals = (0.61078*np.exp(t2m_vals/(t2m_vals+237.3)*17.2694))
+    vpd_vals = ((svp_vals*(1-(rh_vals/100))))*10
+
+    # make lat/lon grid
+    latitudes = lat_var[:]
+    longitudes = lon_var[:]
+    latitudes_2d, longitudes_2d = np.meshgrid(latitudes, longitudes, indexing='ij')
+    latitudes_flat = latitudes_2d.flatten()
+    longitudes_flat = longitudes_2d.flatten()
+
+    # create pandas dataframe
+    df = pd.DataFrame({
+        'time': np.repeat(times, len(latitudes_flat)),
+        'lat': np.tile(latitudes_flat, len(times)),
+        'lon': np.tile(longitudes_flat, len(times)),
+        'temp': t2m_vals.flatten(),
+        'rh': rh_vals.flatten(),
+        'global_rad': grad_vals.flatten(),
+        'direct_rad': dir_vals.flatten(),
+        'diffuse_rad': dif_vals.flatten(),
+        'water_temp': wtemp_vals.flatten(),
+        'wind': wind_vals.flatten(),
+        'vpd': vpd_vals.flatten()
+    })
+    # round all numbers to two decimal places, which is the precision needed by the model
+    df = df.round(2)
+
+    return df

environment.yml

@@ -22,6 +22,9 @@ dependencies:
   - pip=23.3.1
   - boto3=1.34.124
   - scikit-learn=1.5.0
+  - cdsapi=0.7.0
+  - pytz=2024.1
+  - timezonefinder=6.5.2
   - pip:
       - git+https://github.com/isciences/exactextract
       - overturemaps==0.6.0

setup.py

@@ -31,5 +31,7 @@
         "exactextract",
         "overturemaps",
         "scikit-learn>=1.5.0",
+        "cdsapi",
+        "timezonefinder"
     ],
 )