Merge pull request #491 from digitalearthafrica/era5

ERA5 climate data
digitalearthafrica · Feb 19, 2024 · 81adb45 · 81adb45
2 parents 10dcc35 + c8dfab4
commit 81adb45
Show file tree

Hide file tree

Showing 9 changed files with 189 additions and 230 deletions.
diff --git a/Datasets/Climate_Data_ERA5_AWS.ipynb b/Datasets/Climate_Data_ERA5_AWS.ipynb
diff --git a/Datasets/Digital_Elevation_Models.ipynb b/Datasets/Digital_Elevation_Models.ipynb
diff --git a/Datasets/Landsat_Surface_Temperature.ipynb b/Datasets/Landsat_Surface_Temperature.ipynb
diff --git a/Datasets/Soil_Moisture.ipynb b/Datasets/Soil_Moisture.ipynb
diff --git a/Real_world_examples/Irrigated_cropping.ipynb b/Real_world_examples/Irrigated_cropping.ipynb
diff --git a/Real_world_examples/Temperature_trends.ipynb b/Real_world_examples/Temperature_trends.ipynb
diff --git a/Tools/deafrica_tools/load_era5.py b/Tools/deafrica_tools/load_era5.py
@@ -21,20 +21,14 @@
 
 
 ERA5_VARS = ['air_pressure_at_mean_sea_level',
- 'air_temperature_at_2_metres',
- 'air_temperature_at_2_metres_1hour_Maximum',
- 'air_temperature_at_2_metres_1hour_Minimum',
- 'dew_point_temperature_at_2_metres',
- 'eastward_wind_at_100_metres',
- 'eastward_wind_at_10_metres',
- 'integral_wrt_time_of_surface_direct_downwelling_shortwave_flux_in_air_1hour_Accumulation',
- 'lwe_thickness_of_surface_snow_amount',
- 'northward_wind_at_100_metres',
- 'northward_wind_at_10_metres',
- 'precipitation_amount_1hour_Accumulation',
- 'sea_surface_temperature',
- 'snow_density',
- 'surface_air_pressure']
+                 'air_temperature_at_2_metres',
+                 'eastward_wind_at_10_metres',
+                 'northward_wind_at_10_metres',
+                 'total_precipitation_6hr',
+                 'total_precipitation_12hr',
+                 'total_precipitation_24hr',
+                 'sea_surface_temperature',
+                 'surface_pressure']
 
 
 def load_era5(
@@ -101,39 +95,75 @@ def load_era5(
     # Loop through month and year to access ERA5 zarr
     month = date_from.astype('datetime64[M]')
     while month <= date_to.astype('datetime64[M]'):
-        url = f"s3://era5-pds/zarr/{month.astype(object).year:04}/{month.astype(object).month:02}/data/{var}.zarr"
-        ds = xr.open_zarr(fsspec.get_mapper(url, anon=True, 
+
+        try:
+
+            url = f"s3://era5-pds/zarr/{month.astype(object).year:04}/{month.astype(object).month:02}/data/{var}.zarr"
+            ds = xr.open_zarr(fsspec.get_mapper(url, anon=True, 
                                             client_kwargs={'region_name':'us-east-1'}),
-                          consolidated=True)
+                          consolidated=True)    
 
-        # re-order along longitude to go from -180 to 180 if needed
-        if min(lon) < 0:
-            ds = ds.assign_coords({"lon": (((ds.lon + 180) % 360) - 180)})
-            ds = ds.reindex({"lon": np.sort(ds.lon)})
-
-        if lat_range is None:
-            # find the nearest lat lon boundary points
-            test = ds.sel(lat=list(lat), lon=list(lon), method="nearest")
-            # define the lat/lon grid
-            lat_range = slice(test.lat.max().values, test.lat.min().values)
-            lon_range = slice(test.lon.min().values, test.lon.max().values)
-
-        if "time0" in ds.dims:
-            ds = ds.rename({"time0": "time"})
-        if "time1" in ds.dims:
-            ds = ds.rename(
-                {"time1": "time"}
-            )  # This should INTENTIONALLY error if both times are defined
-
-        output = ds[[var]].sel(lat=lat_range, lon=lon_range, time=slice(date_from, date_to)).resample(time=resample).reduce(reduce_func)
-        output.attrs = ds.attrs
-        for v in output.data_vars:
-            output[v].attrs = ds[v].attrs
-
-        datasets.append(output)
-        month += np.timedelta64(1,'M')
+            # re-order along longitude to go from -180 to 180 if needed
+            if min(lon) < 0:
+                ds = ds.assign_coords({"lon": (((ds.lon + 180) % 360) - 180)})
+                ds = ds.reindex({"lon": np.sort(ds.lon)})
+
+            if lat_range is None:
+                # find the nearest lat lon boundary points
+                test = ds.sel(lat=list(lat), lon=list(lon), method="nearest")
+                # define the lat/lon grid
+                lat_range = slice(test.lat.max().values, test.lat.min().values)
+                lon_range = slice(test.lon.min().values, test.lon.max().values)
+
+            if "time0" in ds.dims:
+                ds = ds.rename({"time0": "time"})
+            if "time1" in ds.dims:
+                ds = ds.rename(
+                    {"time1": "time"}
+                )  # This should INTENTIONALLY error if both times are defined
+
+            output = ds[[var]].sel(lat=lat_range, lon=lon_range, time=slice(date_from, date_to)).resample(time=resample).reduce(reduce_func)
+            output.attrs = ds.attrs
+            for v in output.data_vars:
+                output[v].attrs = ds[v].attrs
+
+            datasets.append(output)
+            month += np.timedelta64(1,'M')
+            return assign_crs(xr.combine_by_coords(datasets), 'EPSG:4326')
 
-    return assign_crs(xr.combine_by_coords(datasets), 'EPSG:4326')
+        except:
+            ERA5_dict = {'air_pressure_at_mean_sea_level':'mean_sea_level_pressure',
+                 'air_temperature_at_2_metres':'2m_temperature',
+                 'eastward_wind_at_10_metres':'10m_u_component_of_wind',
+                 'northward_wind_at_10_metres':'10m_v_component_of_wind',
+                 'total_precipitation_6hr': 'total_precipitation_6hr',
+                 'total_precipitation_12hr': 'total_precipitation_12hr',
+                 'total_precipitation_24hr':'total_precipitation_24hr',
+                 'sea_surface_temperature':'sea_surface_temperature',
+                 'surface_pressure':'surface_pressure'}
+
+            ds = xr.open_zarr(
+                    'gs://gcp-public-data-arco-era5/ar/1959-2022-wb13-6h-0p25deg-chunk-1.zarr-v2', 
+                    chunks={'time': 48},
+                    consolidated=True)[ERA5_dict[var]].sel(time=slice(date_from, date_to))
+
+            ds = ds.assign_coords(
+                    longitude=(((ds.longitude + 180) % 360) - 180)).sortby('longitude')
+
+            if lat_range is None:
+                # find the nearest lat lon boundary points
+                test = ds.sel(latitude=list(lat), longitude=list(lon), method="nearest")
+                # define the lat/lon grid
+                lat_range = slice(test.latitude.max().values, test.latitude.min().values)
+                lon_range = slice(test.longitude.min().values, test.longitude.max().values)
+
+            ds = ds.sel(latitude=lat_range, longitude = lon_range)
+
+            ds = ds.resample(time=resample).reduce(reduce_func)
+
+            return assign_crs(ds.to_dataset().rename({ERA5_dict[var]: var,
+                                                     'longitude':'lon',
+                                                     'latitude':'lat'}), 'EPSG:4326')
 
 
 # # older version of scripts to download and use netcdf

diff --git a/Tools/pyproject.toml b/Tools/pyproject.toml
@@ -6,7 +6,8 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "deafrica-tools"
-version = "2.3.4"
+version = "2.3.6"
+
 description = "Functions and algorithms for analysing Digital Earth Africa data."
 authors = [{name = "Digital Earth Africa", email = "[email protected]"}]
 maintainers = [{name = "Digital Earth Africa", email = "[email protected]"}]
@@ -73,6 +74,7 @@ dependencies= [
     "shapely",
     "tqdm",
     "xarray",
+    "gcsfs",
     # Code style
     "black",
     "flake8==6.0.0", # remove constraint in the next sandbox image upgrade

diff --git a/Use_cases/Okavango/3_Historical_rainfall.ipynb b/Use_cases/Okavango/3_Historical_rainfall.ipynb