Download notebook Run this notebook

How to work with ERA5-Land on Earth Data Hub: storm Daniel flood, Greece 2023#

Earth Data Hub (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 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. Then you can go to Earth Data Hub 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: ⤵

[ ]:

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 and Dask 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:

[ ]:

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

<xarray.Dataset> 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<chunksize=(2880, 64, 64), meta=np.ndarray>
    d2m                  (valid_time, latitude, longitude) float32 14TB dask.array<chunksize=(2880, 64, 64), meta=np.ndarray>
    e                    (valid_time, latitude, longitude) float32 14TB dask.array<chunksize=(2880, 64, 64), meta=np.ndarray>
    es                   (valid_time, latitude, longitude) float32 14TB dask.array<chunksize=(2880, 64, 64), meta=np.ndarray>
    evabs                (valid_time, latitude, longitude) float32 14TB dask.array<chunksize=(2880, 64, 64), meta=np.ndarray>
    evaow                (valid_time, latitude, longitude) float32 14TB dask.array<chunksize=(2880, 64, 64), meta=np.ndarray>
    ...                   ...
    swvl4                (valid_time, latitude, longitude) float32 14TB dask.array<chunksize=(2880, 64, 64), meta=np.ndarray>
    t2m                  (valid_time, latitude, longitude) float32 14TB dask.array<chunksize=(2880, 64, 64), meta=np.ndarray>
    tp                   (valid_time, latitude, longitude) float32 14TB dask.array<chunksize=(2880, 64, 64), meta=np.ndarray>
    tsn                  (valid_time, latitude, longitude) float32 14TB dask.array<chunksize=(2880, 64, 64), meta=np.ndarray>
    u10                  (valid_time, latitude, longitude) float32 14TB dask.array<chunksize=(2880, 64, 64), meta=np.ndarray>
    v10                  (valid_time, latitude, longitude) float32 14TB dask.array<chunksize=(2880, 64, 64), meta=np.ndarray>
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:

[ ]:

xr.set_options(keep_attrs=True)

tp = ds.tp
tp_greece = tp.sel(**{"latitude": slice(41, 34), "longitude": slice(19, 28)})
tp_greece

<xarray.DataArray 'tp' (valid_time: 669096, latitude: 70, longitude: 90)> Size: 17GB
dask.array<getitem, shape=(669096, 70, 90), dtype=float32, chunksize=(2880, 49, 64), chunktype=numpy.ndarray>
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.

[ ]:

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

<xarray.DataArray 'tp' (valid_time: 22320, latitude: 70, longitude: 90)> Size: 562MB
dask.array<getitem, shape=(22320, 70, 90), dtype=float32, chunksize=(2871, 49, 64), chunktype=numpy.ndarray>
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:

[ ]:

import costing

costing.estimate_download_size(tp_greece_september_1993_2023)

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:

[ ]:

%%time

tp_greece_september_1993_2023_computed = tp_greece_september_1993_2023.compute()

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.

[ ]:

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

<xarray.DataArray 'tp' (latitude: 70, longitude: 90)> 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:

[ ]:

import display
import matplotlib.pyplot as plt
import cartopy.feature as cfeature

from cartopy import crs

[ ]:

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)

<cartopy.mpl.gridliner.Gridliner at 0x7ff46ed22c00>
../../_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:

[ ]:

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

<xarray.DataArray 'tp' (latitude: 70, longitude: 90)> 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:

[ ]:

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)
../../_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.

[ ]:

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

<xarray.DataArray 'tp' (valid_time: 271728)> Size: 1MB
dask.array<getitem, shape=(271728,), dtype=float32, chunksize=(2880,), chunktype=numpy.ndarray>
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
[ ]:

import costing

costing.estimate_download_size(tp_greece_inland_1993_2023)

Estimated download size: 0.844 GiB (864.0 MiB)

[ ]:

%%time

tp_greece_inland_1993_2023_computed = tp_greece_inland_1993_2023.compute()

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:

[ ]:

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"))

[ ]:

%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.

[ ]:

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});
../../_images/notebooks_insula-code-lab_tutorial-era5-land-floods-precipitation-anomaly-in-greece_36_0.png