Year-round North Atlantic-European Weather Regimes

Following the example in Part 4 of github.com/cmgrams/wr_data_package_era5, compute the regime index for the seven-regime classification of Grams (2026, in review).

[1]:

import numpy as np
import pandas as pd
import pooch
import xarray as xr

from earthkit.meteo import regimes

Get the regime classification data

Retrieve data from https://zenodo.org/records/18154492.

[2]:

files = pooch.retrieve(
    url="doi:10.5281/zenodo.18154492/wr_data_package_V1.1.zip",
    known_hash="dc942ff2a1b3da6dedd3b0b2fadda017fff9e8fc10228ace31c2209a9be7dc62",
    processor=pooch.Unzip(),
)


def get_file(name):
    for file in files:
        if file.endswith(name):
            return file
    raise FileNotFoundError(name)

Regime pattern normalisation based on day-of-year

Load the file with normalisation weights from the repository.

[3]:

mod_ds = xr.open_dataset(get_file("wr_data/EOFs_WRs.nc"))

Create a lookup table for the normalisation weights and define a function that finds the appropriate weight for a given date.

[4]:

mod = 1.0 / mod_ds["normwgt"].to_series()
mod.index = pd.Index(zip(mod.index.month, mod.index.day, mod.index.hour))


def pattern_normalisation_weight(date):
    date = np.asarray(date)
    shp = date.shape
    date = pd.to_datetime(date.flatten())
    if isinstance(date, pd.Timestamp):
        return mod.loc[(date.month, date.day, date.hour)]
    else:
        idx = list(zip(date.month, date.day, date.hour))
        return mod.loc[idx].values.reshape(shp)

Load the regime patterns

Patterns are stored in a NetCDF file in the repository.

[5]:

pattern_ds = xr.open_dataset(get_file("wr_data/Normed_Z0500-patterns_EOFdomain.nc"))

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.

[6]:

patterns = regimes.ModulatedPatterns(
    labels=pattern_ds.attrs["ClassNames"].split(),
    grid={
        "grid": [0.5, 0.5],
        "area": [90, -80, 30, 40],  # 30-90°N, 80°W-40°E
    },
    base_patterns=pattern_ds["Z0500_mean"].values,
    modulator=pattern_normalisation_weight,
)

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.

Project the test field onto the patterns

Load the test field provided with the dataset. This field contains pre-processed Z500 anomalies ready for projection.

[7]:

ds_field = xr.open_dataset(get_file("wr_data_package_V1.1/example_data/Z0500_20250601_00.nc")).squeeze()

# Reconstruct spatial coordinates
ds_field = ds_field.assign_coords({
    "lat": np.linspace(ds_field.attrs["domymin"], ds_field.attrs["domymax"], ds_field["Z0"].shape[0]),
    "lon": np.linspace(ds_field.attrs["domxmin"], ds_field.attrs["domxmax"], ds_field["Z0"].shape[1]),
})

# Make time coordinate a dimension so it can be used by the pattern generator
ds_field = ds_field.expand_dims("time", axis=0)

# Extract values for the EUR-ATL region
ds_field = ds_field.sel({"lat": slice(30, 90), "lon": slice(-80, 40)})
ds_field

[7]:

<xarray.Dataset> Size: 118kB
Dimensions:  (time: 1, lat: 121, lon: 241)
Coordinates:
  * time     (time) datetime64[ns] 8B 2025-06-01
  * lat      (lat) float32 484B 30.0 30.5 31.0 31.5 32.0 ... 88.5 89.0 89.5 90.0
  * lon      (lon) float32 964B -80.0 -79.5 -79.0 -78.5 ... 38.5 39.0 39.5 40.0
    lev      int32 4B 500
Data variables:
    Z0       (time, lat, lon) float32 117kB -7.056 -6.082 ... -76.57 -76.57
Attributes:
    domxmin:  -180.0
    domxmax:  179.5
    domymin:  -90.0
    domymax:  90.0
    date:     20250601_00
    lev:      500
    lpass:    [240 120  24]
    history:  Wed Sep 10 10:11:44 2025: ncks -O -v Z0 Z_N161_20250601_00.nc Z...
    source:   Centre source file: /run/media/cgrams/LSDP_2025_1/imk-tnp-chgr/...
    NCO:      netCDF Operators version 5.1.9 (Homepage = http://nco.sf.net, C...

Create area-based weights for the grid points to use in the projection.

[8]:

weights = np.cos(np.deg2rad(ds_field.coords["lat"]))

Project onto the regime patterns, supply the valid date of the field for the modulator function to select the normalisation weight.

[9]:

# Map the time coordinate of ds_field to the date argument of pattern_normalisation_weight
projections = regimes.project(ds_field["Z0"], patterns, weights, date="time")

Compute the regime index

Standardise the projections to obtain the regime index. Reference values are given in the notebook:

weather regime indices for  20250601_00
0.8218320408050166 AT
1.1691867614026132 ZO
1.0489664646400954 ScTr
-0.6948883613466961 AR
-0.5738999234610754 EuBL
-0.9323144170607468 ScBL
-0.6369122243737739 GL

Read the standardisation parameters (mean and standard deviation) from the text file in the repo.

[ ]:

def read_wri_std_parameters(f):
    f.readline()
    name = f.readline().strip().split()
    mean = f.readline().strip().split()[1:]
    std = f.readline().strip().split()[1:]
    return (
        xr.DataArray([float(v) for v in mean], coords={"pattern": name}),
        xr.DataArray([float(v) for v in std], coords={"pattern": name}),
    )


with open(get_file("wr_data/WRI_std_params.txt"), "r") as f:
    norm_mean, norm_std = read_wri_std_parameters(f)

Compute the regime index.

[11]:

regimes.regime_index(projections, norm_mean, norm_std)

[11]:

<xarray.DataArray 'IWR' (pattern: 7, time: 1)> Size: 56B
array([[ 0.8218322 ],
       [ 1.16918671],
       [ 1.04896637],
       [-0.6948882 ],
       [-0.57389995],
       [-0.93231443],
       [-0.63691209]])
Coordinates:
  * pattern  (pattern) <U4 112B 'AT' 'ZO' 'ScTr' 'AR' 'EuBL' 'ScBL' 'GL'
  * time     (time) datetime64[ns] 8B 2025-06-01
    lev      int32 4B 500