A1 – Preparing Reference Data for Verification#

This annex shows how to download the two reference datasets used to benchmark AIFS forecasts in the user-guide notebooks:

  1. MeteoSwiss automatic-station observations — hourly 2 m temperature and 10 m wind speed, downloaded via meteora and stored in the stationbench format.

  2. ERA5 reanalysis — 2 m temperature and 10 m wind speed from the Copernicus CDS API, also stored in stationbench format.

Both outputs are saved as NetCDF files and loaded lazily with xarray in the main notebooks.

Prerequisites: a valid ~/.cdsapirc file for ERA5 downloads (see the CDS API documentation).

[ ]:

import datetime as dt
import pathlib
import tempfile

import cdsapi
import numpy as np
import xarray as xr
from meteora import clients, units, utils
from shapely import geometry

Configuration#

Adjust the parameters below to match your experiment.

[ ]:

# forecast window
start_date = dt.datetime(2025, 6, 20, 0, tzinfo=dt.UTC)
lead_time = 240  # hours
end_date = start_date + dt.timedelta(hours=lead_time)

# spatial bounding box (lat/lon)
lat_slice = (44.5, 48.5)
lon_slice = (4.5, 11.5)

# output directory
data_dir = pathlib.Path("../data/stationbench")
data_dir.mkdir(parents=True, exist_ok=True)

# output file paths
experiment = "heatwave-2025-jun-ens"
stations_filepath = data_dir / f"{experiment}-stations.nc"
era5_filepath = data_dir / f"{experiment}-era5.nc"

1. MeteoSwiss station observations#

We use meteora to download hourly 2 m temperature and 10 m wind speed from MeteoSwiss automatic stations within the bounding box defined above. The observations are then converted from SI units (K, m/s) and written to a stationbench-ready NetCDF file with dimensions (time, station_id).

[ ]:

if not stations_filepath.exists():
    client = clients.MeteoSwissClient(
        region=geometry.box(*lon_slice, *lat_slice), crs="epsg:4326"
    )
    ts_df = client.get_ts_df(
        variables=["temperature", "wind_speed"],
        start=start_date.replace(tzinfo=None),
        end=end_date.replace(tzinfo=None),
    )
    ts_df = units.convert_units(ts_df, {"temperature": "K"})
    ts_ds = utils.long_to_stationbench(ts_df, client.stations_gdf)
    ts_ds.to_netcdf(stations_filepath)
    print(f"Saved {stations_filepath}")
else:
    print(f"{stations_filepath} already exists; skipping")

ts_ds = xr.open_dataset(stations_filepath)
ts_ds

Saved ../data/stationbench/heatwave-2025-jun-ens-stations.nc

<xarray.Dataset> Size: 4MB
Dimensions:         (time: 1441, station_id: 157)
Coordinates:
  * time            (time) datetime64[ns] 12kB 2025-06-20 ... 2025-06-30
  * station_id      (station_id) <U3 2kB 'ABO' 'AEG' 'AIG' ... 'WFJ' 'WYN' 'ZER'
    longitude       (station_id) float64 1kB ...
    latitude        (station_id) float64 1kB ...
Data variables:
    2m_temperature  (time, station_id) float64 2MB ...
    10m_wind_speed  (time, station_id) float64 2MB ...

2. ERA5 reanalysis#

We download 2 m temperature and 10 m wind components from the ERA5 single-level reanalysis via the Copernicus CDS API for every 6-hourly valid time within the forecast window. The 10 m wind speed scalar is computed from the U and V components, and the dataset is stored in the stationbench format with dimensions (time, prediction_timedelta, latitude, longitude).

[ ]:

if not era5_filepath.exists():
    start_date_naive = start_date.replace(tzinfo=None)
    init_time = np.datetime64(start_date_naive, "ns")
    valid_dts = [
        start_date_naive + dt.timedelta(hours=h) for h in range(6, lead_time + 1, 6)
    ]
    valid_times = np.array([np.datetime64(vdt, "ns") for vdt in valid_dts])
    lead_times = valid_times - init_time

    years = sorted({vdt.strftime("%Y") for vdt in valid_dts})
    months = sorted({vdt.strftime("%m") for vdt in valid_dts})
    days = sorted({vdt.strftime("%d") for vdt in valid_dts})
    hours = sorted({vdt.strftime("%H:00") for vdt in valid_dts})

    with tempfile.TemporaryDirectory() as tmp_dir:
        tmp_filepath = pathlib.Path(tmp_dir) / "era5_raw.nc"
        cdsapi.Client().retrieve(
            "reanalysis-era5-single-levels",
            {
                "product_type": ["reanalysis"],
                "variable": [
                    "2m_temperature",
                    "10m_u_component_of_wind",
                    "10m_v_component_of_wind",
                ],
                "year": years,
                "month": months,
                "day": days,
                "time": hours,
                "data_format": "netcdf",
                "download_format": "unarchived",
            },
            str(tmp_filepath),
        )

        ds = xr.open_dataset(tmp_filepath)
        time_dim = "valid_time" if "valid_time" in ds.dims else "time"
        ds = ds.sel({time_dim: valid_times})

        # rename to short names if needed
        rename_map = {}
        if "2m_temperature" in ds:
            rename_map["2m_temperature"] = "t2m"
        if "10m_u_component_of_wind" in ds:
            rename_map["10m_u_component_of_wind"] = "u10"
        if "10m_v_component_of_wind" in ds:
            rename_map["10m_v_component_of_wind"] = "v10"
        if rename_map:
            ds = ds.rename(rename_map)

        era5_ds = xr.Dataset(
            {
                "2t": (
                    ["time", "prediction_timedelta", "latitude", "longitude"],
                    ds["t2m"].values[np.newaxis],
                ),
                "10si": (
                    ["time", "prediction_timedelta", "latitude", "longitude"],
                    np.sqrt(ds["u10"].values ** 2 + ds["v10"].values ** 2)[np.newaxis],
                ),
            },
            coords={
                "time": [init_time],
                "prediction_timedelta": lead_times,
                "latitude": ds["latitude"].values,
                "longitude": ds["longitude"].values,
            },
        )
        era5_ds.to_netcdf(era5_filepath)
    print(f"Saved {era5_filepath}")
else:
    print(f"{era5_filepath} already exists; skipping")

era5_ds = xr.open_dataset(era5_filepath)
era5_ds

2026-03-18 11:58:39,011 INFO [2025-12-11T00:00:00] Please note that a dedicated catalogue entry for this dataset, post-processed and stored in Analysis Ready Cloud Optimized (ARCO) format (Zarr), is available for optimised time-series retrievals (i.e. for retrieving data from selected variables for a single point over an extended period of time in an efficient way). You can discover it [here](https://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels-timeseries?tab=overview)
2026-03-18 11:58:39,015 INFO Request ID is be7d0364-571e-48b8-a97b-6851aa88506a
2026-03-18 11:58:39,065 INFO status has been updated to accepted
2026-03-18 11:58:47,382 INFO status has been updated to running
2026-03-18 11:59:11,693 INFO status has been updated to successful

Saved ../data/stationbench/heatwave-2025-jun-ens-era5.nc

<xarray.Dataset> Size: 332MB
Dimensions:               (time: 1, prediction_timedelta: 40, latitude: 721,
                           longitude: 1440)
Coordinates:
  * time                  (time) datetime64[ns] 8B 2025-06-20
  * prediction_timedelta  (prediction_timedelta) timedelta64[ns] 320B 06:00:0...
  * latitude              (latitude) float64 6kB 90.0 89.75 ... -89.75 -90.0
  * longitude             (longitude) float64 12kB 0.0 0.25 0.5 ... 359.5 359.8
Data variables:
    2t                    (time, prediction_timedelta, latitude, longitude) float32 166MB ...
    10si                  (time, prediction_timedelta, latitude, longitude) float32 166MB ...