Projection-based weather regimes

docs/examples/index.rst

@@ -10,3 +10,4 @@ Here is a list of example notebooks to illustrate how to use earthkit-data.
     :maxdepth: 1
 
     return_period.ipynb
+    seven_weather_regimes.ipynb

docs/examples/seven_weather_regimes.ipynb

@@ -0,0 +1,330 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "9485a729",
+   "metadata": {},
+   "source": [
+    "# Year-round North Atlantic-European Weather Regimes\n",
+    "\n",
+    "Following the example in Part 4 of [github.com/cmgrams/wr_data_package_era5](https://github.com/cmgrams/wr_data_package_era5/blob/main/wr_data_package_V1.0/scripts_first_steps/WR_read_example.ipynb), compute the regime index for the seven-regime classification of [Grams (2026, in review)](https://doi.org/10.5194/egusphere-2025-6385)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "38c63573",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pooch\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import xarray as xr\n",
+    "\n",
+    "from earthkit.meteo import regimes"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1fb3f386",
+   "metadata": {},
+   "source": [
+    "## Get the regime classification data\n",
+    "\n",
+    "Retrieve data from https://zenodo.org/records/18154492."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "94c5d0a6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "files = pooch.retrieve(\n",
+    "    url=\"doi:10.5281/zenodo.18154492/wr_data_package_V1.1.zip\",\n",
+    "    known_hash=\"dc942ff2a1b3da6dedd3b0b2fadda017fff9e8fc10228ace31c2209a9be7dc62\",\n",
+    "    processor=pooch.Unzip()\n",
+    ")\n",
+    "\n",
+    "def get_file(name):\n",
+    "    for file in files:\n",
+    "        if file.endswith(name):\n",
+    "            return file\n",
+    "    raise FileNotFoundError(name)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0c18029b",
+   "metadata": {},
+   "source": [
+    "## Regime pattern normalisation based on day-of-year\n",
+    "\n",
+    "Load the file with normalisation weights from the repository."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "f4fdd30f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mod_ds = xr.open_dataset(get_file(\"wr_data/EOFs_WRs.nc\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b73380a3",
+   "metadata": {},
+   "source": [
+    "Create a lookup table for the normalisation weights and define a function that finds the appropriate weight for a given date."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "8e820f2d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mod = 1. / mod_ds[\"normwgt\"].to_series()\n",
+    "mod.index = pd.Index(zip(mod.index.month, mod.index.day, mod.index.hour))\n",
+    "\n",
+    "def pattern_normalisation_weight(date):\n",
+    "    date = np.asarray(date)\n",
+    "    shp = date.shape\n",
+    "    date = pd.to_datetime(date.flatten())\n",
+    "    if isinstance(date, pd.Timestamp):\n",
+    "        return mod.loc[(date.month, date.day, date.hour)]\n",
+    "    else:\n",
+    "        idx = list(zip(date.month, date.day, date.hour))\n",
+    "        return mod.loc[idx].values.reshape(shp)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "651d73e0",
+   "metadata": {},
+   "source": [
+    "## Load the regime patterns\n",
+    "\n",
+    "Patterns are stored in a NetCDF file in the repository."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "418019a3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pattern_ds = xr.open_dataset(get_file(\"wr_data/Normed_Z0500-patterns_EOFdomain.nc\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fdba08fc",
+   "metadata": {},
+   "source": [
+    "Seven base patterns. The normalisation of projections is implemented via a modulation of the regime patterns, i.e., we vary the amplitude of the patterns."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "14ab5c58",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "patterns = regimes.ModulatedRegimePatterns(\n",
+    "    regimes=pattern_ds.attrs[\"ClassNames\"].split(),\n",
+    "    grid={\n",
+    "        \"grid\": [0.5, 0.5],\n",
+    "        \"area\": [90, -80, 30, 40]  # 30-90°N, 80°W-40°E\n",
+    "    },\n",
+    "    patterns=pattern_ds[\"Z0500_mean\"].values,\n",
+    "    modulator=pattern_normalisation_weight\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "69a9f27d",
+   "metadata": {},
+   "source": [
+    "**Note:** the regime data from this repository uses ascending order for the latitude coordinate. The latitude ordering of fields projected onto these patterns must match."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "18c7b6ca",
+   "metadata": {},
+   "source": [
+    "## Project the test field onto the patterns"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8a2a0526",
+   "metadata": {},
+   "source": [
+    "Load the test field provided with the dataset. This field contains pre-processed Z500 anomalies ready for projection."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "f0b8eeba",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_field = xr.open_dataset(get_file(\"wr_data_package_V1.1/example_data/Z0500_20250601_00.nc\")).squeeze()\n",
+    "\n",
+    "# Extract values for the EUR-ATL region and create associated coordinates\n",
+    "# following the example in the reference notebook\n",
+    "field = ds_field[\"Z0\"].values[240:361,200:441]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a491d706",
+   "metadata": {},
+   "source": [
+    "Create area-based weights for the grid points to use in the projection."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "6b0a13c0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lat = np.linspace(ds_field.attrs[\"domymin\"], ds_field.attrs[\"domymax\"], ds_field[\"Z0\"].shape[0])\n",
+    "\n",
+    "weights = np.cos(np.deg2rad(lat))\n",
+    "weights = np.repeat(weights[240:361], field.shape[1]).reshape(field.shape)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ea43918a",
+   "metadata": {},
+   "source": [
+    "Project onto the regime patterns, supply the valid date of the field for the modulator function to select the normalisation weight."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "d016fba8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "projections = regimes.project(field, patterns, weights, date=np.datetime64(\"2025-06-01 00:00\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9b316d2c",
+   "metadata": {},
+   "source": [
+    "## Compute the regime index\n",
+    "\n",
+    "Standardise the projections to obtain the regime index.\n",
+    "Reference values are given in the notebook:\n",
+    "\n",
+    "    weather regime indices for  20250601_00\n",
+    "    0.8218320408050166 AT\n",
+    "    1.1691867614026132 ZO\n",
+    "    1.0489664646400954 ScTr\n",
+    "    -0.6948883613466961 AR\n",
+    "    -0.5738999234610754 EuBL\n",
+    "    -0.9323144170607468 ScBL\n",
+    "    -0.6369122243737739 GL\n",
+    "\n",
+    "Read the standardisation parameters (mean and standard deviation) from the text file in the repo."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "id": "342bdb6e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def read_wri_std_parameters(path):\n",
+    "    f.readline()\n",
+    "    name = f.readline().strip().split()\n",
+    "    mean = f.readline().strip().split()[1:]\n",
+    "    std = f.readline().strip().split()[1:]\n",
+    "    return (\n",
+    "        {n: float(v) for n, v in zip(name, mean)},\n",
+    "        {n: float(v) for n, v in zip(name, std)}\n",
+    "    )\n",
+    "    \n",
+    "with open(get_file(\"wr_data/WRI_std_params.txt\"), \"r\") as f:\n",
+    "    norm_mean, norm_std = read_wri_std_parameters(f)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9485431d",
+   "metadata": {},
+   "source": [
+    "Compute the regime index."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "913f8d12",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'AT': np.float32(0.8218321),\n",
+       " 'ZO': np.float32(1.169187),\n",
+       " 'ScTr': np.float32(1.0489665),\n",
+       " 'AR': np.float32(-0.6948884),\n",
+       " 'EuBL': np.float32(-0.5738999),\n",
+       " 'ScBL': np.float32(-0.9323146),\n",
+       " 'GL': np.float32(-0.6369122)}"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "regimes.standardise(projections, norm_mean, norm_std)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "earthkit-meteo",
+   "language": "python",
+   "name": "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.14.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

src/earthkit/meteo/regimes/__init__.py

@@ -0,0 +1,13 @@
+# (C) Copyright 2021 ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+
+from .index import project
+from .index import standardise
+from .patterns import ConstantRegimePatterns
+from .patterns import ModulatedRegimePatterns
+from .patterns import RegimePatterns

src/earthkit/meteo/regimes/array/__init__.py

@@ -0,0 +1,9 @@
+# (C) Copyright 2021 ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+
+from .index import *