fix era5 issue

city_metrix/layers/era_5_hottest_day.py

@@ -6,6 +6,7 @@
 import cdsapi
 import os
 import xarray as xr
+import glob
 
 from .layer import Layer
 
@@ -71,8 +72,10 @@ def hourly_mean_temperature(image):
             utc_times.append(utc_time_hourly)
 
         utc_dates = list(set([dt.date() for dt in utc_times]))
-
-        dataarray_list = []
+
+        # {"dataType": "an"(analysis)/"fc"(forecast)/"pf"(perturbed forecast)}
+        an_list = []
+        fc_list = []
         c = cdsapi.Client()
         for i in range(len(utc_dates)):
             c.retrieve(
@@ -91,23 +94,44 @@ def hourly_mean_temperature(image):
                              '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',
+                    'data_format': 'grib',
                     'download_format': 'unarchived'
                 },
-                f'download_{i}.nc')
+                f'download_{i}.grib')
+
+            # {"dataType": "an"(analysis)/"fc"(forecast)/"pf"(perturbed forecast)}
+            with xr.open_dataset(f'download_{i}.grib',backend_kwargs={"filter_by_keys": {"dataType": "an"}}) as ds:
+                # Subset times for the day
+                times = [time.astype('datetime64[s]').astype(datetime).replace(tzinfo=pytz.UTC) for time in ds['time'].values]
+                indices = [i for i, value in enumerate(times) if value in utc_times]
+                subset_ds = ds.isel(time=indices).load()
 
-            with xr.open_dataset(f'download_{i}.nc') as dataarray:
+            an_list.append(subset_ds)
+
+            with xr.open_dataset(f'download_{i}.grib',backend_kwargs={"filter_by_keys": {"dataType": "fc"}}) as ds:
+                # reduce dimension
+                ds = ds.assign_coords(datetime=ds.time + ds.step)
+                ds = ds.stack(new_time=("time", "step"))
+                ds = ds.swap_dims({"new_time": "datetime"}).drop_vars(["time", "step", "new_time"])
+                ds = ds.rename(datetime="time")
                 # 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]
+                times = [time.astype('datetime64[s]').astype(datetime).replace(tzinfo=pytz.UTC) for time in ds['time'].values]
                 indices = [i for i, value in enumerate(times) if value in utc_times]
-                subset_dataarray = dataarray.isel(valid_time=indices).load()
+                subset_ds = ds.isel(time=indices).load()
+
+            fc_list.append(subset_ds)
+
+            # Remove local files
+            for file in glob.glob(f'download_{i}.grib*'):
+                os.remove(file)
 
-            dataarray_list.append(subset_dataarray)
+        an_data = xr.concat(an_list, dim='time')
+        fc_data = xr.concat(fc_list, dim='time')
+        fc_data = fc_data.sel(time=~fc_data.indexes['time'].duplicated())
+        fc_data = fc_data.transpose(*an_data.dims)
 
-            # Remove local file
-            os.remove(f'download_{i}.nc')
+        data =  xr.merge([an_data, fc_data], join="outer")
 
-        data = xr.concat(dataarray_list, dim='valid_time')
         # xarray.Dataset to xarray.DataArray
         data = data.to_array()
 

city_metrix/metrics/era_5_met_preprocessing.py

@@ -19,8 +19,8 @@ def era_5_met_preprocessing(zones: GeoDataFrame) -> GeoSeries:
     v10_var = era_5_data.sel(variable='v10').values
     sst_var = era_5_data.sel(variable='sst').values
     cdir_var = era_5_data.sel(variable='cdir').values
-    sw_var = era_5_data.sel(variable='msdrswrfcs').values
-    lw_var = era_5_data.sel(variable='msdwlwrfcs').values
+    sw_var = era_5_data.sel(variable='avg_sdirswrfcs').values
+    lw_var = era_5_data.sel(variable='avg_sdlwrfcs').values
     d2m_var = era_5_data.sel(variable='d2m').values
     time_var = era_5_data['valid_time'].values
     lat_var = era_5_data['latitude'].values

environment.yml

@@ -27,5 +27,6 @@ dependencies:
   - timezonefinder=6.5.2
   - scikit-image=0.24.0
   - exactextract=0.2.0
+  - cfgrib=0.9.15.0
   - pip:
       - overturemaps==0.6.0

setup.py

@@ -27,6 +27,7 @@
         "cdsapi",
         "timezonefinder",
         "scikit-image>=0.24.0",
-        "exactextract>=0.2.0"
+        "exactextract>=0.2.0",
+        "cfgrib"
     ],
 )