ERA5 data on the HEALPix grid

ERA5 data on the HEALPix grid#

ERA51 data for selected variables has been restructured and remapped onto the HEALPix grid for the years 2010 to 2023. This dataset, which we call HERA5, is accessible through the internal TCO catalog and can be read and used by everyone. HERA5 provides data at different temporal resolutions and enables an open and fast access to the data from everywhere and without a local copy. This notebook shows the access to the data and some basic usage examples.

This tutorial makes use of the easygems python package which hosts several convenience functions to work with and display Earth-system data on the HEALPix grid. If you don’t have the package installed, you can run the code cell below to install it via pip:

%pip install easygems

Loading the data from the catalog#

First, let’s import all packages that are required to run the full notebook.

import intake
import healpy as hp
import cmocean
import numpy as np
import xarray as xr

import cartopy.crs as ccrs
import cartopy.feature as cf
import matplotlib.pyplot as plt

import easygems.healpix as egh

Matplotlib is building the font cache; this may take a moment.

Now we can open the catalog and select the hierarchical HEALPix ERA5 data HERA5.

cat = intake.open_catalog("https://tcodata.mpimet.mpg.de/internal.yaml")
list(cat)

['BCO', 'METEOR', 'ORCESTRA', 'HERA5', 'HIFS', 'airinfo']

cat.HERA5.to_dask().pipe(egh.attach_coords)

/usr/share/miniconda3/envs/orcestra_book/lib/python3.12/site-packages/intake_xarray/base.py:21: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  'dims': dict(self._ds.dims),

<xarray.Dataset> Size: 66GB
Dimensions:    (time: 164, cell: 196608, level: 29)
Coordinates:
    latitude   (cell) float32 786kB dask.array<chunksize=(196608,), meta=np.ndarray>
  * level      (level) int64 232B 50 70 100 125 150 175 ... 900 925 950 975 1000
    longitude  (cell) float32 786kB dask.array<chunksize=(196608,), meta=np.ndarray>
  * time       (time) datetime64[ns] 1kB 2010-01-16T04:15:00 ... 2023-08-16T0...
    crs        int64 8B 0
    lat        (cell) float64 2MB 0.2984 0.5968 0.5968 ... -0.5968 -0.2984
    lon        (cell) float64 2MB 45.0 45.35 44.65 45.0 ... 315.4 314.6 315.0
Dimensions without coordinates: cell
Data variables: (12/73)
    100u       (time, cell) float32 129MB dask.array<chunksize=(24, 4096), meta=np.ndarray>
    100v       (time, cell) float32 129MB dask.array<chunksize=(24, 4096), meta=np.ndarray>
    10si       (time, cell) float32 129MB dask.array<chunksize=(24, 4096), meta=np.ndarray>
    10u        (time, cell) float32 129MB dask.array<chunksize=(24, 4096), meta=np.ndarray>
    10v        (time, cell) float32 129MB dask.array<chunksize=(24, 4096), meta=np.ndarray>
    2d         (time, cell) float32 129MB dask.array<chunksize=(24, 4096), meta=np.ndarray>
    ...         ...
    isor       (cell) float32 786kB dask.array<chunksize=(196608,), meta=np.ndarray>
    lsm        (cell) float32 786kB dask.array<chunksize=(196608,), meta=np.ndarray>
    sdfor      (cell) float32 786kB dask.array<chunksize=(196608,), meta=np.ndarray>
    sdor       (cell) float32 786kB dask.array<chunksize=(196608,), meta=np.ndarray>
    slor       (cell) float32 786kB dask.array<chunksize=(196608,), meta=np.ndarray>
    z_sfc      (cell) float32 786kB dask.array<chunksize=(196608,), meta=np.ndarray>
Attributes:
    acknowledgment:  Contains modified Copernicus Climate Change Service info...
    contact:         lukas.kluft@mpimet.mpg.de
    creator:         Lukas Kluft
    description:     Selected variables from ERA5, restructured and saved on ...
    institution:     Max Planck Institute for Meteorology
    source:          Post-processed dataset based on the ERA5 mirror located ...
    title:           HERA5 - HEALPixelation of ERA5

Some notes on the dataset#

  • Dimensions

    • the dimensions are time, cell and (pressure-)level and the coordinate reference system crs. Here, cell refers to the cells in the HEALPix grid, which is a pixelation of a sphere. Get to know more about HEALPix here. The zoom level of the HERA5 data is \(z = 7\) and the number of cells respectively \(12 \cdot z^{4} = 196608\)

    • you can distinguish between 2d and 3d variables by their dimensions, the former having only (time, cell) and the latter having (time, level, cell).

  • Temporal hierarchy

    • The dataset is constructed hierarchically in time, which means that you can directly access either monthly aggregated values P1M, daily P1D, or the original hourly data PT1H. By default, the catalog entry points to the monthly values which are great to do a fast analysis of longterm and coarse features such as calculating the global mean surface temperature over 10 years. If you need higher temporal resolutions, you can add the respective acronym when loading the catalog: ds = cat.HERA5(time="PT1H").to_dask().

    • In principle, one could also implement spatial hierarchies. However, the original ERA5 data is already rather coarse such that there is no need for coarser spatial (zoom) levels in the present case.

  • Spatial resolution

    • The spatial resolution is given by the zoom level \(z = 7\) which is closest to the original ERA5 model output and roughly corresponds to 0.35 x 0.7° effective grid resolution.

  • Variables

    • Most variables are from the ERA5 “analysis” product. Only few variables such as total precipitation and boundary layer height are added from the “forecast” product. Each variable has a type attribute which will let you know whether it is from the analysis or forecast product.

Plotting examples#

The below plotting examples show some typical usecases and the respective healpy functions that are frequently used in the analysis.

2m temperature at the BCO in Aug-Sep 2020#

We would like to plot the 2m air temperature at one location, the BCO, in August and September 2020. For this rather short time period, we would like to work with the hourly data. Therefore, we need to specify the time hierarchy when addressing the data in the catalog.

era5 = cat.HERA5(time="PT1H").to_dask().pipe(egh.attach_coords)

/usr/share/miniconda3/envs/orcestra_book/lib/python3.12/site-packages/intake_xarray/base.py:21: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  'dims': dict(self._ds.dims),

Next, we select the cell index that is the nearest neighbor to the BCO location using healpy.ang2pix().

i_bco = hp.ang2pix(
    egh.get_nside(era5),
    -59.42875, # longitude at BCO
    13.16264, # latitude at BCO
    nest=egh.get_nest(era5),
    lonlat=True,
)

Now, we can plot the time series by selecting the respective cell as well as the time range.

era5["2t"].sel(cell=i_bco, time=slice("2020-08", "2020-09")).plot()

[<matplotlib.lines.Line2D at 0x7f5604ddcb00>]
_images/e3b7d43acbb0c27a50982d6df8b04305d3dff7de9cc6f975be7e3cdfcf4b0ae6.png

Evolution of the moisture field at BCO#

Next, we extend our analysis and use a 3d variable to plot a time-height diagramm of moisture. Again, we select the BCO cell, a time range and now also a range in height level where we leave out the upper most layers (<100hPa) as we are interested in the troposphere only.

era5["q"].sel(
    cell=i_bco,
    time="2020-08",
    level=slice(100, 1000),
).plot(
    x="time",
    yincrease=False,
    cmap="cmo.dense",
    vmin=0,
)

<matplotlib.collections.QuadMesh at 0x7f55ef624770>
_images/e21eeba50d0968dabba53aaced0382daba1a2273383f101297ea25c72d47e462.png

Select ORCESTRA region#

We want to see how the RH profile changes with latitude when averaged over the ORCESTRA campaign region on the 20th of August 2020. The easygems package provides the the isel_extent function which takes a lat/lon extent ([W, E, S, N]) as input and returns the corresponding indices of the global HEALPix grid. These indices can then be used to select the region of interest.

is_orcestra = egh.isel_extent(era5, [-60, -10, -5, 20])

fig, ax = plt.subplots()
era5["r"].sel(
    cell=is_orcestra,
    time="2020-08-20",
    level=slice(100, 1000),
).mean(
    "time"
).groupby(
    era5.lat.sel(cell=is_orcestra)
).mean().plot(
    x="lat",
    yincrease=False,
    cmap="cmo.dense",
    ax=ax,
)
ax.set_title("Relative humidity profile on 2020-08-20 (ORCESTRA region)")

Text(0.5, 1.0, 'Relative humidity profile on 2020-08-20 (ORCESTRA region)')
_images/28d5da7f2ec13f1e4929d60f4c004caa8c97cf3efe96be394ad328393d4afba4.png

A snapshot of global SST#

using the healpy.mollview() function

hp.mollview(
    era5["sst"].sel(time="2020-08-01T12:00:00", method="nearest"),
    title="global SST on 2020-01-08 12:00",
    nest=egh.get_nest(era5),
    flip="geo",
    min=273,
    max=303,
    cmap="cmo.thermal",
    unit=era5["sst"].units,
)

/usr/share/miniconda3/envs/orcestra_book/lib/python3.12/site-packages/xarray/core/indexing.py:1620: PerformanceWarning: Slicing with an out-of-order index is generating 107 times more chunks
  return self.array[key]
_images/286a95040d0044bd4f5a6979a92dc722b1d1756567dd3828f15e1f23de2e964b.png

Plotting the Atlantic basin only#

The healpy package only provides functions to plot the whole globe. However, in real-life analyses it is often necessary to constrain a plot to certain regions. The function easygems.healpix_show() allows to plot two-dimensional data onto a cartopy GeoAxis with arbitrary projection and extent.

era5 = cat.HERA5(time="P1M").to_dask().pipe(egh.attach_coords)
var = era5["tp"]

fig, ax = plt.subplots(figsize=(10, 6), subplot_kw={"projection": ccrs.PlateCarree()})
ax.set_extent([-65, -5, -10, 25])
ax.coastlines(lw=0.8)
im = egh.healpix_show(
    var.sel(time="2020-08").values,
    method="linear",
    cmap="cmo.rain",
    vmin=0,
)
fig.colorbar(im, label=f"{var.long_name} / {var.units}", shrink=0.7)

/usr/share/miniconda3/envs/orcestra_book/lib/python3.12/site-packages/intake_xarray/base.py:21: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  'dims': dict(self._ds.dims),

<matplotlib.colorbar.Colorbar at 0x7f55a1eb5e80>

/usr/share/miniconda3/envs/orcestra_book/lib/python3.12/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)
_images/930d0a1797d6bef624cc88167563a94abff48e704a11b09a74a79ab4dbbdc286.png

Further reading#

1

Contains modified Copernicus Climate Change Service information 2020. Neither the European Commission nor ECMWF is responsible for any use that may be made of the Copernicus information or data it contains.

Contents