Note

To run this notebook locally, follow the setup instructions in the GitHub repository.

Download notebook

Use case: ERA5 on a small area#

Earth Data Hub offers an innovative and super-efficient way to access data.

Here we present how to best access the service in the simplest use case.

Setup the environment#

If you haven’t done it already, follow the Getting started notebook to setup your environment and DestinE credentials.

Download 30 years of an ERA5 variable on a small area#

Our use case is to compute various climatological averages of one ERA5 variable over a long period of time, e.g. the surface temperature over 30 years. Note that we will be using the ERA5 single levels dataset that has more than 100 variable and every variable is almost 3TiB in size.

The best practice for downloading data form the Earth Data Hub comprise the following steps:

  1. open the dataset using the code snippet found on the ERA5 dataset page

  2. select the variable

  3. select the area of interest (optionally aligning it on chunk boundaries)

  4. select the time interval of interest

  5. download the data to memory (with .persist() or .compute())

  6. save the data or compute the result in memory

Open the dataset#

The following assumes you set up your EDH API key in your netrc file.

[1]:

import xarray as xr

era5_single_levels_dataset = xr.open_dataset(
    "https://data.earthdatahub.destine.eu/era5/era5-single-levels-atmosphere-v0.zarr",
    storage_options={"client_kwargs": {"trust_env": True}},
    chunks={},
    engine="zarr",
)
era5_single_levels_dataset

[1]:

<xarray.Dataset> Size: 407TB
Dimensions:     (valid_time: 760440, latitude: 721, longitude: 1440)
Coordinates:
  * valid_time  (valid_time) datetime64[ns] 6MB 1940-01-01 ... 2026-09-30T23:...
  * latitude    (latitude) float64 6kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude   (longitude) float64 12kB 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
Data variables: (12/129)
    alnid       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(1080, 60, 60), meta=np.ndarray>
    alnip       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(1080, 60, 60), meta=np.ndarray>
    aluvd       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(1080, 60, 60), meta=np.ndarray>
    aluvp       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(1080, 60, 60), meta=np.ndarray>
    anor        (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(1080, 60, 60), meta=np.ndarray>
    asn         (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(1080, 60, 60), meta=np.ndarray>
    ...          ...
    viiwn       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(1080, 60, 60), meta=np.ndarray>
    vilwd       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(1080, 60, 60), meta=np.ndarray>
    vilwe       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(1080, 60, 60), meta=np.ndarray>
    vilwn       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(1080, 60, 60), meta=np.ndarray>
    z           (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(1080, 60, 60), meta=np.ndarray>
    zust        (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(1080, 60, 60), meta=np.ndarray>
Attributes:
    GRIB_edition:            1
    GRIB_centre:             ecmf
    GRIB_centreDescription:  European Centre for Medium-Range Weather Forecasts
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             European Centre for Medium-Range Weather Forecasts
    history:                 2026-08-19T13:34 GRIB to CDM+CF via cfgrib-0.9.1...

Select and prepare the variable#

Note that all operation are lazy, that is are not applied to the whole 2.9TiB of data, but are just recorded for later use. Download and computations are only done when requested.

[2]:

t2m_world = era5_single_levels_dataset.t2m
t2m_world

[2]:

<xarray.DataArray 't2m' (valid_time: 760440, latitude: 721, longitude: 1440)> Size: 3TB
dask.array<open_dataset-t2m, shape=(760440, 721, 1440), dtype=float32, chunksize=(1080, 60, 60), chunktype=numpy.ndarray>
Coordinates:
  * valid_time  (valid_time) datetime64[ns] 6MB 1940-01-01 ... 2026-09-30T23:...
  * latitude    (latitude) float64 6kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude   (longitude) float64 12kB 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
Attributes: (12/32)
    GRIB_paramId:                             167
    GRIB_dataType:                            an
    GRIB_numberOfPoints:                      1038240
    GRIB_typeOfLevel:                         surface
    GRIB_stepUnits:                           1
    GRIB_stepType:                            instant
    ...                                       ...
    GRIB_totalNumber:                         0
    GRIB_units:                               K
    long_name:                                2 metre temperature
    units:                                    K
    standard_name:                            unknown
    last_ingested_valid_time_index:           760439

Select the area of interest#

For example we are interested in the area of the Alps.

This a very convenient example as:

  • it is a large enough area to be more interesting than a time-series of a single point,

  • it is small enough to fit a single Zarr chunk, so the notebook can be run even with a slow internet connection,

  • when plotted it is easy to identify even without adding coastlines and country borders that will make the notebook more complex.

[3]:

aoi_selection = {
    "latitude": slice(49, 45.5),  # NOTE: ERA5 has a descending latitude coordinate
    "longitude": slice(5, 14.5),
}
t2m_aoi = t2m_world.sel(aoi_selection)
t2m_aoi

[3]:

<xarray.DataArray 't2m' (valid_time: 760440, latitude: 15, longitude: 39)> Size: 2GB
dask.array<getitem, shape=(760440, 15, 39), dtype=float32, chunksize=(1080, 15, 39), chunktype=numpy.ndarray>
Coordinates:
  * valid_time  (valid_time) datetime64[ns] 6MB 1940-01-01 ... 2026-09-30T23:...
  * latitude    (latitude) float64 120B 49.0 48.75 48.5 ... 46.0 45.75 45.5
  * longitude   (longitude) float64 312B 5.0 5.25 5.5 5.75 ... 14.0 14.25 14.5
Attributes: (12/32)
    GRIB_paramId:                             167
    GRIB_dataType:                            an
    GRIB_numberOfPoints:                      1038240
    GRIB_typeOfLevel:                         surface
    GRIB_stepUnits:                           1
    GRIB_stepType:                            instant
    ...                                       ...
    GRIB_totalNumber:                         0
    GRIB_units:                               K
    long_name:                                2 metre temperature
    units:                                    K
    standard_name:                            unknown
    last_ingested_valid_time_index:           760439

Note that the data size is now much smaller, we went from 2.9TiB for global data to 1.7GiB of our small area of interest.

We plot the map of the temperature at one time to double check that the selection is correct:

[4]:

t2m_aoi.sel(valid_time="2020-01-01T00:00:00").plot(figsize=(6, 3))

[4]:

<matplotlib.collections.QuadMesh at 0x13622cad0>
../../_images/notebooks_edh-learning_10-use-case-era5-small-area_9_1.png

Optional: align the area of interest on chunk boundaries#

This a pro move and you can skip it.

When accessing the data in Zarr you always download whole chunks, even if you are only interested in part of them. In the case of ERA5 single level the spatial chunks are (60, 60) in size and even if you can read above that your DataArray in only 1.7GiB of data you are most probably downloading more data and then throwing a part of it away.

You can use the following (not very nice) code to grow your area of interest to the boundaries of the Zarr chunks, so you use all the data you download. Note that now you also have the size of the data that would be downloaded is it was not compressed (the data is in fact compressed so actual download is smaller).

[5]:

def align_indexer(step, indexer):
    if not isinstance(indexer, slice):
        return indexer
    assert indexer.step is None
    start = indexer.start // step * step if indexer.start else indexer.start
    stop = ((indexer.stop - 1) // step + 1) * step if indexer.stop else indexer.stop
    return slice(start, stop)


query_results = xr.core.indexing.map_index_queries(
    t2m_world, indexers=aoi_selection, method=None, tolerance=None
)
print(query_results.dim_indexers)

aoi_iselection = {
    dim: align_indexer(60, query_results.dim_indexers[dim])
    for dim in query_results.dim_indexers
}
print(aoi_iselection)

t2m_aoi = t2m_world.isel(aoi_iselection)
t2m_aoi

{'latitude': slice(164, 179, None), 'longitude': slice(20, 59, None)}
{'latitude': slice(120, 180, None), 'longitude': slice(0, 60, None)}

[5]:

<xarray.DataArray 't2m' (valid_time: 760440, latitude: 60, longitude: 60)> Size: 11GB
dask.array<getitem, shape=(760440, 60, 60), dtype=float32, chunksize=(1080, 60, 60), chunktype=numpy.ndarray>
Coordinates:
  * valid_time  (valid_time) datetime64[ns] 6MB 1940-01-01 ... 2026-09-30T23:...
  * latitude    (latitude) float64 480B 60.0 59.75 59.5 ... 45.75 45.5 45.25
  * longitude   (longitude) float64 480B 0.0 0.25 0.5 0.75 ... 14.25 14.5 14.75
Attributes: (12/32)
    GRIB_paramId:                             167
    GRIB_dataType:                            an
    GRIB_numberOfPoints:                      1038240
    GRIB_typeOfLevel:                         surface
    GRIB_stepUnits:                           1
    GRIB_stepType:                            instant
    ...                                       ...
    GRIB_totalNumber:                         0
    GRIB_units:                               K
    long_name:                                2 metre temperature
    units:                                    K
    standard_name:                            unknown
    last_ingested_valid_time_index:           760439

So, the full time-series for a whole chunk is really 10GiB, not 1.7GiB.

Let’s have a look at the area of the whole chunk.

[6]:

t2m_aoi.sel(valid_time="2020-01-01T00:00:00").plot()

[6]:

<matplotlib.collections.QuadMesh at 0x1362dbd90>
../../_images/notebooks_edh-learning_10-use-case-era5-small-area_13_1.png

Great! A single 60x60 chunk gets a big part of central Europe.

Select the time interval of interest#

We take one of the typical 30 year periods to compute climatologies and finally we have a definition of the data we want to download.

[7]:

t2m_aoi_toi = t2m_aoi.sel(valid_time=slice("1981", "2010"))
t2m_aoi_toi

[7]:

<xarray.DataArray 't2m' (valid_time: 262968, latitude: 60, longitude: 60)> Size: 4GB
dask.array<getitem, shape=(262968, 60, 60), dtype=float32, chunksize=(1080, 60, 60), chunktype=numpy.ndarray>
Coordinates:
  * valid_time  (valid_time) datetime64[ns] 2MB 1981-01-01 ... 2010-12-31T23:...
  * latitude    (latitude) float64 480B 60.0 59.75 59.5 ... 45.75 45.5 45.25
  * longitude   (longitude) float64 480B 0.0 0.25 0.5 0.75 ... 14.25 14.5 14.75
Attributes: (12/32)
    GRIB_paramId:                             167
    GRIB_dataType:                            an
    GRIB_numberOfPoints:                      1038240
    GRIB_typeOfLevel:                         surface
    GRIB_stepUnits:                           1
    GRIB_stepType:                            instant
    ...                                       ...
    GRIB_totalNumber:                         0
    GRIB_units:                               K
    long_name:                                2 metre temperature
    units:                                    K
    standard_name:                            unknown
    last_ingested_valid_time_index:           760439

From the representation above we learn a few things:

  1. the uncompressed data to be downloaded, e.g. 3.5GiB (compression depends on the dataset and the variable)

  2. the number of chunks to be downloaded, e.g. 245

Download the data to memory#

Finally we are ready to download only the data that we are interested in, in memory.

The best practice is to call the .compute() method load the data into a numpy array in memory.

This operation is the slow one. It can take up to 10 minutes.

It depends on the download speed of your internet connection and on the load on the Earth Data Hub. The closer you are to the data the faster it is, and this is one of the reason the EDH is best suited to be used from within the DestinE platform.

[8]:

%%time

t2m_aoi_toi_data = t2m_aoi_toi.compute()

CPU times: user 20.7 s, sys: 6.2 s, total: 26.9 s
Wall time: 3min 55s

Perform any computation#

Prepare the raw data#

[9]:

t2m_aoi_toi_data_C = t2m_aoi_toi_data - 273.15
t2m_aoi_toi_data_C.attrs["units"] = "°C"
t2m_aoi_toi_data_C.attrs["long_name"] = "surface temperature"
t2m_aoi_toi_data_C

[9]:

<xarray.DataArray 't2m' (valid_time: 262968, latitude: 60, longitude: 60)> Size: 4GB
array([[[ 4.100006 ,  4.100006 ,  3.850006 , ...,  0.3500061,
          0.3500061,  0.3500061],
        [ 3.850006 ,  3.850006 ,  3.850006 , ...,  0.3500061,
          0.3500061,  0.3500061],
        [ 3.850006 ,  4.100006 ,  4.350006 , ...,  0.6000061,
          0.6000061,  0.6000061],
        ...,
        [ 3.350006 ,  3.100006 ,  2.600006 , ..., -4.649994 ,
         -5.649994 , -6.399994 ],
        [ 3.350006 ,  2.850006 ,  2.600006 , ...,  0.6000061,
         -0.6499939, -1.6499939],
        [ 3.100006 ,  2.850006 ,  2.600006 , ...,  6.100006 ,
          4.850006 ,  1.6000061]],

       [[ 4.100006 ,  4.100006 ,  4.100006 , ...,  0.1000061,
          0.1000061,  0.3500061],
        [ 4.100006 ,  4.100006 ,  4.100006 , ...,  0.1000061,
          0.1000061,  0.3500061],
        [ 3.850006 ,  3.850006 ,  3.850006 , ...,  0.6000061,
          0.6000061,  0.8500061],
...
        [ 5.600006 ,  5.850006 ,  6.100006 , ..., -2.149994 ,
         -2.649994 , -4.149994 ],
        [ 5.600006 ,  5.600006 ,  5.850006 , ..., -1.8999939,
         -4.149994 , -5.649994 ],
        [ 6.100006 ,  6.350006 ,  6.350006 , ...,  2.100006 ,
         -0.3999939, -3.649994 ]],

       [[ 5.100006 ,  5.100006 ,  5.100006 , ..., -1.3999939,
         -1.1499939, -0.8999939],
        [ 5.100006 ,  5.100006 ,  5.100006 , ..., -0.3999939,
         -0.3999939, -0.1499939],
        [ 5.100006 ,  5.100006 ,  5.100006 , ...,  0.1000061,
          0.3500061,  0.6000061],
        ...,
        [ 5.100006 ,  5.350006 ,  5.600006 , ..., -2.399994 ,
         -2.399994 , -3.649994 ],
        [ 5.100006 ,  5.350006 ,  5.600006 , ..., -1.8999939,
         -3.899994 , -5.399994 ],
        [ 5.350006 ,  5.600006 ,  5.850006 , ...,  1.8500061,
         -0.6499939, -3.649994 ]]], shape=(262968, 60, 60), dtype=float32)
Coordinates:
  * valid_time  (valid_time) datetime64[ns] 2MB 1981-01-01 ... 2010-12-31T23:...
  * latitude    (latitude) float64 480B 60.0 59.75 59.5 ... 45.75 45.5 45.25
  * longitude   (longitude) float64 480B 0.0 0.25 0.5 0.75 ... 14.25 14.5 14.75
Attributes: (12/32)
    GRIB_paramId:                             167
    GRIB_dataType:                            an
    GRIB_numberOfPoints:                      1038240
    GRIB_typeOfLevel:                         surface
    GRIB_stepUnits:                           1
    GRIB_stepType:                            instant
    ...                                       ...
    GRIB_totalNumber:                         0
    GRIB_units:                               K
    long_name:                                surface temperature
    units:                                    °C
    standard_name:                            unknown
    last_ingested_valid_time_index:           760439

Monthly climatology#

[10]:

%%time

t2m_aoi_toi_month_mean = t2m_aoi_toi_data_C.groupby("valid_time.month").mean()
t2m_aoi_toi_month_mean

CPU times: user 1.42 s, sys: 1.56 s, total: 2.97 s
Wall time: 3.3 s

[10]:

<xarray.DataArray 't2m' (month: 12, latitude: 60, longitude: 60)> Size: 173kB
array([[[ 5.9818087 ,  5.9669843 ,  5.9524713 , ..., -4.6589036 ,
         -4.6404414 , -4.5902905 ],
        [ 5.966629  ,  5.9526863 ,  5.939012  , ..., -4.094842  ,
         -4.065144  , -3.9923122 ],
        [ 6.005804  ,  5.9744687 ,  5.9492054 , ..., -3.5102892 ,
         -3.5077865 , -3.4893947 ],
        ...,
        [ 5.372132  ,  5.2229385 ,  4.908522  , ..., -0.54323304,
         -0.6738021 , -0.7514413 ],
        [ 5.3958364 ,  5.315269  ,  4.998812  , ...,  0.5302754 ,
         -0.5924382 , -0.66884124],
        [ 5.5873632 ,  5.500552  ,  5.268542  , ...,  3.7213507 ,
          2.4584892 ,  0.1901969 ]],

       [[ 5.435061  ,  5.423069  ,  5.411557  , ..., -4.5158525 ,
         -4.515467  , -4.454127  ],
        [ 5.4172597 ,  5.404431  ,  5.390594  , ..., -3.9671335 ,
         -3.9808495 , -3.9163103 ],
        [ 5.449186  ,  5.420084  ,  5.3954844 , ..., -3.480913  ,
         -3.50844   , -3.4638371 ],
...
        [ 8.523972  ,  8.35806   ,  8.028974  , ...,  4.6578    ,
          4.5535216 ,  4.527012  ],
        [ 8.590157  ,  8.537646  ,  8.20685   , ...,  5.6508393 ,
          4.6596746 ,  4.6317544 ],
        [ 8.833831  ,  8.77931   ,  8.573055  , ...,  8.549518  ,
          7.447578  ,  5.4716263 ]],

       [[ 6.4842925 ,  6.4615417 ,  6.43884   , ..., -3.8836484 ,
         -3.8491683 , -3.776545  ],
        [ 6.469015  ,  6.4461217 ,  6.426031  , ..., -3.2930243 ,
         -3.2093337 , -3.113703  ],
        [ 6.519353  ,  6.482286  ,  6.4501595 , ..., -2.6695938 ,
         -2.6275818 , -2.6314178 ],
        ...,
        [ 5.94075   ,  5.8095894 ,  5.5204873 , ...,  0.6526263 ,
          0.48733467,  0.37250245],
        [ 5.992099  ,  5.9452276 ,  5.653638  , ...,  1.6936718 ,
          0.6220236 ,  0.508619  ],
        [ 6.1915736 ,  6.135068  ,  5.917882  , ...,  4.9268003 ,
          3.711836  ,  1.4372185 ]]], shape=(12, 60, 60), dtype=float32)
Coordinates:
  * month      (month) int64 96B 1 2 3 4 5 6 7 8 9 10 11 12
  * latitude   (latitude) float64 480B 60.0 59.75 59.5 ... 45.75 45.5 45.25
  * longitude  (longitude) float64 480B 0.0 0.25 0.5 0.75 ... 14.25 14.5 14.75
Attributes: (12/32)
    GRIB_paramId:                             167
    GRIB_dataType:                            an
    GRIB_numberOfPoints:                      1038240
    GRIB_typeOfLevel:                         surface
    GRIB_stepUnits:                           1
    GRIB_stepType:                            instant
    ...                                       ...
    GRIB_totalNumber:                         0
    GRIB_units:                               K
    long_name:                                surface temperature
    units:                                    °C
    standard_name:                            unknown
    last_ingested_valid_time_index:           760439
[11]:

t2m_aoi_toi_month_mean.plot(vmin=-30, col="month", col_wrap=3)

[11]:

<xarray.plot.facetgrid.FacetGrid at 0x136444d70>
../../_images/notebooks_edh-learning_10-use-case-era5-small-area_22_1.png

Daily climatological mean and quintile extremes#

[12]:

%%time

t2m_aoi_toi_doy_mean = t2m_aoi_toi_data_C.groupby("valid_time.dayofyear").mean()
t2m_aoi_toi_doy_mean = t2m_aoi_toi_doy_mean.sel(dayofyear=slice(0, 365))
t2m_aoi_toi_doy_mean

CPU times: user 1.28 s, sys: 519 ms, total: 1.8 s
Wall time: 1.82 s

[12]:

<xarray.DataArray 't2m' (dayofyear: 365, latitude: 60, longitude: 60)> Size: 5MB
array([[[ 5.904207  ,  5.870874  ,  5.839971  , ..., -5.089555  ,
         -5.082957  , -5.0829525 ],
        [ 5.8865    ,  5.862889  ,  5.8344164 , ..., -4.3419957 ,
         -4.3058844 , -4.2820997 ],
        [ 5.927126  ,  5.8861523 ,  5.8514304 , ..., -3.6145732 ,
         -3.6137052 , -3.6147463 ],
        ...,
        [ 5.1840544 ,  4.9882197 ,  4.6434255 , ..., -0.07325779,
         -0.28298035, -0.45034602],
        [ 5.2830133 ,  5.1284995 ,  4.779885  , ...,  1.1416657 ,
         -0.02152167, -0.18749389],
        [ 5.514266  ,  5.3503757 ,  5.0739846 , ...,  4.3979406 ,
          3.1308913 ,  0.9215325 ]],

       [[ 5.8250422 ,  5.791014  ,  5.7660136 , ..., -6.0864205 ,
         -6.0947537 , -6.0570803 ],
        [ 5.828862  ,  5.8003893 ,  5.784417  , ..., -5.4367695 ,
         -5.3806934 , -5.331736  ],
        [ 5.8896275 ,  5.8524737 ,  5.819487  , ..., -4.756582  ,
         -4.7326236 , -4.76995   ],
...
        [ 5.352464  ,  5.2437835 ,  5.005241  , ..., -0.36804944,
         -0.4121467 , -0.467355  ],
        [ 5.4253807 ,  5.404895  ,  5.1607966 , ...,  0.6833377 ,
         -0.46388277, -0.5246467 ],
        [ 5.683368  ,  5.6465626 ,  5.4590616 , ...,  3.9409862 ,
          2.644079  ,  0.34861723]],

       [[ 6.063586  ,  6.038239  ,  6.00942   , ..., -6.568361  ,
         -6.481903  , -6.380514  ],
        [ 6.0705314 ,  6.0451837 ,  6.032684  , ..., -5.656569  ,
         -5.4834795 , -5.3484106 ],
        [ 6.1399765 ,  6.1094213 ,  6.0799065 , ..., -4.6952934 ,
         -4.568905  , -4.5409536 ],
        ...,
        [ 5.254198  ,  5.155934  ,  4.9451675 , ..., -0.24929945,
         -0.3888828 , -0.5569388 ],
        [ 5.2965593 ,  5.3062816 ,  5.0916963 , ...,  0.91354203,
         -0.30034113, -0.46023697],
        [ 5.509407  ,  5.5073237 ,  5.38649   , ...,  4.1562657 ,
          2.8041487 ,  0.5128533 ]]], shape=(365, 60, 60), dtype=float32)
Coordinates:
  * dayofyear  (dayofyear) int64 3kB 1 2 3 4 5 6 7 ... 360 361 362 363 364 365
  * latitude   (latitude) float64 480B 60.0 59.75 59.5 ... 45.75 45.5 45.25
  * longitude  (longitude) float64 480B 0.0 0.25 0.5 0.75 ... 14.25 14.5 14.75
Attributes: (12/32)
    GRIB_paramId:                             167
    GRIB_dataType:                            an
    GRIB_numberOfPoints:                      1038240
    GRIB_typeOfLevel:                         surface
    GRIB_stepUnits:                           1
    GRIB_stepType:                            instant
    ...                                       ...
    GRIB_totalNumber:                         0
    GRIB_units:                               K
    long_name:                                surface temperature
    units:                                    °C
    standard_name:                            unknown
    last_ingested_valid_time_index:           760439
[13]:

point = {
    "latitude": 46.5,
    "longitude": 11.8,
    "method": "nearest",
}

t2m_point_doy_mean = t2m_aoi_toi_doy_mean.sel(**point)

t2m_point_doy_mean.plot()

[13]:

[<matplotlib.lines.Line2D at 0x133f6c6e0>]
../../_images/notebooks_edh-learning_10-use-case-era5-small-area_25_1.png 
[14]:

%%time

t2m_aoi_toi_doy_quantile = t2m_aoi_toi_data_C.groupby("valid_time.dayofyear").quantile(
    [0.9, 0.5, 0.1]
)
t2m_aoi_toi_doy_quantile = t2m_aoi_toi_doy_quantile.sel(dayofyear=slice(0, 365))
t2m_aoi_toi_doy_quantile

CPU times: user 35.6 s, sys: 406 ms, total: 36 s
Wall time: 36.1 s

[14]:

<xarray.DataArray 't2m' (dayofyear: 365, latitude: 60, longitude: 60,
                         quantile: 3)> Size: 32MB
array([[[[  8.1000061,   6.1000061,   3.3500061],
         [  8.1000061,   5.8500061,   3.3500061],
         [  8.1000061,   5.8500061,   3.3500061],
         ...,
         [  2.3500061,  -4.1499939, -12.9249939],
         [  2.1000061,  -4.1499939, -12.4249939],
         [  2.1000061,  -4.1499939, -12.4249939]],

        [[  8.1000061,   5.8500061,   3.6000061],
         [  8.1000061,   5.8500061,   3.3500061],
         [  8.1000061,   5.8500061,   3.3500061],
         ...,
         [  2.8500061,  -3.3999939, -11.1499939],
         [  2.6000061,  -3.3999939, -10.6499939],
         [  2.8500061,  -3.3999939, -10.3999939]],

        [[  8.1000061,   5.9750061,   3.6000061],
         [  8.1000061,   5.8500061,   3.6000061],
         [  8.1000061,   5.8500061,   3.3500061],
         ...,
...
         [  6.3500061,  -0.1499939,  -6.1499939],
         [  6.3500061,  -0.3999939,  -6.8999939],
         [  6.6000061,  -0.3999939,  -7.1499939]],

        [[ 11.3500061,   6.6000061,  -2.1499939],
         [ 11.1000061,   6.6000061,  -1.8999939],
         [ 10.6000061,   6.3500061,  -1.6749939],
         ...,
         [  7.1000061,   1.1000061,  -4.8999939],
         [  6.8500061,   0.1000061,  -7.6499939],
         [  6.6000061,  -0.1499939,  -7.6749939]],

        [[ 11.6000061,   6.6000061,  -1.6499939],
         [ 11.1250061,   6.8500061,  -1.4249939],
         [ 10.8500061,   6.8500061,  -1.3999939],
         ...,
         [  9.3500061,   4.2250061,   0.0750061],
         [  8.8750061,   3.1000061,  -3.1499939],
         [  7.8500061,   0.8500061,  -6.8999939]]]],
      shape=(365, 60, 60, 3))
Coordinates:
  * dayofyear  (dayofyear) int64 3kB 1 2 3 4 5 6 7 ... 360 361 362 363 364 365
  * latitude   (latitude) float64 480B 60.0 59.75 59.5 ... 45.75 45.5 45.25
  * longitude  (longitude) float64 480B 0.0 0.25 0.5 0.75 ... 14.25 14.5 14.75
  * quantile   (quantile) float64 24B 0.9 0.5 0.1
Attributes: (12/32)
    GRIB_paramId:                             167
    GRIB_dataType:                            an
    GRIB_numberOfPoints:                      1038240
    GRIB_typeOfLevel:                         surface
    GRIB_stepUnits:                           1
    GRIB_stepType:                            instant
    ...                                       ...
    GRIB_totalNumber:                         0
    GRIB_units:                               K
    long_name:                                surface temperature
    units:                                    °C
    standard_name:                            unknown
    last_ingested_valid_time_index:           760439
[15]:

from matplotlib import pyplot as plt

_, ax = plt.subplots()

t2m_point_doy_quantile = t2m_aoi_toi_doy_quantile.sel(**point)

ax.fill_between(
    t2m_point_doy_quantile.dayofyear,
    t2m_point_doy_quantile.sel(quantile=0.9),
    t2m_point_doy_quantile.sel(quantile=0.1),
    alpha=0.3,
)

t2m_point_doy_quantile.sel(quantile=0.5).plot(ax=ax, alpha=0.5)
t2m_point_doy_mean.plot(ax=ax)

[15]:

[<matplotlib.lines.Line2D at 0x133f3c440>]
../../_images/notebooks_edh-learning_10-use-case-era5-small-area_27_1.png