Atmosphere — pressure-level temperature & jet-stream wind¶

Pull mid-tropospheric temperature and wind for a North Atlantic / European box at the 500 hPa pressure level — the classical altitude for synoptic charts and mid-latitude jet-stream analysis. This notebook exercises the pressure-level branch of the package: ERA5 pressure-levels NetCDFs are 4-D (time, level, lat, lon), so you pass level= on AggregationConfig and the aggregator pins the level via pyramids.netcdf.NetCDF.sel before reducing.

Domain context. 500 hPa lies near the steering level for surface weather; jet-stream cores typically sit between 250 and 300 hPa with the strongest winds; 500 hPa temperature is a clean proxy for whether an air mass is anomalously cold or warm aloft. Pulling temperature and the two horizontal wind components at 500 hPa lets you reconstruct the geostrophic flow over a region.

Setup¶

Consolidate the imports and create the output directories. earthlens provides the unified EarthLens entry point plus the ECMWF Catalog / Variable models and the aggregate_netcdf aggregator; pyramids provides Dataset for reading the GeoTIFFs back in for plotting. OUT holds the raw ERA5 NetCDFs and AGG the aggregated GeoTIFFs.

In [ ]:

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
from pyramids.dataset import Dataset

from earthlens import AggregationConfig, EarthLens, aggregate_netcdf
from earthlens.ecmwf import Catalog, Variable

OUT = Path("data/era5-500hpa")
OUT.mkdir(parents=True, exist_ok=True)

AGG = OUT / "aggregated"
AGG.mkdir(parents=True, exist_ok=True)

from pathlib import Path import matplotlib.pyplot as plt import numpy as np from pyramids.dataset import Dataset from earthlens import AggregationConfig, EarthLens, aggregate_netcdf from earthlens.ecmwf import Catalog, Variable OUT = Path("data/era5-500hpa") OUT.mkdir(parents=True, exist_ok=True) AGG = OUT / "aggregated" AGG.mkdir(parents=True, exist_ok=True)

Step 1 — pressure-level catalog entries¶

Pressure-level variables live on reanalysis-era5-pressure-levels. Each row carries cds_pressure_level so the request automatically forwards a pressure_level field; the catalog's per-row default is ['1000'] for the bundled ERA5 entries.

In [ ]:

cat = Catalog()
for code in ("temperature", "u-component-of-wind", "v-component-of-wind"):
    spec = cat.get_variable("reanalysis-era5-pressure-levels", code)
    print(
        f"{code:25s}  nc={spec.nc_variable:5s}  units={spec.units:5s}  "
        f"levels={spec.cds_pressure_level}  is_flux={spec.is_flux}"
    )

cat = Catalog() for code in ("temperature", "u-component-of-wind", "v-component-of-wind"): spec = cat.get_variable("reanalysis-era5-pressure-levels", code) print( f"{code:25s} nc={spec.nc_variable:5s} units={spec.units:5s} " f"levels={spec.cds_pressure_level} is_flux={spec.is_flux}" )

Step 2 — daily means at 500 hPa over Europe¶

One week of daily aggregates over the North Atlantic / Europe box (30°–65°N, 50°W–30°E). The catalog's bundled pressure-level rows ship with cds_pressure_level=['1000']; we override at the catalog row by constructing custom variable specs in code if we needed to. For this demo we'll lean on the request-shape default (which submits the level from the catalog row), then slice to 500 hPa at aggregation time via level=500 on AggregationConfig.

Reality check: the catalog default of 1000 hPa means the retrieved NetCDF only contains that single level. To get 500 hPa, pass extras={"pressure_level": ["500"]} directly on the catalog row by monkey-patching, OR — more simply for this demo — write a custom Variable and call aggregate_netcdf directly. The latter keeps the demo readable.

In [ ]:

# Build per-variable Variable instances with pressure_level=500.
# We patch the Catalog at retrieval time so `_api` sees the
# right cds_pressure_level field.
PATCHED = {
    "temperature": Variable(
        cds_dataset="reanalysis-era5-pressure-levels",
        cds_variable="temperature",
        nc_variable="t",
        units="K",
        product_type=["reanalysis"],
        cds_pressure_level=["500"],
    ),
    "u-component-of-wind": Variable(
        cds_dataset="reanalysis-era5-pressure-levels",
        cds_variable="u_component_of_wind",
        nc_variable="u",
        units="m s**-1",
        product_type=["reanalysis"],
        cds_pressure_level=["500"],
    ),
    "v-component-of-wind": Variable(
        cds_dataset="reanalysis-era5-pressure-levels",
        cds_variable="v_component_of_wind",
        nc_variable="v",
        units="m s**-1",
        product_type=["reanalysis"],
        cds_pressure_level=["500"],
    ),
}
PATCHED  # one Variable per code, 500 hPa pinned

# Build per-variable Variable instances with pressure_level=500. # We patch the Catalog at retrieval time so `_api` sees the # right cds_pressure_level field. PATCHED = { "temperature": Variable( cds_dataset="reanalysis-era5-pressure-levels", cds_variable="temperature", nc_variable="t", units="K", product_type=["reanalysis"], cds_pressure_level=["500"], ), "u-component-of-wind": Variable( cds_dataset="reanalysis-era5-pressure-levels", cds_variable="u_component_of_wind", nc_variable="u", units="m s**-1", product_type=["reanalysis"], cds_pressure_level=["500"], ), "v-component-of-wind": Variable( cds_dataset="reanalysis-era5-pressure-levels", cds_variable="v_component_of_wind", nc_variable="v", units="m s**-1", product_type=["reanalysis"], cds_pressure_level=["500"], ), } PATCHED # one Variable per code, 500 hPa pinned

Step 3 — submit the retrieve¶

We use ECMWF._api(var_info) directly (the lower-level entry point) so the patched Variable instances drive the request without going through the catalog. For multi-variable retrieval through the facade with the default catalog values, see cds_quickstart.ipynb.

Build the request¶

Construct the EarthLens ECMWF backend once for the whole week and box; the patched Variable instances (not the catalog defaults) will drive each retrieve.

In [ ]:

ecmwf = EarthLens(
    data_source="ecmwf",
    start="2022-01-01",
    end="2022-01-07",
    cadence="daily",
    dataset="reanalysis-era5-pressure-levels",
    variables=list(PATCHED),
    aoi=[-50.0, 30.0, 30.0, 65.0],
    path=str(OUT),
)

ecmwf = EarthLens( data_source="ecmwf", start="2022-01-01", end="2022-01-07", cadence="daily", dataset="reanalysis-era5-pressure-levels", variables=list(PATCHED), aoi=[-50.0, 30.0, 30.0, 65.0], path=str(OUT), )

Submit each retrieve¶

Call ecmwf._api(var_info) once per patched Variable to fetch the three 500 hPa fields. Each request goes to the CDS queue and may take a few minutes.

In [ ]:

for code, var_info in PATCHED.items():
    print(f"Retrieving {code} at 500 hPa...")
    ecmwf._api(var_info)
print("done")

for code, var_info in PATCHED.items(): print(f"Retrieving {code} at 500 hPa...") ecmwf._api(var_info) print("done")

Step 4 — aggregate to a daily mean with level=500¶

The retrieved NetCDF carries shape (time, level=1, lat, lon). The aggregator detects the pressure-level dim and, when level=500 is set, calls nc.sel(pressure_level=500) to pin the slice before reducing. Because we only requested level 500, this is a 1-element pin — harmless. (For multi-level retrievals, this is the knob that lets you aggregate one level at a time without splitting the input file.)

In [ ]:

for code, var_info in PATCHED.items():
    nc_path = OUT / f"{var_info.cds_variable}_reanalysis-era5-pressure-levels.nc"
    aggregate_netcdf(
        nc_path,
        var_info,
        AggregationConfig(
            freq="1D",
            op="auto",
            out_dir=AGG,
            level=500,
        ),
    )

sorted(p.name for p in AGG.glob("*.tif"))

for code, var_info in PATCHED.items(): nc_path = OUT / f"{var_info.cds_variable}_reanalysis-era5-pressure-levels.nc" aggregate_netcdf( nc_path, var_info, AggregationConfig( freq="1D", op="auto", out_dir=AGG, level=500, ), ) sorted(p.name for p in AGG.glob("*.tif"))

Step 5 — plot 500 hPa temperature with wind vectors¶

Pick the first day of the week. Overlay arrows for the (u, v) wind on top of the temperature field — the canonical synoptic chart shape.

Read the first day's fields¶

Load the earliest temperature, u, and v GeoTIFFs into arrays. The read_file(...) / read_array() chain is split so each Dataset is read on its own line.

In [ ]:

t_path = str(min(AGG.glob("temperature_1D_*.tif")))
u_path = str(min(AGG.glob("u_component_of_wind_1D_*.tif")))
v_path = str(min(AGG.glob("v_component_of_wind_1D_*.tif")))

t_air = Dataset.read_file(t_path).read_array()
u_wind = Dataset.read_file(u_path).read_array()
v_wind = Dataset.read_file(v_path).read_array()

t_path = str(min(AGG.glob("temperature_1D_*.tif"))) u_path = str(min(AGG.glob("u_component_of_wind_1D_*.tif"))) v_path = str(min(AGG.glob("v_component_of_wind_1D_*.tif"))) t_air = Dataset.read_file(t_path).read_array() u_wind = Dataset.read_file(u_path).read_array() v_wind = Dataset.read_file(v_path).read_array()

Plot temperature with wind vectors¶

Render the 500 hPa temperature field in °C and overlay subsampled (u, v) arrows on top — the canonical synoptic chart shape.

In [ ]:

fig, ax = plt.subplots(figsize=(10, 5))
im = ax.imshow(t_air - 273.15, cmap="RdBu_r", origin="upper", aspect="auto")
fig.colorbar(im, ax=ax, label="500 hPa T [°C]")

# Subsample for arrow density.
step = 8
yy, xx = np.mgrid[0 : t_air.shape[0] : step, 0 : t_air.shape[1] : step]
ax.quiver(
    xx, yy, u_wind[::step, ::step], -v_wind[::step, ::step], color="black", scale=600
)
ax.set_title("ERA5 500 hPa temperature + wind, 2022-01-01 — N. Atlantic / Europe")
ax.set_xlabel("lon (pixels)")
ax.set_ylabel("lat (pixels)")
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(10, 5)) im = ax.imshow(t_air - 273.15, cmap="RdBu_r", origin="upper", aspect="auto") fig.colorbar(im, ax=ax, label="500 hPa T [°C]") # Subsample for arrow density. step = 8 yy, xx = np.mgrid[0 : t_air.shape[0] : step, 0 : t_air.shape[1] : step] ax.quiver( xx, yy, u_wind[::step, ::step], -v_wind[::step, ::step], color="black", scale=600 ) ax.set_title("ERA5 500 hPa temperature + wind, 2022-01-01 — N. Atlantic / Europe") ax.set_xlabel("lon (pixels)") ax.set_ylabel("lat (pixels)") plt.tight_layout() plt.show()

Notes¶

• Why we built Variable by hand. The bundled catalog rows ship cds_pressure_level=['1000'] (the default level the catalog was audited at). For a different level you either edit the YAML row in-place, monkey-patch the catalog at runtime, or — as here — construct Variable instances directly. The first option is best for recurring workflows; the third is fine for a one-off notebook.
• level=500 on AggregationConfig. This is the production pressure-level path. The aggregator reads nc.dimension_names, detects pressure_level, and calls nc.sel(pressure_level=500) before reducing. If you forget to pass level=, the aggregator raises a ValueError naming the dim and pointing at the AggregationConfig.level field.
• Multi-level retrievals. Build a Variable with cds_pressure_level=['200', '500', '850'], retrieve once, then call aggregate_netcdf three times with different level= values to produce per-level GeoTIFF stacks without re-downloading.
• Wind-direction sign convention. ERA5 wind components are earth-relative; positive u is eastward, positive v is northward. In imshow(origin='upper') the vertical axis points downward, so we pass -V to quiver to keep arrows pointing in the right direction on screen.