# How to work with ERA5-Land on Earth Data Hub: storm Daniel flood, Greece 2023 [Earth Data Hub](https://earthdatahub.destine.eu) (EDH) offers an innovative and super-efficient way to access earth related data. This notebook will provide you guidance on how to access and use the [ERA5-Land](https://earthdatahub.destine.eu/collections/era5/datasets/reanalysis-era5-land) dataset on Earth Data Hub. ## Goal of this tutorial The first goal of this tutorial is to compute the total precipitation observed during the Storm Daniel event, from 6 to 7 September 2023, in Greece, and compare it with the average 1993-2022 precipitation in the same area. The second goal of this tutorial is to compare the 2023 cumulative precipitation on a specific interland location (39.25 N, 21.9 E) with the cumulative 1993-2022 precipitation in the same location. ## What you will learn: - how to access the dataset - select and reduce the data - plot the results --- To access datasets on Earth Data Hub you need to instruct your tools (xarray, Zarr, etc.) to use EDH personal access token when downloading the data. To obtain a personal access token you first need to register to the [Destination Earth platform](https://platform.destine.eu). Then you can go to [Earth Data Hub account settings](https://earthdatahub.destine.eu/account-settings) where you can find your default personal access token or create others. After retrieving your personal access token, please cut and paste it below: ⤵ ```python PAT = "your_EDH_personal_access_token" #e.g. PAT="edh_pat_44bbb7e9192a4c6bb47ddf07d07564eee5d17de8dfc48f7118f88e3bc4a4157f8fe2403f5aa0a2d53441b6922ea9a33a" ``` ## Working with EDH data Datasets on Earth Data Hub are often very large and remotely hosted. Typical use imply a selection of the data followed by one or more reduction steps to be performed in a local or distributed Dask environment. The structure of a workflow that uses EDH data typically looks like this: 1. Data access 2. Data selection 3. Data reduction (optional) 4. Data download 5. Further operations and visualisation [Xarray](https://docs.xarray.dev/en/stable) and [Dask](https://www.dask.org) work together following a lazy principle. This means that when you access and manipulate a Zarr store the data is in not immediately downloaded and loaded in memory. Instead, Dask constructs a task graph that represents the operations to be performed. A smart user will first reduce the amount of data that needs to be downloaded and explicitly call `compute()` on it. Once the `compute()` operation is complete the data is loaded into memory and available for subsequent fast processing. ### 1. Data access To access the data, only the dataset metadata must be downloaded. Xarray does this automatically when you access a Zarr dataset: ```python import xarray as xr ds = xr.open_dataset( f"https://edh:{PAT}@data.earthdatahub.destine.eu/era5/reanalysis-era5-land-no-antartica-v0.zarr", chunks={}, engine="zarr", ) ds ``` ```text Size: 709TB Dimensions: (valid_time: 669096, latitude: 1472, longitude: 3600) Coordinates: * valid_time (valid_time) datetime64[ns] 5MB 1950-01-01 ... 2026-... * latitude (latitude) float64 12kB 90.0 89.9 89.8 ... -57.0 -57.1 * longitude (longitude) float64 29kB 0.0 0.1 0.2 ... 359.8 359.9 depthBelowLandLayer float64 8B ... number int64 8B ... surface float64 8B ... Data variables: (12/50) asn (valid_time, latitude, longitude) float32 14TB dask.array d2m (valid_time, latitude, longitude) float32 14TB dask.array e (valid_time, latitude, longitude) float32 14TB dask.array es (valid_time, latitude, longitude) float32 14TB dask.array evabs (valid_time, latitude, longitude) float32 14TB dask.array evaow (valid_time, latitude, longitude) float32 14TB dask.array ... ... swvl4 (valid_time, latitude, longitude) float32 14TB dask.array t2m (valid_time, latitude, longitude) float32 14TB dask.array tp (valid_time, latitude, longitude) float32 14TB dask.array tsn (valid_time, latitude, longitude) float32 14TB dask.array u10 (valid_time, latitude, longitude) float32 14TB dask.array v10 (valid_time, latitude, longitude) float32 14TB dask.array Attributes: Conventions: CF-1.7 GRIB_centre: ecmf GRIB_centreDescription: European Centre for Medium-Range Weather Forecasts GRIB_edition: 1 GRIB_subCentre: 0 history: 2024-10-29T11:44 GRIB to CDM+CF via cfgrib-0.9.1... institution: European Centre for Medium-Range Weather Forecasts ``` ⚠ At this point, no data has been downloaded yet, nor loaded in memory. ### Storm Daniel precipitation (6-7 September 2023) vs average September precipitation (1993-2022) ### 2. Data selection First, we perform a geographical selection corresponding to the Greece area: ```python xr.set_options(keep_attrs=True) tp = ds.tp tp_greece = tp.sel(**{"latitude": slice(41, 34), "longitude": slice(19, 28)}) tp_greece ``` ```text Size: 17GB dask.array Coordinates: * valid_time (valid_time) datetime64[ns] 5MB 1950-01-01 ... 2026-... * latitude (latitude) float64 560B 40.9 40.8 40.7 ... 34.1 34.0 * longitude (longitude) float64 720B 19.0 19.1 19.2 ... 27.8 27.9 depthBelowLandLayer float64 8B ... number int64 8B ... surface float64 8B ... Attributes: (12/30) GRIB_NV: 0 GRIB_Nx: 3600 GRIB_Ny: 1472 GRIB_cfName: unknown GRIB_cfVarName: tp GRIB_dataType: fc ... ... GRIB_totalNumber: 0 GRIB_typeOfLevel: surface GRIB_units: m long_name: Total precipitation standard_name: unknown units: m ``` Then, we select only the data belonging to the month of September for the 1993-2023 time period (30 years). We will use this to extract the Storm Daniel event and the 1993-2022 September mean. ```python tp_greece_september = tp_greece[tp_greece.valid_time.dt.month.isin([9])] tp_greece_september_1993_2023 = tp_greece_september.sel(valid_time=slice('1993', '2023')) tp_greece_september_1993_2023 ``` ```text Size: 562MB dask.array Coordinates: * valid_time (valid_time) datetime64[ns] 179kB 1993-09-01 ... 202... * latitude (latitude) float64 560B 40.9 40.8 40.7 ... 34.1 34.0 * longitude (longitude) float64 720B 19.0 19.1 19.2 ... 27.8 27.9 depthBelowLandLayer float64 8B ... number int64 8B ... surface float64 8B ... Attributes: (12/30) GRIB_NV: 0 GRIB_Nx: 3600 GRIB_Ny: 1472 GRIB_cfName: unknown GRIB_cfVarName: tp GRIB_dataType: fc ... ... GRIB_totalNumber: 0 GRIB_typeOfLevel: surface GRIB_units: m long_name: Total precipitation standard_name: unknown units: m ``` At this point, the data should be small enough to fit in memory. We can check the size of it with the `costing.py` module: ```python import costing costing.estimate_download_size(tp_greece_september_1993_2023) ``` ```text Estimated download size: 0.475 GiB (486.0 MiB) ``` ### 3. Data download This is the phase where we explicitly trigger the download of the data. Remember to assign the return of the `compute()` function to a new variable, so that the data is kept in memory. We can measure the time it takes: ```python %%time tp_greece_september_1993_2023_computed = tp_greece_september_1993_2023.compute() ``` ```text CPU times: user 26.1 s, sys: 15.8 s, total: 41.9 s Wall time: 46.6 s ``` ### 4. Data reduction With the data loaded in memory we can easily compute the total precipitation for the Storm Daniel event. Notice that for each day of the year the total daily precipitation is given at T00:00:00. ```python tp_greece_storm_daniel = tp_greece_september_1993_2023_computed.sel(valid_time=["2023-09-06", "2023-09-07"]).sum("valid_time") # only takes the data at T00:00:00 timestamps tp_greece_storm_daniel ``` ```text Size: 25kB array([[0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.], ..., [0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.]], shape=(70, 90), dtype=float32) Coordinates: * latitude (latitude) float64 560B 40.9 40.8 40.7 ... 34.1 34.0 * longitude (longitude) float64 720B 19.0 19.1 19.2 ... 27.8 27.9 depthBelowLandLayer float64 8B 100.0 number int64 8B 0 surface float64 8B 0.0 Attributes: (12/30) GRIB_NV: 0 GRIB_Nx: 3600 GRIB_Ny: 1472 GRIB_cfName: unknown GRIB_cfVarName: tp GRIB_dataType: fc ... ... GRIB_totalNumber: 0 GRIB_typeOfLevel: surface GRIB_units: m long_name: Total precipitation standard_name: unknown units: m ``` ### 5. Visualisation Finally, we can plot the Storm Daniel event on a map: ```python import display import matplotlib.pyplot as plt import cartopy.feature as cfeature from cartopy import crs ``` ```python tp_greece_storm_daniel_mm = tp_greece_storm_daniel * 1000 # convert to [mm] tp_greece_storm_daniel_mm.attrs["units"] = "mm" plt.style.use("bmh") ax = display.map( tp_greece_storm_daniel_mm, projection=crs.Miller(), vmax=400, title="Storm Daniel total precipitation, 6-7 September 2023" ); ax.gridlines(draw_labels=True, alpha = 0.2) ``` ```text ``` ![image](https://earthdatahub.destine.eu/docs/_images/notebooks_insula-code-lab_tutorial-era5-land-floods-precipitation-anomaly-in-greece_22_1.png) We want to compare the total precipitation observed during Storm Daniel with the average total precipitation observed in September between 1993 and 2023 (30 years). To do this we compute the 1993_2022 september mean: ```python tp_greece_september_1993_2022_mean_mm = 1000*tp_greece_september_1993_2023_computed[tp_greece_september_1993_2023_computed['valid_time'].dt.hour == 0].sel(valid_time=slice("1993", "2022")).sum("valid_time")/30 tp_greece_september_1993_2022_mean_mm ``` ```text Size: 25kB array([[0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.], ..., [0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.]], shape=(70, 90), dtype=float32) Coordinates: * latitude (latitude) float64 560B 40.9 40.8 40.7 ... 34.1 34.0 * longitude (longitude) float64 720B 19.0 19.1 19.2 ... 27.8 27.9 depthBelowLandLayer float64 8B 100.0 number int64 8B 0 surface float64 8B 0.0 Attributes: (12/30) GRIB_NV: 0 GRIB_Nx: 3600 GRIB_Ny: 1472 GRIB_cfName: unknown GRIB_cfVarName: tp GRIB_dataType: fc ... ... GRIB_totalNumber: 0 GRIB_typeOfLevel: surface GRIB_units: m long_name: Total precipitation standard_name: unknown units: m ``` Finally, we can plot the Storm Daniel event and the average September 1993-2022 precipitation side by side: ```python axs = display.maps( [tp_greece_storm_daniel_mm, tp_greece_september_1993_2022_mean_mm], projection=crs.Miller(), vmax=400, axs_set=[ {"title": "Storm Daniel precipitation, 6-7 September 2023"}, {"title": "Average September precipitation, 1993-2022"}, ], ) for ax in axs: ax.gridlines(draw_labels=True, alpha = 0.2) ``` ![image](https://earthdatahub.destine.eu/docs/_images/notebooks_insula-code-lab_tutorial-era5-land-floods-precipitation-anomaly-in-greece_26_0.png) ### Accumulated precipitation: 2023 vs 1993-2022 In this section we will compare the accumulated 2023 precipitation on an specific inland location (39.25 N, 21.9 E) with the accumulated precipitation of each year between 1993 and 2022 in the same location. We will also highlight the 1993-2022 accumulated mean. We follow the same logic as before for selecting the data, estimating the download cost, computing the selection and plotting the result. ```python tp_greece_inland_1993_2023 = ds.tp.sel(**{"latitude": 39.25, "longitude": 21.90, "method": "nearest"}).sel(valid_time=slice("1993", "2023")) tp_greece_inland_1993_2023 ``` ```text Size: 1MB dask.array Coordinates: * valid_time (valid_time) datetime64[ns] 2MB 1993-01-01 ... 2023-... depthBelowLandLayer float64 8B ... latitude float64 8B 39.2 longitude float64 8B 21.9 number int64 8B ... surface float64 8B ... Attributes: (12/30) GRIB_NV: 0 GRIB_Nx: 3600 GRIB_Ny: 1472 GRIB_cfName: unknown GRIB_cfVarName: tp GRIB_dataType: fc ... ... GRIB_totalNumber: 0 GRIB_typeOfLevel: surface GRIB_units: m long_name: Total precipitation standard_name: unknown units: m ``` ```python import costing costing.estimate_download_size(tp_greece_inland_1993_2023) ``` ```text Estimated download size: 0.844 GiB (864.0 MiB) ``` ```python %%time tp_greece_inland_1993_2023_computed = tp_greece_inland_1993_2023.compute() ``` ```text CPU times: user 10.5 s, sys: 5.18 s, total: 15.6 s Wall time: 10.8 s ``` Now that the data is loaded in memory, we can easily extract the precipitation values for each year: ```python import datetime tp_greece_inland_2023 = tp_greece_inland_1993_2023_computed[tp_greece_inland_1993_2023_computed['valid_time'].dt.hour == 0].sel(valid_time="2023") tp_greece_inland_1993_2022 = tp_greece_inland_1993_2023_computed[tp_greece_inland_1993_2023_computed['valid_time'].dt.hour == 0].sel(valid_time=slice("1993", "2022")) ``` ```python %load_ext autoreload %autoreload 2 ``` Using the `display.compare()` method we can plot the accumulated precipitation for all the years between 1993 and 2023, highlighting the $10^{th}$, $50^{th}$, $90^{th}$ percentiles and the 2023 data. ```python plt.style.use('fivethirtyeight') ax = display.compare( tp_greece_inland_2023, tp_greece_inland_1993_2022, time="valid_time", ylim=[-100, 1750], ) ax.set_title('Accumulated precipitation [mm], Greece (39.25 N, 21.90 E)', loc='left', fontsize = 14) plt.legend(['10$^{th}$ percetile', '90$^{th}$ percetile', '50$^{th}$ percetile', '2023'], prop={'size': 10}); ``` ![image](https://earthdatahub.destine.eu/docs/_images/notebooks_insula-code-lab_tutorial-era5-land-floods-precipitation-anomaly-in-greece_36_0.png)