From f81968d42a703f0e42eb75acd526fadd3ba567e9 Mon Sep 17 00:00:00 2001 From: Shawn Corvec <26825570+scorvec@users.noreply.github.com> Date: Wed, 28 Oct 2026 07:10:32 -0400 Subject: [PATCH] Ensembles: control members, dew point from relative humidity, case summaries, CRPS tests - open: a scalar member coordinate on ensemble data is the field's member, so an ensemble's control, which cfgrib opens apart from the perturbed members, is no longer dropped from the ensemble scores - variables: rh2m, and td2m derived from t2m and rh2m (Bolton 1980) for models that publish no dew point - ensemble.summarise: one row per case with the members' count, mean, spread, min, max, fair CRPS (and with observation error), the observation's rank and event probabilities, for databases that keep every case - significance.ensemble_tests: the paired, autocorrelation-inflated test on per-run mean CRPS Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_0148RnpHupoNFzgsLvayCDr9 --- README.md | 2 +- src/mverify/ensemble.py | 81 ++++++++++++++++++++++++++++ src/mverify/io/open.py | 15 +++++- src/mverify/significance.py | 36 ++++++++----- src/mverify/variables.py | 14 ++++- tests/test_ensemble.py | 104 ++++++++++++++++++++++++++++++++++++ tests/test_open.py | 1 + 7 files changed, 237 insertions(+), 16 deletions(-) create mode 100644 src/mverify/ensemble.py create mode 100644 tests/test_ensemble.py diff --git a/README.md b/README.md index e67601c..69e5211 100644 --- a/README.md +++ b/README.md @@ -23,7 +23,7 @@ Source: [scorvec/nwp_verification](https://github.com/scorvec/nwp_verification). |---|---|---| | surface stations (GHCNh) | 2 m temperature and dew point, surface (station-level) pressure, sea-level pressure, 10 m wind (u, v, speed, vector), 24 h precipitation at gauges | bias, MAE, RMSE, error SD, correlation; for wind speed the distribution too: QQ quantiles, a hybrid Weibull fit and exceedance frequencies, against the observed (`distribution.py`) | | radiosondes | height, temperature and wind at 1000-100 hPa, specific humidity at 1000-300 hPa | same, plus vector wind RMSE | -| ensembles | any of the above | fair CRPS, ensemble-mean RMSE, spread, spread/skill | +| ensembles | any of the above (a dew point derived from 2 m relative humidity where a model publishes only that) | fair CRPS, ensemble-mean RMSE, spread, spread/skill; paired CRPS tests between ensembles (`significance.ensemble_tests`); `ensemble.summarise` keeps one row per case (mean, spread, CRPS, the observation's rank, event probabilities) so a database of every case need not keep every member | | precipitation analyses (Stage IV over CONUS, GPCP globally) | 24 h precipitation | ETS, frequency bias, POD, FAR, FSS at 50-400 km, bias, RMSE; per-cell maps of bias and RMSE, the bias tested cell by cell (`precip/maps.py`) | Every score also has an observation-error-aware version (`rmse_oc`, `crps_oe`, `ssr_oe`). Radiosonde diff --git a/src/mverify/ensemble.py b/src/mverify/ensemble.py new file mode 100644 index 0000000..9c856c6 --- /dev/null +++ b/src/mverify/ensemble.py @@ -0,0 +1,81 @@ +"""Ensemble pairs summarised to one row per case, for a database that keeps every case of every run. + +A case is one observation against one run of an ensemble: model, init, lead, valid time, station (or cell) and +variable. ``summarise`` turns the members' pairs (mverify.match, one row per member) into one row per case with what +the ensemble scores need, so a window of runs is scored from these rows without the members: + + * n_members, mean, spread (the members' sample standard deviation), min and max + * crps, the fair CRPS (``mverify.scores.crps_ensemble``), and crps_oe, the same with each member perturbed by + N(0, sigma_o) noise (``mverify.scores.ensemble``; NaN where sigma_o is unknown) + * below and ties: members below the observation and equal to it, the observation's rank for rank histograms + (a tie, as with no rain, is split at random when the histogram is drawn) + * p_ge_: the share of members at or above each threshold of the variable (``THRESHOLDS``), for Brier scores and + reliability against the observed exceedance + +Every other column of the pairs that is the same for all members of a case (the station's position and height, the +terrain difference, sigma_o ...) is kept as it is. A case missing some members keeps the members it has +(n_members says how many); the scores of mverify.scores.ensemble over the cases with every member are these rows' +crps and spread exactly. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +CASE = ["model", "init", "lead_h", "valid_time", "station_id", "variable"] +CELL = ["model", "init", "lead_h", "valid_time", "lat", "lon", "variable"] # a case on an analysis grid +# Event thresholds in the variable's canonical unit: 24 h precipitation (mm) at gauges and on the analysis grid, 10 m +# wind speed (m s-1, about 20, 30 and 40 kt). +THRESHOLDS = {"pr24": (1, 5, 10, 25), "tp": (1, 5, 10, 25), "ws10": (10, 15, 20)} + + +def _crps(x: np.ndarray, m: np.ndarray, obs: np.ndarray) -> np.ndarray: + """Fair CRPS per row of members ``x`` (NaN-padded on the right after sorting), ``m`` of them in each row.""" + term1 = np.nansum(np.abs(x - obs[:, None]), axis=1) / m + i = np.arange(x.shape[1])[None, :] + w = np.where(i < m[:, None], 2 * i - m[:, None] + 1, 0) + pair_sum = 2 * np.nansum(x * w, axis=1) + with np.errstate(invalid="ignore", divide="ignore"): + return np.where(m > 1, term1 - pair_sum / (2 * m * (m - 1)), term1) + + +def summarise( + pairs: pd.DataFrame, keys: list[str] = CASE, seed: int = 0, thresholds: dict[str, tuple] = THRESHOLDS +) -> pd.DataFrame: + """One row per case of the members' ``pairs`` (columns ``keys`` - CASE for stations, CELL for grid cells -, + member, fcst, obs, sigma_o where known, and any others); ``seed`` fixes the observation-error noise of crps_oe.""" + ens = pairs.dropna(subset=["member"]) + if ens.empty: + return pd.DataFrame() + ens = ens.sort_values([*keys, "member"], kind="stable") + gid = ens.groupby(keys, sort=False, dropna=False).ngroup().to_numpy() + first = np.r_[0, np.flatnonzero(np.diff(gid)) + 1] + m = np.diff(np.r_[first, len(gid)]) + pos = np.arange(len(gid)) - np.repeat(first, m) + x = np.full((len(first), int(m.max())), np.nan) + x[gid, pos] = ens["fcst"].to_numpy(float) + + meta = [c for c in ens.columns if c not in ("member", "fcst")] + out = ens.iloc[first][meta].reset_index(drop=True) + obs = out["obs"].to_numpy(float) + sigma = out["sigma_o"].to_numpy(float) if "sigma_o" in out else np.full(len(out), np.nan) + known = np.isfinite(sigma) + noisy = x + np.random.default_rng(seed).standard_normal(x.shape) * np.where(known, sigma, 0.0)[:, None] + srt = np.sort(x, axis=1) # NaN last + with np.errstate(invalid="ignore"): + out = out.assign( + n_members=m.astype("int16"), + mean=np.nanmean(x, axis=1), + spread=np.where(m > 1, np.nanstd(x, axis=1, ddof=1), np.nan), + min=np.nanmin(x, axis=1), + max=np.nanmax(x, axis=1), + crps=_crps(srt, m, obs), + crps_oe=np.where(known, _crps(np.sort(noisy, axis=1), m, obs), np.nan), + below=(x < obs[:, None]).sum(axis=1).astype("int16"), + ties=(x == obs[:, None]).sum(axis=1).astype("int16"), + ) + for t in sorted({t for ts in thresholds.values() for t in ts}): + applies = out["variable"].isin([v for v, ts in thresholds.items() if t in ts]).to_numpy() + out[f"p_ge_{t:g}"] = np.where(applies, (x >= t).sum(axis=1) / m, np.nan) + return out diff --git a/src/mverify/io/open.py b/src/mverify/io/open.py index 58afb62..7cbb23a 100644 --- a/src/mverify/io/open.py +++ b/src/mverify/io/open.py @@ -454,8 +454,8 @@ def _iter_dataset(ds: xr.Dataset, path: Path, init_hint, fmt: str, static: dict) if not stacked: init, lead = _field_time(sub, mode, idim, ldim, init_c, lead_c) if init is not None: - fields.append(Field(vname, init, lead, None, grid, sub.transpose(ydim, xdim), unit, - source=f"{path.name}:{fmt}", accum=_accum(vname, sub, lead))) + fields.append(Field(vname, init, lead, _scalar_member(sub, mdim), grid, sub.transpose(ydim, xdim), + unit, source=f"{path.name}:{fmt}", accum=_accum(vname, sub, lead))) continue init, lead = _field_time(sub.isel({mdim: 0}), mode, idim, ldim, init_c, lead_c) if init is None: @@ -469,6 +469,17 @@ def _iter_dataset(ds: xr.Dataset, path: Path, init_hint, fmt: str, static: dict) return fields +def _scalar_member(da: xr.DataArray, mdim: str | None) -> int | None: + """The member of a field whose member coordinate is a scalar: cfgrib opens an ensemble's control ("cf") apart + from its perturbed members ("pf"), so the control's ``number`` is no dimension. A deterministic GRIB field carries + ``number`` 0 too, and stays memberless (None).""" + if mdim is None or mdim not in da.coords or da[mdim].ndim: + return None + if "GRIB_dataType" in da.attrs and da.attrs["GRIB_dataType"] not in ("cf", "pf"): + return None + return int(da[mdim].values) + + def _accum(vname: str, da: xr.DataArray, lead: pd.Timedelta) -> tuple[pd.Timedelta, pd.Timedelta] | None: """Lead-time window of an accumulated field: GRIB startStep when known, else since the init. diff --git a/src/mverify/significance.py b/src/mverify/significance.py index c9cd19a..9d11029 100644 --- a/src/mverify/significance.py +++ b/src/mverify/significance.py @@ -79,22 +79,34 @@ def fdr(p: np.ndarray, q: float = 0.10) -> np.ndarray: return out -def tests(pairs: pd.DataFrame) -> pd.DataFrame: - """p values for every ordered pair of models: columns ``BY`` + ``model``, ``reference``, ``runs`` (paired - runs), ``runs_eff`` (the effective runs the test counts, ``effective_runs``) and ``p``. - - ``pairs`` carries a ``region`` column (``scores.assign_regions``) and should be the - common sample, so that both models of a pair are judged on the same observations. - """ +def _paired(errors: pd.DataFrame, col: str) -> pd.DataFrame: + """The paired test of every ordered pair of models on the per-run mean of ``col``.""" cols = [*BY, "model", "reference", "runs", "runs_eff", "p"] - if pairs.empty or pairs["model"].nunique() < 2: + if errors.empty or errors["model"].nunique() < 2: return pd.DataFrame(columns=cols) - mse = (squared_errors(pairs).groupby(["model", *BY, "init"], dropna=False)["se"].mean() - .unstack("model").sort_index(level="init")) + per_run = (errors.groupby(["model", *BY, "init"], dropna=False)[col].mean() + .unstack("model").sort_index(level="init")) rows = [] - for a, b in permutations(mse.columns, 2): - d = (mse[a] - mse[b]).dropna() + for a, b in permutations(per_run.columns, 2): + d = (per_run[a] - per_run[b]).dropna() for key, g in d.groupby(level=BY, sort=True): d_run = g.to_numpy() rows.append((*key, a, b, len(g), effective_runs(d_run), paired_test(d_run))) return pd.DataFrame(rows, columns=cols) + + +def tests(pairs: pd.DataFrame) -> pd.DataFrame: + """p values for every ordered pair of models: columns ``BY`` + ``model``, ``reference``, ``runs`` (paired + runs), ``runs_eff`` (the effective runs the test counts, ``effective_runs``) and ``p``. + + ``pairs`` carries a ``region`` column (``scores.assign_regions``) and should be the + common sample, so that both models of a pair are judged on the same observations. + """ + return _paired(squared_errors(pairs), "se") + + +def ensemble_tests(cases: pd.DataFrame, stat: str = "crps") -> pd.DataFrame: + """As ``tests``, for ensembles: every ordered pair of models on the per-run mean ``stat`` (``crps`` or + ``crps_oe``) of their cases (``mverify.ensemble.summarise``), which carry ``region`` and should be the common + sample of the ensembles compared.""" + return _paired(cases, stat) diff --git a/src/mverify/variables.py b/src/mverify/variables.py index a16a0d9..4ead5e1 100644 --- a/src/mverify/variables.py +++ b/src/mverify/variables.py @@ -70,6 +70,11 @@ def _upper(base: str, level: int, label: str, unit: str, grib: tuple[str, ...], standard_names=("dew_point_temperature",), aliases=("d2m", "2d", "td2m", "TD2", "DPT_2maboveground", "dewpoint_temperature_2m"), lo=-90, hi=40, sigma_o=1.6), + # Some models publish 2 m relative humidity instead of the dew point (GEFS, GEPS); the dew point is derived from it + # and the 2 m temperature (DERIVED). Never observed itself. + Var("rh2m", "2 m relative humidity", "%", "surface", grib=("2r", "r2"), grib_level="heightAboveGround", + standard_names=("relative_humidity",), aliases=("r2", "2r", "rh2m", "RH2", "RH_2maboveground", "hurs"), + lo=0, hi=105), Var("msl", "mean sea-level pressure", "hPa", "surface", grib=("prmsl", "msl", "mslet"), grib_level="meanSea", standard_names=("air_pressure_at_mean_sea_level", "air_pressure_at_sea_level"), aliases=("msl", "prmsl", "PRMSL_meansealevel", "MSLP", "psl", "slp", "mslp"), lo=870, hi=1090, @@ -137,8 +142,15 @@ def _upper(base: str, level: int, label: str, unit: str, grib: tuple[str, ...], CATALOG: dict[str, Var] = {v.name: v for v in [*SURFACE, *UPPER, PRECIP, GAUGE]} +def dew_point(t, rh): + """Dew point (degC) from temperature (degC) and relative humidity (%), over water (Bolton 1980), the inverse of + the vapour pressure in ``specific_humidity``.""" + g = np.log(np.clip(rh, 0.1, 100.0) / 100.0) + 17.67 * t / (t + 243.5) + return 243.5 * g / (17.67 - g) + + # Quantities a model may only provide through components: name -> (components, function). -DERIVED = {"ws10": (("u10", "v10"), lambda u, v: (u**2 + v**2) ** 0.5)} +DERIVED = {"ws10": (("u10", "v10"), lambda u, v: (u**2 + v**2) ** 0.5), "td2m": (("t2m", "rh2m"), dew_point)} # Base names of pressure-level fields as they appear in netCDF / zarr files. UPPER_ALIASES = { diff --git a/tests/test_ensemble.py b/tests/test_ensemble.py new file mode 100644 index 0000000..2e26fc4 --- /dev/null +++ b/tests/test_ensemble.py @@ -0,0 +1,104 @@ +import numpy as np +import pandas as pd +import pytest + +from mverify import scores, significance +from mverify import variables as V +from mverify.ensemble import CELL, summarise + +INIT = pd.Timestamp("2026-10-01") + + +def _pairs(rng, cases=40, members=11, var="t2m", model="ENS"): + rows = [] + for c in range(cases): + obs = rng.normal(10, 3) + for k in range(members): + rows.append({"model": model, "init": INIT, "lead_h": 24.0, "valid_time": INIT + pd.Timedelta("24h"), + "station_id": f"S{c:03d}", "variable": var, "member": k, "fcst": obs + rng.normal(0.5, 2), + "obs": obs, "sigma_o": 1.4, "lat": 50.0, "lon": c * 1.0}) + return pd.DataFrame(rows) + + +def test_summary_rows_give_the_ensemble_scores(): + p = _pairs(np.random.default_rng(1)) + s = summarise(p) + ref = scores.ensemble(p, ["model"]).iloc[0] + assert len(s) == 40 and (s["n_members"] == 11).all() + assert s["crps"].mean() == pytest.approx(ref["crps"]) + assert np.sqrt((s["spread"] ** 2).mean()) == pytest.approx(ref["spread"]) + assert np.sqrt(((s["mean"] - s["obs"]) ** 2).mean()) == pytest.approx(ref["rmse_mean"]) + assert s["crps_oe"].notna().all() and set(s.columns) >= {"lat", "lon", "sigma_o", "below", "ties", "min", "max"} + + +def test_ranks_thresholds_and_a_missing_member(): + base = {"model": "ENS", "init": INIT, "lead_h": 24.0, "valid_time": INIT + pd.Timedelta("24h"), "sigma_o": np.nan} + rain = [dict(base, station_id="A", variable="pr24", member=k, fcst=f, obs=0.0) + for k, f in enumerate([0.0, 0.0, 2.0, 6.0, 12.0])] + temp = [dict(base, station_id="B", variable="t2m", member=k, fcst=f, obs=1.5) + for k, f in enumerate([1.0, 2.0, 3.0, 4.0])] # member 4 missing + s = summarise(pd.DataFrame(rain + temp)).set_index("station_id") + a, b = s.loc["A"], s.loc["B"] + assert (a["below"], a["ties"], a["n_members"]) == (0, 2, 5) + assert (a["p_ge_1"], a["p_ge_5"], a["p_ge_10"], a["p_ge_25"]) == (0.6, 0.4, 0.2, 0.0) + assert np.isnan(a["p_ge_15"]) and np.isnan(b["p_ge_1"]) + assert b["n_members"] == 4 and b["below"] == 1 + assert b["crps"] == pytest.approx(scores.crps_ensemble(np.array([[1.0, 2.0, 3.0, 4.0]]), np.array([1.5]))[0]) + assert np.isnan(b["crps_oe"]) # no observation error known + + +def test_grid_cells_are_cases_by_position(): + p = _pairs(np.random.default_rng(2), cases=3, members=5, var="tp").drop(columns="station_id") + s = summarise(p, keys=CELL) + assert len(s) == 3 and list(s["lon"]) == [0.0, 1.0, 2.0] + + +def test_dew_point_inverts_the_vapour_pressure(): + t, rh = np.array([-10.0, 5.0, 25.0]), np.array([40.0, 75.0, 100.0]) + td = V.dew_point(t, rh) + e = lambda x: 6.112 * np.exp(17.67 * x / (x + 243.5)) # specific_humidity's vapour pressure + assert td[2] == pytest.approx(25.0) and e(td) == pytest.approx(rh / 100 * e(t)) + assert V.DERIVED["td2m"][0] == ("t2m", "rh2m") + + +def test_ensemble_tests_find_the_better_ensemble(): + rng = np.random.default_rng(3) + rows = [] + for d in range(20): + for model, crps in (("A", 1.0), ("B", 1.5)): + for c in range(30): + rows.append({"model": model, "network": "GHCNH", "variable": "t2m", "region": "Global", + "lead_h": 120.0, "init": INIT + pd.Timedelta(days=d), "crps": crps + rng.normal(0, 0.3)}) + t = significance.ensemble_tests(pd.DataFrame(rows)).set_index(["model", "reference"]) + assert t.loc[("A", "B"), "p"] < 0.01 and t.loc[("A", "B"), "runs"] == 20 + assert t.loc[("B", "A"), "p"] == pytest.approx(t.loc[("A", "B"), "p"]) + + +def test_grib_control_and_perturbed_members(tmp_path): + eccodes = pytest.importorskip("eccodes") + pytest.importorskip("cfgrib") + from mverify.io.open import open_forecast + + p = tmp_path / "ens.grib2" + with open(p, "wb") as fh: + for number in (0, 1, 2): # the control (cf) and two perturbed members (pf) + h = eccodes.codes_grib_new_from_samples("regular_ll_sfc_grib2") + eccodes.codes_set_key_vals(h, {"Ni": 36, "Nj": 19, "latitudeOfFirstGridPointInDegrees": 90, + "longitudeOfFirstGridPointInDegrees": 0, + "latitudeOfLastGridPointInDegrees": -90, + "longitudeOfLastGridPointInDegrees": 350, + "iDirectionIncrementInDegrees": 10, "jDirectionIncrementInDegrees": 10, + "dataDate": int(INIT.strftime("%Y%m%d")), "dataTime": 0, + "productDefinitionTemplateNumber": 1, + "typeOfProcessedData": 3 if number == 0 else 4, + "perturbationNumber": number, "numberOfForecastsInEnsemble": 3, + "stepUnits": 1, "endStep": 24}) + eccodes.codes_set(h, "shortName", "2t") + eccodes.codes_set_key_vals(h, {"typeOfFirstFixedSurface": 103, "scaledValueOfFirstFixedSurface": 2, + "scaleFactorOfFirstFixedSurface": 0}) + eccodes.codes_set_values(h, np.full(36 * 19, 280.0 + number)) + eccodes.codes_write(h, fh) + eccodes.codes_release(h) + fields = sorted(open_forecast(p).fields, key=lambda f: f.member) + assert [f.member for f in fields] == [0, 1, 2] + assert [round(float(f.values().mean()), 2) for f in fields] == [6.85, 7.85, 8.85] diff --git a/tests/test_open.py b/tests/test_open.py index d34a197..08c593e 100644 --- a/tests/test_open.py +++ b/tests/test_open.py @@ -101,6 +101,7 @@ def test_grib_round_trip(tmp_path, init): assert fc["t2m"].values().mean() == pytest.approx(16.85) assert fc["msl"].values().mean() == pytest.approx(1013.0) assert fc["t2m"].lead_h == 6 and fc["t2m"].init == init + assert fc["t2m"].member is None # a deterministic message's number 0 is no member def test_unrecognised_bytes_are_refused(tmp_path):