Ensembles: control members, dew point from relative humidity, case summaries, CRPS tests

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

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_<t>: 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

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.
 

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)

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 = {

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]

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):