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