{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# A1 – Preparing Reference Data for Verification\n", "\n", "This annex shows how to download the two reference datasets used to benchmark AIFS forecasts in the user-guide notebooks:\n", "\n", "1. **MeteoSwiss automatic-station observations** — hourly 2 m temperature and 10 m wind speed, downloaded via [meteora](https://github.com/martibosch/meteora) and stored in the [stationbench](https://github.com/martibosch/stationbench) format.\n", "2. **ERA5 reanalysis** — 2 m temperature and 10 m wind speed from the [Copernicus CDS API](https://cds.climate.copernicus.eu), also stored in stationbench format.\n", "\n", "Both outputs are saved as NetCDF files and loaded lazily with xarray in the main notebooks.\n", "\n", "> **Prerequisites**: a valid `~/.cdsapirc` file for ERA5 downloads (see the [CDS API documentation](https://cds.climate.copernicus.eu/how-to-api))." ] }, { "cell_type": "code", "execution_count": null, "id": "1", "metadata": {}, "outputs": [], "source": [ "import datetime as dt\n", "import pathlib\n", "import tempfile\n", "\n", "import cdsapi\n", "import numpy as np\n", "import xarray as xr\n", "from meteora import clients, units, utils\n", "from shapely import geometry" ] }, { "cell_type": "markdown", "id": "2", "metadata": {}, "source": [ "## Configuration\n", "\n", "Adjust the parameters below to match your experiment." ] }, { "cell_type": "code", "execution_count": null, "id": "3", "metadata": {}, "outputs": [], "source": [ "# forecast window\n", "start_date = dt.datetime(2025, 6, 20, 0, tzinfo=dt.UTC)\n", "lead_time = 240 # hours\n", "end_date = start_date + dt.timedelta(hours=lead_time)\n", "\n", "# spatial bounding box (lat/lon)\n", "lat_slice = (44.5, 48.5)\n", "lon_slice = (4.5, 11.5)\n", "\n", "# output directory\n", "data_dir = pathlib.Path(\"../data/stationbench\")\n", "data_dir.mkdir(parents=True, exist_ok=True)\n", "\n", "# output file paths\n", "experiment = \"heatwave-2025-jun-ens\"\n", "stations_filepath = data_dir / f\"{experiment}-stations.nc\"\n", "era5_filepath = data_dir / f\"{experiment}-era5.nc\"" ] }, { "cell_type": "markdown", "id": "4", "metadata": {}, "source": [ "## 1. MeteoSwiss station observations\n", "\n", "We use [meteora](https://github.com/martibosch/meteora) to download hourly 2 m temperature and 10 m wind speed from MeteoSwiss automatic stations within the bounding box defined above.\n", "The observations are then converted from SI units (K, m/s) and written to a stationbench-ready NetCDF file with dimensions `(time, station_id)`." ] }, { "cell_type": "code", "execution_count": null, "id": "5", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "c62d9124a50d4c6090e651431011763f", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Stations: 0%| | 0/158 [00:00\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset> Size: 4MB\n",
       "Dimensions:         (time: 1441, station_id: 157)\n",
       "Coordinates:\n",
       "  * time            (time) datetime64[ns] 12kB 2025-06-20 ... 2025-06-30\n",
       "  * station_id      (station_id) <U3 2kB 'ABO' 'AEG' 'AIG' ... 'WFJ' 'WYN' 'ZER'\n",
       "    longitude       (station_id) float64 1kB ...\n",
       "    latitude        (station_id) float64 1kB ...\n",
       "Data variables:\n",
       "    2m_temperature  (time, station_id) float64 2MB ...\n",
       "    10m_wind_speed  (time, station_id) float64 2MB ...
" ], "text/plain": [ " Size: 4MB\n", "Dimensions: (time: 1441, station_id: 157)\n", "Coordinates:\n", " * time (time) datetime64[ns] 12kB 2025-06-20 ... 2025-06-30\n", " * station_id (station_id) \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset> Size: 332MB\n",
       "Dimensions:               (time: 1, prediction_timedelta: 40, latitude: 721,\n",
       "                           longitude: 1440)\n",
       "Coordinates:\n",
       "  * time                  (time) datetime64[ns] 8B 2025-06-20\n",
       "  * prediction_timedelta  (prediction_timedelta) timedelta64[ns] 320B 06:00:0...\n",
       "  * latitude              (latitude) float64 6kB 90.0 89.75 ... -89.75 -90.0\n",
       "  * longitude             (longitude) float64 12kB 0.0 0.25 0.5 ... 359.5 359.8\n",
       "Data variables:\n",
       "    2t                    (time, prediction_timedelta, latitude, longitude) float32 166MB ...\n",
       "    10si                  (time, prediction_timedelta, latitude, longitude) float32 166MB ...
" ], "text/plain": [ " Size: 332MB\n", "Dimensions: (time: 1, prediction_timedelta: 40, latitude: 721,\n", " longitude: 1440)\n", "Coordinates:\n", " * time (time) datetime64[ns] 8B 2025-06-20\n", " * prediction_timedelta (prediction_timedelta) timedelta64[ns] 320B 06:00:0...\n", " * latitude (latitude) float64 6kB 90.0 89.75 ... -89.75 -90.0\n", " * longitude (longitude) float64 12kB 0.0 0.25 0.5 ... 359.5 359.8\n", "Data variables:\n", " 2t (time, prediction_timedelta, latitude, longitude) float32 166MB ...\n", " 10si (time, prediction_timedelta, latitude, longitude) float32 166MB ..." ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "if not era5_filepath.exists():\n", " start_date_naive = start_date.replace(tzinfo=None)\n", " init_time = np.datetime64(start_date_naive, \"ns\")\n", " valid_dts = [\n", " start_date_naive + dt.timedelta(hours=h) for h in range(6, lead_time + 1, 6)\n", " ]\n", " valid_times = np.array([np.datetime64(vdt, \"ns\") for vdt in valid_dts])\n", " lead_times = valid_times - init_time\n", "\n", " years = sorted({vdt.strftime(\"%Y\") for vdt in valid_dts})\n", " months = sorted({vdt.strftime(\"%m\") for vdt in valid_dts})\n", " days = sorted({vdt.strftime(\"%d\") for vdt in valid_dts})\n", " hours = sorted({vdt.strftime(\"%H:00\") for vdt in valid_dts})\n", "\n", " with tempfile.TemporaryDirectory() as tmp_dir:\n", " tmp_filepath = pathlib.Path(tmp_dir) / \"era5_raw.nc\"\n", " cdsapi.Client().retrieve(\n", " \"reanalysis-era5-single-levels\",\n", " {\n", " \"product_type\": [\"reanalysis\"],\n", " \"variable\": [\n", " \"2m_temperature\",\n", " \"10m_u_component_of_wind\",\n", " \"10m_v_component_of_wind\",\n", " ],\n", " \"year\": years,\n", " \"month\": months,\n", " \"day\": days,\n", " \"time\": hours,\n", " \"data_format\": \"netcdf\",\n", " \"download_format\": \"unarchived\",\n", " },\n", " str(tmp_filepath),\n", " )\n", "\n", " ds = xr.open_dataset(tmp_filepath)\n", " time_dim = \"valid_time\" if \"valid_time\" in ds.dims else \"time\"\n", " ds = ds.sel({time_dim: valid_times})\n", "\n", " # rename to short names if needed\n", " rename_map = {}\n", " if \"2m_temperature\" in ds:\n", " rename_map[\"2m_temperature\"] = \"t2m\"\n", " if \"10m_u_component_of_wind\" in ds:\n", " rename_map[\"10m_u_component_of_wind\"] = \"u10\"\n", " if \"10m_v_component_of_wind\" in ds:\n", " rename_map[\"10m_v_component_of_wind\"] = \"v10\"\n", " if rename_map:\n", " ds = ds.rename(rename_map)\n", "\n", " era5_ds = xr.Dataset(\n", " {\n", " \"2t\": (\n", " [\"time\", \"prediction_timedelta\", \"latitude\", \"longitude\"],\n", " ds[\"t2m\"].values[np.newaxis],\n", " ),\n", " \"10si\": (\n", " [\"time\", \"prediction_timedelta\", \"latitude\", \"longitude\"],\n", " np.sqrt(ds[\"u10\"].values ** 2 + ds[\"v10\"].values ** 2)[np.newaxis],\n", " ),\n", " },\n", " coords={\n", " \"time\": [init_time],\n", " \"prediction_timedelta\": lead_times,\n", " \"latitude\": ds[\"latitude\"].values,\n", " \"longitude\": ds[\"longitude\"].values,\n", " },\n", " )\n", " era5_ds.to_netcdf(era5_filepath)\n", " print(f\"Saved {era5_filepath}\")\n", "else:\n", " print(f\"{era5_filepath} already exists; skipping\")\n", "\n", "era5_ds = xr.open_dataset(era5_filepath)\n", "era5_ds" ] } ], "metadata": { "kernelspec": { "display_name": "Python (Pixi)", "language": "python", "name": "pixi-kernel-python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.12" }, "pixi-kernel": { "environment": "user-guide" } }, "nbformat": 4, "nbformat_minor": 5 }