From 8bbed593817572d20049c981d1618c4d102f2015 Mon Sep 17 00:00:00 2001 From: Francesco Zanetta Date: Thu, 5 Nov 2026 16:35:35 +0100 Subject: [PATCH 1/7] add verify_ensemble with fair crps, spread, rank and probability counts --- src/verification/__init__.py | 399 ++++++++++++++-- tests/unit/data/verify_golden_gridded.nc | Bin 0 -> 25619 bytes tests/unit/data/verify_golden_stations.nc | Bin 0 -> 75790 bytes .../data/verify_golden_stations_tolerant.nc | Bin 0 -> 57675 bytes tests/unit/test_verify_ensemble.py | 428 ++++++++++++++++++ tests/unit/test_verify_regression.py | 162 +++++++ 6 files changed, 961 insertions(+), 28 deletions(-) create mode 100644 tests/unit/data/verify_golden_gridded.nc create mode 100644 tests/unit/data/verify_golden_stations.nc create mode 100644 tests/unit/data/verify_golden_stations_tolerant.nc create mode 100644 tests/unit/test_verify_ensemble.py create mode 100644 tests/unit/test_verify_regression.py diff --git a/src/verification/__init__.py b/src/verification/__init__.py index 61bd24a3..cd781630 100644 --- a/src/verification/__init__.py +++ b/src/verification/__init__.py @@ -382,7 +382,7 @@ def _obs_only_statistics( continue LOG.info("Parameter %s only in obs: computing obs statistics only", param) obs_stats = _compute_statistics( - obs[param].where(masks).where(station_masks), + _apply_masks(obs[param], masks, station_masks), prefix=param + ".", source=obs_label, dim=dim, @@ -433,6 +433,54 @@ def _create_station_group_masks( ) +def _region_and_station_masks( + obs: xr.Dataset, regions: list[dict], holdout_stations: list[str] | None +) -> tuple[xr.DataArray, xr.DataArray]: + """Region masks (region, *spatial) and station group masks (station_group[, values]) + on the points of *obs*, shared by verify() and verify_ensemble().""" + masks = ShapefileSpatialAggregationMasks(regions=regions).get_masks( + lon=obs["longitude"], lat=obs["latitude"] + ) + station_masks = _create_station_group_masks( + obs["values"] if "values" in obs.dims else None, holdout_stations + ) + if holdout_stations: + LOG.info( + "Station group masks created: %d holdout, %d holdin stations", + int(station_masks.sel(station_group="holdout").sum()), + int(station_masks.sel(station_group="holdin").sum()), + ) + return masks, station_masks + + +def _apply_masks( + data: xr.DataArray, masks: xr.DataArray, station_masks: xr.DataArray +) -> xr.DataArray: + """NaN outside each region and station group; adds the "region" and + "station_group" dims (the latter always present: ["all"], or ["all", + "holdout", "holdin"] with holdout stations).""" + return data.where(masks).where(station_masks) + + +def _too_many_missing( + fcst_missing: xr.DataArray, + obs_param: xr.DataArray, + dim: list[str], + max_missing_fraction: float, +) -> xr.DataArray: + """True where the scores must be NaN because the forecast is missing at too many points. + + *fcst_missing* flags the points without a forecast, *obs_param* is the truth after + _apply_masks(). Among the obs-valid in-region points, the fraction where the forecast + is missing is compared to *max_missing_fraction*. Normalising by obs availability + avoids penalising parameters with fewer stations. + """ + missing_fraction = fcst_missing.where(obs_param.notnull()).mean( + dim=dim, skipna=True + ) + return missing_fraction > max_missing_fraction + + def _bbox_polygon(lon_min, lon_max, lat_min, lat_max) -> Polygon: return Polygon( [ @@ -532,22 +580,10 @@ def verify( ) fcst_aligned, obs_aligned = xr.align(fcst, obs, join="inner", copy=False) - region_polygons = ShapefileSpatialAggregationMasks(regions=regions) - masks = region_polygons.get_masks( - lon=obs_aligned["longitude"], lat=obs_aligned["latitude"] + masks, station_masks = _region_and_station_masks( + obs_aligned, regions, holdout_stations ) - station_masks = _create_station_group_masks( - obs_aligned["values"] if "values" in obs_aligned.dims else None, - holdout_stations, - ) - if holdout_stations: - LOG.info( - "Station group masks created: %d holdout, %d holdin stations", - int(station_masks.sel(station_group="holdout").sum()), - int(station_masks.sel(station_group="holdin").sum()), - ) - scores = [] for param in fcst_aligned.data_vars: if param not in obs_aligned.data_vars: @@ -565,20 +601,11 @@ def verify( thresholds, ) - # Apply all region masks at once via broadcast — adds "region" as leading dim - fcst_param = fcst_aligned[param].where(masks) - obs_param = obs_aligned[param].where(masks) - - # Apply station group masks — adds "station_group" dim, always present - # (["all"], or ["all", "holdout", "holdin"] when holdout_stations is set) - fcst_param = fcst_param.where(station_masks) - obs_param = obs_param.where(station_masks) - # Missing fraction: among obs-valid in-region points, fraction where fcst is missing. - # Normalising by obs availability avoids penalising parameters with fewer stations. - missing_fraction = ( - fcst_param.isnull().where(obs_param.notnull()).mean(dim=dim, skipna=True) + fcst_param = _apply_masks(fcst_aligned[param], masks, station_masks) + obs_param = _apply_masks(obs_aligned[param], masks, station_masks) + too_many_missing = _too_many_missing( + fcst_param.isnull(), obs_param, dim, max_missing_fraction ) - too_many_missing = missing_fraction > max_missing_fraction score = _compute_scores( fcst_param, @@ -625,3 +652,319 @@ def verify( LOG.info("Computed metrics in %.2f seconds", time.time() - start) LOG.info("Metrics dataset: \n%s", out) return out + + +# --------------------------------------------------------------------------- +# Ensemble verification +# --------------------------------------------------------------------------- + + +def _np_grouped_bincount( + index: np.ndarray, + weight: np.ndarray, + groups: np.ndarray, + nbins: int, + n_core: int, +) -> np.ndarray: + """Weighted counts of *index* per bin, summed over the points of each group. + + *index* (int) and *weight* have shape (..., *spatial) with ``n_core`` spatial + axes; *groups* is a boolean (G, P) array over the P flattened spatial points. + Returns (..., G, nbins). Indices outside [0, nbins) must carry weight 0. + """ + lead = index.shape[: index.ndim - n_core] + n_lead = int(np.prod(lead, dtype=int)) + n_points = int(np.prod(index.shape[index.ndim - n_core :], dtype=int)) + weight = np.broadcast_to(weight, index.shape).reshape(n_lead, n_points) + index = np.clip(index.reshape(n_lead, n_points), 0, nbins - 1) + # One bincount per group over all leading entries at once: offset each row + flat = index + nbins * np.arange(n_lead)[:, None] + out = np.empty((n_lead, groups.shape[0], nbins)) + for g, sel in enumerate(groups): + out[:, g, :] = np.bincount( + flat[:, sel].ravel(), + weights=weight[:, sel].ravel(), + minlength=n_lead * nbins, + ).reshape(n_lead, nbins) + return out.reshape(*lead, groups.shape[0], nbins) + + +def _np_rank_counts( + min_rank: np.ndarray, + max_rank: np.ndarray, + weight: np.ndarray, + groups: np.ndarray, + n_ranks: int, + n_core: int, +) -> np.ndarray: + """Rank histogram counts with ties split evenly between ranks min..max. + + *weight* is 1 / (max - min + 1) at valid points, 0 elsewhere. Each point adds its + weight to every rank in [min, max]: done as +weight at min and -weight at max + 1, + then a cumulative sum over the ranks, so no (rank, points) array is built. + """ + nbins = n_ranks + 1 + up = _np_grouped_bincount(min_rank, weight, groups, nbins, n_core) + down = _np_grouped_bincount(max_rank + 1, weight, groups, nbins, n_core) + return np.cumsum(up - down, axis=-1)[..., :n_ranks] + + +def _group_reduce( + func, + args: list[xr.DataArray], + groups: xr.DataArray, + dim: list[str], + out_dims: dict[str, int], + **kwargs, +) -> xr.DataArray: + """Apply a grouped numpy reduction over the spatial *dim* of *args*. + + *groups* is the boolean (region, station_group, *dim) mask; the result has the + non-spatial dims of *args*, then region, station_group and *out_dims*. + """ + regions = groups["region"] + station_groups = groups["station_group"] + group_np = ( + groups.transpose("region", "station_group", *dim) + .values.astype(bool) + .reshape(regions.size * station_groups.size, -1) + ) + bin_dims = list(out_dims) + bin_shape = tuple(out_dims.values()) + + def _reduce(*arrays): + out = func(*arrays, groups=group_np, n_core=len(dim), **kwargs) + return out.reshape( + *out.shape[:-2], regions.size, station_groups.size, *bin_shape + ) + + result = xr.apply_ufunc( + _reduce, + *args, + input_core_dims=[dim] * len(args), + output_core_dims=[["region", "station_group", *bin_dims]], + dask="parallelized", + output_dtypes=[float], + dask_gufunc_kwargs={ + "output_sizes": { + "region": regions.size, + "station_group": station_groups.size, + **out_dims, + }, + "allow_rechunk": True, + }, + ) + return result.assign_coords( + region=regions.values, station_group=station_groups.values + ) + + +def _spatial_means( + fields: dict[str, xr.DataArray], + valid: xr.DataArray, + groups: xr.DataArray, + dim: list[str], +) -> dict[str, xr.DataArray]: + """Mean of each (member-free) field over the valid points of each group.""" + zero = xr.zeros_like(valid, dtype=int) + count = _group_reduce( + _np_grouped_bincount, + [zero, valid.astype(float)], + groups, + dim, + {"_bin": 1}, + nbins=1, + ).squeeze("_bin", drop=True) + count = count.where(count > 0) + means = {} + for name, field in fields.items(): + total = _group_reduce( + _np_grouped_bincount, + [zero, field.where(valid, 0.0)], + groups, + dim, + {"_bin": 1}, + nbins=1, + ).squeeze("_bin", drop=True) + means[name] = total / count + return means + + +def _threshold_events( + data: xr.DataArray, thresholds: dict[str, list[float]], dim_name: str +) -> xr.DataArray: + """Boolean events per threshold, with the threshold dim and labels of verify()'s + contingency tables (e.g. ``gt_0p1``; see _binary_confusion_matrix).""" + events = [] + for op_txt, values in thresholds.items(): + try: + op_fn = getattr(op, op_txt) + except AttributeError: + raise AttributeError(f"operator {op_txt} is not available") + labels = [f"{op_txt}_{str(v).replace('.', 'p')}" for v in values] + vals_da = xr.DataArray(values, dims=dim_name, coords={dim_name: labels}) + events.append(op_fn(data, vals_da)) + return xr.concat(events, dim=dim_name) + + +def verify_ensemble( + fcst: xr.Dataset, + obs: xr.Dataset, + fcst_label: str, + obs_label: str, + regions: list[dict], + dim: list[str] | None = None, + threshold_dict: dict[str, dict[str, list[float]]] | None = None, + num_workers: int | None = None, + holdout_stations: list[str] | None = None, + max_missing_fraction: float = 0.0, + member_dim: str = "member", +) -> xr.Dataset: + """Ensemble scores of *fcst*, which has a *member_dim* dimension (N >= 2 members). + + Arguments as for verify(); the result is meant to be merged into the Dataset of + ``verify(fcst.mean(member_dim), ...)``, and has the same dims as its continuous + scores (``source`` = *fcst_label*, ``forecast_reference_time``, ``step``, + ``region``, ``station_group``), with, per parameter ``p``: + + - ``p.CRPS``: fair CRPS (Ferro 2014), spatial mean. + - ``p.ENS_VAR``: member variance (ddof=1), spatial mean. + - ``p.MEMBER_MSE``: mean over members of the member MSE. + - ``p.rank_hist``: counts per ``rank`` (0..N, rank of the truth among the members). + Ties are split evenly between the tied ranks, so counts can be fractional. + - ``p.prob_table`` (only with thresholds for p): counts over the threshold dim of + verify()'s contingency tables (``p.threshold``, same labels), ``n_members_yes`` + (0..N members meeting the threshold) and ``obs`` (0/1, the truth meeting it). + + The Dataset attr ``n_members`` is N. Members are reduced first (one member-free + field per score), then the region and station-group masks are applied. + + A point is missing if any member or the truth is NaN there. As in verify(), the + scores of a region, station group and step are NaN when the fraction of obs-valid + points without a (complete) forecast exceeds *max_missing_fraction*. + Parameters in ``CIRCULAR_PARAMS`` and ``CATEGORICAL_PARAMS`` have NaN CRPS, + ENS_VAR, MEMBER_MSE and rank_hist; their prob_table is computed as usual (on the + raw member values, as verify()'s contingency tables). + """ + start = time.time() + if member_dim not in fcst.dims: + raise ValueError(f"Forecast has no '{member_dim}' dimension.") + n_members = fcst.sizes[member_dim] + if n_members < 2: + raise ValueError(f"Ensemble verification needs >= 2 members, got {n_members}.") + + if num_workers is None: + try: + num_workers = len(os.sched_getaffinity(0)) + except AttributeError: + num_workers = max((os.cpu_count() or 6) - 2, 1) + + if dim is None: + dim = ["x", "y"] if "x" in fcst.dims and "y" in fcst.dims else ["values"] + + if not regions: + raise ValueError( + "At least one region must be specified. " + "Provide an ordered list of region specs via the 'regions' argument." + ) + + fcst_aligned, obs_aligned = xr.align(fcst, obs, join="inner", copy=False) + masks, station_masks = _region_and_station_masks( + obs_aligned, regions, holdout_stations + ) + groups = masks & station_masks # (region, station_group, *dim) + if fcst_aligned.chunks: + # The member reductions (e.g. the sort of the CRPS) need all members in one chunk + fcst_aligned = fcst_aligned.chunk({member_dim: -1}) + + rank = xr.DataArray(np.arange(n_members + 1), dims="rank") + rank = rank.assign_coords(rank=rank) + + results = [] + for param in fcst_aligned.data_vars: + if param not in obs_aligned.data_vars: + LOG.warning("Parameter %s not in obs, skipping", param) + continue + thresholds = ( + threshold_dict.get(param, None) + if isinstance(threshold_dict, dict) + else None + ) + LOG.info( + "Verifying ensemble (%d members) parameter %s for %d regions, thresholds: %s", + n_members, + param, + len(masks.region), + thresholds, + ) + prefix = param + "." + members = fcst_aligned[param] + truth = obs_aligned[param] + + # Member-free point fields; a point with any member missing is missing + fcst_missing = members.isnull().any(member_dim) + valid = ~fcst_missing & truth.notnull() + members = members.where(valid) + truth = truth.where(valid) + + too_many_missing = _too_many_missing( + fcst_missing, + _apply_masks(obs_aligned[param], masks, station_masks), + dim, + max_missing_fraction, + ) + + point_fields = { + f"{prefix}CRPS": scores.probability.crps_for_ensemble( + members, truth, member_dim, method="fair", preserve_dims="all" + ), + f"{prefix}ENS_VAR": members.var(member_dim, ddof=1), + f"{prefix}MEMBER_MSE": ((members - truth) ** 2).mean(member_dim), + } + result = xr.Dataset(_spatial_means(point_fields, valid, groups, dim)) + + # Rank of the truth among the members (0-based): ties span min..max + min_rank = (truth > members).sum(member_dim) + max_rank = (truth >= members).sum(member_dim) + tie_weight = valid / (max_rank - min_rank + 1) + result[f"{prefix}rank_hist"] = _group_reduce( + _np_rank_counts, + [min_rank, max_rank, tie_weight], + groups, + dim, + {"rank": n_members + 1}, + n_ranks=n_members + 1, + ).assign_coords(rank=rank) + + if param in CIRCULAR_PARAMS or param in CATEGORICAL_PARAMS: + # No meaningful CRPS/spread for angles or class codes. full_like keeps + # the shape only, so the (lazy) scores above are never computed. + result = xr.full_like(result, np.nan) + + if thresholds is not None: + threshold_dim = f"{prefix}threshold" + n_yes = _threshold_events( + fcst_aligned[param], thresholds, threshold_dim + ).sum(member_dim) + obs_event = _threshold_events(obs_aligned[param], thresholds, threshold_dim) + joint = n_yes * 2 + obs_event.astype(int) # n_members_yes, obs + table = _group_reduce( + _np_grouped_bincount, + [joint, valid.astype(float)], + groups, + dim, + {"n_members_yes": n_members + 1, "obs": 2}, + nbins=2 * (n_members + 1), + ) + result[f"{prefix}prob_table"] = table.transpose( + ..., "region", "station_group", threshold_dim, "n_members_yes", "obs" + ).assign_coords(n_members_yes=np.arange(n_members + 1), obs=[0, 1]) + + result = result.where(~too_many_missing).expand_dims({"source": [fcst_label]}) + results.append(_merge_metrics([result], num_workers=num_workers)) + + out = xr.merge(results, join="outer", compat="no_conflicts") + out.attrs["n_members"] = n_members + LOG.info("Computed ensemble metrics in %.2f seconds", time.time() - start) + LOG.info("Ensemble metrics dataset: \n%s", out) + return out diff --git a/tests/unit/data/verify_golden_gridded.nc b/tests/unit/data/verify_golden_gridded.nc new file mode 100644 index 0000000000000000000000000000000000000000..62fbdf513b85e8c91a9ae142178aeb7f08dd8092 GIT binary patch literal 25619 zcmeHP34D!L+CMkBA-I+Zf?95)HLWFLA4~e)$RZg_R=bKTxd_P($z_6e!Z%b^zp0`2 znQ57}C><@TrpvF|QPodT#n4YtYir+^LG?S&Ip=@#I%2tPkl)vwANRcf=Pb{2p7WgN zJ=>cpanZ3A%hxXN7Z}J8eo`SgLsG(3!hytR;uGScll^R~{6YoIzn)5(-&j60Z6f0o z6_@v87NS(I!>l~=w=h($$F>;C=g%s!K+2ldJt{FbXM`))lj+K1!90mL`tkpxUAcLg zIqrxSP24U|ddA?0R$=YJS~eM)F}QU^^RVWv!rHJPo{il?UU1{`CLw;sz7|V0%Hn6S zSg1-0i2bXP#I-CIR)tgPzi0rLDxd~E`Zpq)2ufp=UWF5b*@QtdVY)2hZFTkAw>pwa z__J!v#{ctHizc;-%E`&i$aFhByzNyu8`2$Pb|J8;TZttSn2oF7!mBAGl^4}DG9^VQ zjbx>kMaOrG>7Ejw&^;|BDza+~vpq_5e-6~NiB!mNy|2ppPD9>=wJbrh~nRfur3kTqR zHAQNp?6VBjXMa}gCH1fUBpgfHl_BWDmO8O1DKTt}@)(xdj2)qL%R$EK`tyh0O+F4> z|I!je!DbNEQyVN|Xn!{uxJhYRMqFz*L9HOWqTD>oCdO}oO zmoj#O7Hc^{WvpVJ$LZmoF>Oe0&ZrSfrO#l0JsO!uxqUX~VcjPTxCfQ{Fp#s_N=c_E zf$>rZxu5Ql)Zla5S$QR*if%MZTKUbaCPI(V($!Yr{+c{?Y-UzgPiNLB7Yni*UT@Z_ z0&DdX?_GMsz&a7jIc*>%l?aj2qs)DHSOl0T(Or(7Hba5%DA913JSX zARlE8^ujR*`8p`1z_`KZ_byovItYB7Sy`eG4QO+)nE`*UsIXqHucCv0Jte5}SgK;_1igXpP%S_BW@>cs9j4L2ceq9eU)kH3J}?Hp*;;<^&DZGQJ4T~}?^ulvzT-4H z_&%f2!FRGo2jAy4I{3mjLZ9IKf|ehAr)qTYou$#iceX|c-#Ho`eCKL(@SUg8!S{8I z4!&<_bnsoI(ZP4IMhD*|8XbJ!RCGj3@u1>m8h%S*_@g17wD1wF8T&_3^yZ8$*UG=G z;deB=Lc=RHyh_9GDvW;jxYN>eM!7w!6&-%fljCt_$t02To>p$KDz`?bkC27_snI=( zzD}c~-#*jm=(qJ69sTyXMn}JG(CFy5Z5kbXw`+9p-J#LJ_um>Fe0OSe@I9c>!S|p> z2j4>)9efXKbnrc*(ZTnuMhD+>8XbJkYjp6vpwYqiqDBYb3hKoU`vPB^MhD-D8a=m` zthka!hdh-vI^?OM(IHQeMrSo;{dF`t2VXgPy!AcY zP0J6y$~@V$06dk(r5Y&JK&b{wHSjPrAa{5M7wbMG&Bf)eR>v-~P19A=u;UYOKX-iS z&EoTpkF9CZTercyGFsZFCZwh%CdWjDd4}e?@`mPQW%#x^E?;+BPDCU~ji;G3j#nXK z-?uTBa|xx!-)zGld3c~2r$vyilhbi#suK#aX7n$nPDw59W3UmcsrQio>d+4(=GogYO-MaSARGP7N7I-Rl)8s)LO zb3BDfw{FX|W?NhT??3OOj6QE5K7Vj(zYQxzbJaS5R6xzeeil?gB&p|{)5~+Zxjb#( znak7m+v1?{{THrtf4u2L`g0-Y7a5$y!z6^wT5*|!av@tzSigR+oONH z^7j!tj47#`J2dOD!{_pSe|Y~?<>e=gmE)RD`01VReJ(HC$=IY5hg!`FAL$sdsl~at z_9Gp7`-g03SZx_kH|CA*Go{@~N6ddFPfBb$(sA{tkS{v79O+o}?ccARta8Yh(xCU^ zL9ttJkz-F#js0Jx>^G`kTK+%DWj;5SgzWut#>qX#+GSTSWk0>!Xz$+f(T6oT|5}F{ z&3^n={^c|B6PbOa<4gbLIme$E=>RUDGal;ttK-_j# zw~sDzj9A!@vNPItF=N`E5lJeYwk^z0cOMSr{q4~2I5cJo^Ay&f^$2pYcasif3d#>9+i z6?VKd$>`DW)1e>Le(tX2jV|v$S*GLrMp9=-*yM{J8_Q?i_?>IoDkJ-7esHhjtM6LA zVML>WQ|g~GdjDrvZb1GiqxIh4-91~LHWJ!)t#_{1>ARMf?Xf##TvNqWw77-BO;r4t zibFNbRH&oqz~1;S^fO9c#J4S{)dP9jDte%b8>$%TkRR!2?^7xt;xL5~gN_(<#D(_2 zL@34gs)v?6k}ap>wi+|O2zq4g*V~N&p9D;4^y40*%f;LUySD5%E}R{5?5(U#M%tqt z&&M_S%7|R~Leh~1-x>k$_B&aSzTXHAa)-}6{h`r*@TukXcYkDDTUu?+v=!@(vi@z# zT$}&7(P_fWsWXp%Vtjfwu>QpH>x}JVTECNc>=UC>LZ|YNmi1QRLhw|NY znC=@Jn(&;neq2w^^vhW+shq(40t%0Ui2(&YR-0ztQ&>G682RK%p{Rpr^K}}@>q-F* zn)B*NibF=;3pLi;uFB-vVtPYa0DF@@nFGxVd(L!zQ1oSw(2%T9o?&8m!H4&ALL6OO zXZN-jrS(+iH^i!$uUB=APf2C8S5D{YA^bm0K9hqZzpOD%dD@rbphNqpNykKN-$3dG@?ace=~oBD_Vbrs2)`KYMt1 zd-2a}BkD4Ol91W>8CU?Xwj6H|@5*#%mdl;t%oUf$utpw{o0H`l;mJ(rH&9$da`-g_ zHyKW4R@TY+^5DJqbA<7Ge!q^;p<{y@(h=}Gg(buh$l<)Jr|f}Ysw1#D<9}N>V7otS zX@1@EK<66lM5A-hNKNr;aV1Jq8^fiFR%n{Ko5P&&h%7Kom6LjJT2{hjMJ;M5G+3}^ zZftT~BBQ*1!kMVtxHr3v94ZgEHz$`roS^84)w!C;d#4{>3~v|p`myskGg6kZu%Px6 zS*O&vXzm4gK3ZF4Jh8EfOs$F2JC++isiVkMhH_Pwxx(~Z>!*f}icR5lu^(lUm$g>} zMYpm4AZ>@*C4w~a8i3LmXv!giv>5N`msX<=@%ikQ+elMZ8=0laC~w7-@vB)P zztSztOic)w8b-10=H#Pz$%tw8gXzmkZHimG^R+bIU_MOMmJ?$GHSE;zAPuK;ykq^B z|1wLwv_2>D;1IujKfPFIm}hydW{J;f_`HTMX!s(>GmhkqJfFv#>Pc=V))F2lKbo$(ZQl#2YmBT zez9;`7UZ|0G4QQ)=^nQgELIHDkij{*u5@RfCoR`C*p>Nf$7*<- zhM(b>Z*H{hztL=2OG5>jIb>AKljO!+V6)70_t%W+&aWBMU0*Y%JHFc!y9*MgN&ZSw z%#DCCAp*t(2>6Gd{l)z3>3**~Yu3QHUpiRb)1$sKSNOp%t?+|U*RGVzL3*)tkly*+ zVKO&PvjIP^;VBw^LBmryW~}C)dQ^}GP}25opKUG#6q`NtgxUG+cev&BY>(O)p>;k5 zN&Q=-M9DHL@(AO%{WtSMmAc9-a>m72)}`>^*tqzc*RM;|(xFy@%#4j)#Or@d9ZLhX z_-6L>!D;Yf|B^XaFY|LN3Z-vUu7slAsIxfTJXmLIc#ekWYIq*UjP;+ou#(yR55KQ( zRQj67EHiz@W5)FLjv3QeJ7#R>U|#Fm?~&xbzEN*f*(*$H;G{ZTN~X!jsa=**ljq%A zl*zx&>1K0yL&J+Syja6aIA-kZi3`kZ^2c?Pk0vf9Z|j6U4Kuf%u3eZheF|pAbPi$0 z^yydppP$X(JGaJfJ}hZN06X=^wZ-Dy-_gIa9kg|tQGcf4^&0+M!y7n$AtJxBm`l7U z_nC_4+6sL`Xf`HwXqdhgG^f)yf@ku+*~bOs!%qh6JOlgc_4RcWredhCmU6bh zogI=dF!uQ{b#py=bNY?Xst8ZDH2V(?%}!;UGGn{U{s4C*apGC^YV;vu?x!TciW3+s zPGD#X*gERe&XVdR8rPSNa;UbIIti`_^HAER;q4mUq2YgPcqhkg7mxcwI0-RCiv5e3 z;#?rkMU)E8V_oi zUXteY!x}!qaf_p|Ripuwv`d?vdU5f1praoO5A99J|0GUI=o<09;Pbj_d;(W|hp)Jh zfEbq_WRnZ;hXCX76&Q!Fz?f75<9z{)sS+@zO2C*Z0mI(|!@&c?p##Hb1H(rH!zlv~ zNSYE@GR^P%SSjOmCF(;ozpaOBd@ZHW)w`kDcfv5`w(W%|nQMg-x%lQGjD6?7W*`j# z$!TH7`^RZP&yW+0>4x=(Ghb|zh7#dt;TIeD?7!@-Uzr4`mSd)mSa?AyqtbPJujPK; zS18BTi)_4Ne-EbAhXnO$VybQwh-=hT7|S|ZG*NSjH!c0Bsl{g)wq;Z^dBl7DIz*S# z{DP{$ff?I0Tv5aHCz{Op=}$A6G5tvAZNmJ}}Kga%dhHng@pFfuVU|XdW1v2ZrW>p?P3v9vGSjhUS5xd5mpn9_i3L zFf$3S@JfbxETliQv49q#M}W~Iz~~WR^awC|1QE zjB8+wYha9PV2o>EjB8+wYha9PV2o>EjB8+wYha9PV2o>EjB8+wYv7KdC-#{=ZLspR zc9r45%{s(-_hz0eEC5jEvsEK745TZsDdy4`JDoR!C4P)#z_7Y2Ck1-m1 zjo3?~vBe%^K*f$N#u%gj%+5TQW%T}vFXp{(zB`GW`<$7bbLQ;XDQ9PQ$JDPE=#;B) zE;|PYg7iY9JY@|U%Ufhj9OqfKN&R}w>>StId1^fS4@6e0V*JN-LGu*nIF}vK>TvTQ zM#`}_5L7Nx#{hPOJ;_5HSkaKi{=q$C<05({L`L)?`DqrDXh;9`iRjrYGB&1;PlcF> zgwEkz>Qt-rL8Yn{x`lUnzfP4(RjO60PF!fQm(x#eziz9k@6^sME?<9g7)Pjjv*#I5YEv$DO&H&aJuOV%`QM3rRA!>!A(Fkd3GDj|mo^MjjmLjDn} zGAn1^s?_UWWz}r7GS;SE2@!FW{8jYeQf{JhIunUTqwys3j+mX+%}wU>9{owvKWZ}8 zfXtihHICRbqAR`FY0dlXZv4!R7N#?r*`HWzY0A9Oa?_*u)%&55opeI}(W{p=AoE5m z`)WVSLo4Hc7M4J*7Sgq6Z11?tUBuj=eb!C;k{i2ImKr-BRgin=DlhtzrYW45FY0XX zva4j>Xy%Pw?$yT5ODki>ZV}=WR4Jib&xl^#Vxz({XKmiEudaV2UlafOA9KqL7uRlU zVHrm{#v0T6_TDj(34B0gLoHJvJgk2&Q?JOF&JiXbZ=Y%vy{pjwP2S!gYX4J%*P6ao z49AAmv`=7Oemk3Y6|Me{{F^k=++OQyGma{(02823xZuF!Enq_!PqCePX+Kz(#}e|p zbncbFFzciAi48dV#D{})PEPR?_4JM4_2d$K$PJ5%(m3pL)+a|jpue3Y%q|&yVrP{l z=*iA$B<+8Mc?k4n{Mg*rRJgh)+$QIB)k? z@N3kJd%vK>k7DJ$ikge8}~c_>kL3;zMqb#E0C*5+8D#NPNf*miUm{ zUgAUUClViWJ4k%U?I`gfH$>t?Zf}VXxqT!)LL|hpCBtG=%FY%$z0ErKM21@+G zqN5lo@uBA^i4Q$TOMK`vM&d)Cu@WEpOqckOJ450_?o5dfxw9lbj9?;zRB_i4VE!B|hYCkob_hQQ|}HCW#NZn*|?a3w{W?Rl?f@#yCRz z?<#S$^V$9~pZ3y(91wi;hiKP9iH~+2lK5!XVTq4+{UGttuG119a?eP7$UQ6ZA@`ib zhurfLA9Alse8{~f@gX-=;zRCri4VCqBtGOmlK7DOyTpgw#}XfM|B(2Q`$XbHZXRut zl|H^eZeEEGxy}+F{g{iyhd%ivKJ>{i@u82a#P1;{0>x#1A>JU9#D|_GBtG;Jkr4B| z0Q!`a%0r(@5+C-dEb$@NTjE1*6^ReIJ`x{t>qvaatt;^%SCROTt4e&x^_BRL%O~A0 z`$Lda9&#H?e8_Dg@gX-@;zO=B`xidSNHb9mtme{jtlwhu$c_>pb_$XB&^J`#L*Foo z4}CjHeCQi3@gX-x;zMq%#E0BCi4VE)5+8C0Nqop1Eb$?Ch{T87p%Nc*6D2<6j+6M1 zJ6_^L?gWVsxf3NmP<2)eyg;$6ROD zGSnpGAX}ol&+3$ioSFZJ3j%0m$=*WTo`~Zzh2Qul20ye>g3~NiNXY0MR1i&Aqh_i? zz5od8HoPA{|MZ@d7W>JXdrzCGT0Cr6zLlmU`HYi)CZNg;lWBh-BwNgd?o3TNI7c4r z+uAU&egi*(`-!ZSPigj-kY7X#hpgcEQcixdO3qsK89P}HTFlMfuv^4zFlEj3>2h9l zHj$sKD!P-!g{4PGGv6<=xfefK%1GV>`Qm-2(Qmf2%=F31l<}B^-8G|E%nWQ+KUfQu z+L;5VGz+(UlZ7d9Rht(H?0A+0rOVir}wF2=U&TUH^X=NqJvk8+sk^riDA#(q6 zn=mq(a*L2lb-szvZWvS~4}L5kO=fG&ba|25S6Jr%9?{M}A1PSl-RQqU1h+T{i*dVC zzZHoa)ZBNnh^i_5zQzy=OVC{uA}ye_81-h6c2Me-eoe9-$-{gPgRsR`5u$aYtUNzm z<~tYMLCiNUcGAM-G*z{_xuX=&e$zkP$J_zA?>I&DjRnU>br#SMUmha#tcnw!sUDd5 zP$Ta?e4CDwk8b#)7nznoGe3N>r3fFs)=uMf{iDNIPAW6=!?#ElmCx*G5MI>(|RTKY2QTPA>z@O0x52jJS-AZ1)G2fd(Tk z=j?ok1J;2CBQ9*GUJ1Ip_*f4bjJSiDC*zH0+=wt3aqZ{UuUXh+n2-{*k{DC1|vHg zEalXiz>>j;PK(#*-=tZyS+$Ud z9;|>?qs~Yn31OY0B8n=!pryZ97nCIxO4BQXrng#XFz=&>J~C89g<3L;boxDe8NX>c-i5fB6ft2_NwwUh+6y{)&~N7AGM@~^^xTgPXEd3F z7{mmO=J_G!iJ*>JKa;6Va&&&8B~Xm5OeSkC)sZ;XV^cGGI#rU#8yo5tMHk3x3u(^_aUiMD3C!W>r?zsD-f^np=NVn$HG zxC~*T{Y|dK>lNHuQPa;neRquF+hp=m4m}6W0s;e?1vK^#2+__#b$Zg-IS3C+2SN8>NYSZIr(p*Gbs?0fA zSArEQcOS+`-$8`&TwjjO3H0@_xrIKT4SG4~Kh(=yQ*%2Zmkr|DPl4I!uE*$Va6Frq z>+?%W*i*tKB+SBZ+Jr^OGf|Ax@{1j%8?`^KkiRro{qH&DZg6S{8|qC%>0BKVG|aO`&l8I z3;MS=olwXpbyi$DN$YKOB~Q~U6qCVo{JMUjklCLOQ+i!6AIJ5h^`YMX5*_-WOa5D@ z&MJql2Ij8P>5OvzNgwz66VEAsIGwBgW7|tgm(yV_4>e6z#($rZ?MAh9{hxQM7?yXw z&R#_t)70?fKY!)7NuT&E;;C$+3Uc9a`>GgEj6Q zuPUTWr930;DQr)qZ_|^%?8(+X`L5zvs^Q2S z?rF-AzdqXTc=N9ErQt;@2fH+(?>)sMd*FzswScL8kN)~WoBTQq`)KsuPPA*W=bhV? zFWyzY_MVV<#s98SaCxKQ-t@S|RWHT+AHAz+Ljt*{1UFgMJD{7T`qIxQ>iwr^#J++7 zkH*HUUHdn@{E*U!TSNOYe5PVVVD#*G)!|z2!z(HM{`e=x#*{)YhZ^e}Z<&~m!Os;? zFX^tXT*BzD`|qLUD)`1;`}XsAwNQondx9p!tF8%M$~#SsR~N5qFk77!uQocBkkD&# zy!xVL-k$SD#;a35>)tAw>aqU&--d4ZEM9e7Tiy9MEkEO4+Qd!M;?@26&L8VQv45ej zZqFVauRaI}^SM1fUcLS4?n80(_$B*BUB5)-rSv;(L?g)tMf3Ivu9+ z2gXmH-*|q!y8Gh1A7+iBa<_WAx1f5iX`nniO3N?pHEr>XG4X2ub{C!Zw9BVn|7%sA zw8TQ{oa;|K+f8y(6AgyHrX0PYOwazX;BR4PlpS9+x%M>rp4tCfKH5EV@z({+_JVy9 zdRdrC?>r8i5&kk4}(0#j$>OXDVcdpmn)R_9N8#Xj?R*yaVxN;uj zRmHVHi6iAs7f=&wPWC%(=c0z)C}tWnuz>o&zuK&W?>nlb==9tjOhgUOG0La-$)`%W zjthA zHTF+e_1;oLiyI%so9XQ1f9B8AE?&KLaoL57Lxg|O^O^zE{;TETPd_RW{IdVlx;DCj zLhTRxtnDzne~pguD(Jzn`^pcexV>9o*z@pswX;+3`whmW;{`u7UT}ui3wq_FHTk_- z#jDjHg-$IOE9y;5$HiNg^l)t&uO7Pou*j^I@hb4I9*Kj$>_q*pk$1bw)GvS+-6>e2 zGOd4bkqe85Q$GYgx$=PBmNw~pVEw#-0YT z6$`M?gtL_C3Qxu6W>*{6oKkZHlQC&Bhb=iIzRg;WvCdxAT?A2Wy(&Fbh)n4nvOVSN z@h1}QAYpd%PG3HRVnTdp-paupfLP72itz2E?wX3`tKsznD<6>aGSko`HgAYKNnsT@ z_Xkes)rudW-8;t>Aksbj?DfoaEywTML3@X$#k4kTP+}nlNYx*&L3AxPTMl`+t967dSlqjn z`lCJ^2R`aE=OaVNwPZJx$MLAwQh7^!90#$e&r*5V1INA8&xwjq5Bfuthh4x&IkX4m zWPQL#dEm0bPBKOPGBt_vGDSON3OO>B>jk!CH{{FIQh7^!Sr1@K<;#fn$`tmNsp-FD zA4~e%_(V%KJJ{^NEMW5kn;+Qxz;+$j`hnN=1J)UTKj6Xjw#DD69S}p+2i!(dz8;T~ z@MsB-k?>fG>64`i*FCH){?4Q5i`WeVn%>+2_;L39#)1@mkM726o5Kl7e{2X0LK@_L zn2V6r36&4Bm?W=lI*XgqG=&r5NcZBt>IvWcG~bT?d-(x^%tb|c=yq_V7rlpoAS8m4 z@PdNPxiZLavh>{GjE`ziO?r5Cnrtxgz&{>0h>i%0!NCS&0nsG>OhKQpp1-u=9D|W< zs$_|0MMuUQ6Sp)BMz;TwC6N`@?`L6egRwBrAx>r^u;&{b!a3VuEF$t9&H3~hhuPZw zLh8y!alsQ^wr@b8mN>fTOSQ9S9ZrWKznI7`Y_9#K=N;O2Md9KvRmNQGE+YDGaDPU-;4^i91{Ve>?$5C9kzJH5x}SyjWyWXrz6|sN z)py&|scw2aL&7s9JWIl}DbDnLm!08GMOYta^^cp($}@KtZ6HC{#D25X{c6pv);Vu+ z^x_)6`LJd6Cn0N}`B_Zuqj&fv5?(6dWfERaF(K8Sulyu^lpa}9W18m*j#+x^G!Io* zuW?=pvlp1@^Ias&URb8jXD=$#WB!7&xpA$%^!Gwn_uy?ZX+_`SBHuQm!C_j3IGi-U z=3xInbm?*e?mkt32HCuY@ni%NOUb;Gsb$Z`7LBkw(?MjTLkN|Cue zYtNdvCSC;%YxL{Bw%^99V-pul7`-fB&HLl6X{{HfV_wITp<74o>v?-)yb26@cF~LK z^S?;v58gWD+^QXPx4$y^j@@Rui&|jND{cM%uCXMYKQv)<3|is0 zUyH`BNXHQec2Av6^J~<(S8Fh(z@S@~#JNW;Oy?)ew*Rp=ZC}1!J}E&f<5l2Q%N|!% zDL=XIqUMF?#;d^A*=6i!8*FH@%^Ro@UNuoZK>3>Yi!S7#gbzvBJobt5eBWKvr|-?+ z4nV(q(t;jh1&@pK3;(@7OdC42<>XCvvb~#?5+SA&7hNs9jEt_~ud?NoVYPuYW$Uc% zp4kUgG%Wzdd6)V&*JfH-x=_4IrRn|unuJp&d|kpfC_Y{6aM(MX1NyCX zV~@=ca(C^eR=dtW82qQpVa0Q;`@Yk6&nwT$uSl!gX0x&({EM6`Iv-YM{yzTtkeHK7 zzb}7T^1RGm<%i`Fs{{;JY2b=3Gf8;>X+C!XeYru>>Chvb=l z^ty5-@?pO0#r7*Lr@i<6_zEY?{5Rlx!cgSbihTm2OcCOkZg_d{uVA8gIp|=#T&XH@qvfVJ3pRjHF zrue?r_b~CW->wGDzj#MET6tl!(^IY}6+7%Ly-U5X^ho$&V6|&^lx9cUWN#Cms%-1* z*Z%XV_mpoQv?X2d-%z$!?KQ8X_X8zhbN-ZVTdpY6mX>HackwkP^|V8)BVG@bw1DQ# zYs~spDdy^SaO2KoWk&YERAY^s%HH??j4SSNNqIJGdy4ne%SyNTlXq9?ep4Bry7+o% zlQYUmugJT_+oURw+BHpD8U3SjaoW`A0e#YxI_<+hXc&~LI4&zys!F>HO2H5(r_M!h zE1%biZaCj3Rhgqs-Qbk(mf~Fd{IqL+SC!5??kVG`K6f2+PHJ)IrZRtec)_xXR}|%1 z&eW&Q$;#oe&%F=)bz6D#Uh&V$)i|g8(kRczaqhn;f$NXod|dF3;yh$j?#`}9ly1|G z&9<*{L)qLsU{*rhkIL|q=VNabzNaK4Opj?_u_J`u%bln zM)tRrKW^oE)_v6@#pl$e64eeTD z=TxO+gAz*j0aul%v(bls9e+)^K5}M6XuxmEvyc*j_m5mv8XqgrG`Qb6<J{b2{PiWOgdI@EJ8xg}e*GKD40T_&{g+ad zt#iA~{rTE8<-6LhgOW(HGQ8d4?iZ$AQHEVSHGRzTYf2xVE%TE-Qk1W2ck`USkRl>vvuWn-Xo`y!TDz>A+PB%1-@S_%|OlUiG(bjdS@o zjX=NfsuKTho;Ix(b8EaxUEg9h&*1Z=p5n#`mNr=%Ph_u`*UyWBB-~iSO(e`;A3rV6 zX(y&JmhX9Ed0R|c)o&4fo3yfYhHze8aR>3OD@>d8=V$5n1aYeh^I+VY0^ZxcY6zc| z*SE5D&g4s(15X3y@>^k26y zAzQKsEHFQegWIo0t{$13(y6-sn)9sAGfX1?-M{Aw8nLs>(xbX#-P8N6PeyQV3MW~3 z3z3i-OP|mM9PVD@s#QbVRA%zu;WKJEXk)VOUTD?Q(WCiA(Amnmnb>%K*4B!X_uuPg z>ap=DE8e|FKk2t);(+K=_0wInggmrnnXs%Kq$^2-ZefSR%lYunf}$RCzQFKzl!Kpv5AEdRP!87*t~2mKQ4h)iLmuh{m5+noCfGrus0ZaT zTYp(mA1LYrg}`-Uemwji>*Mm;4fKe~liOkn-8rIraJM ztDKqMFS~V6@dD|KF8`Xmfm1lqs(aNl0WLJ(j{cLs7eHB?Rv&TG+idp?Ue~*^8T+E8 z-UcBO4wZ12gxMD>Gi|d9UtG$eHJ`N%;Ynr=mIfD^URdA$%5j`;ZoieTW7#)U>HlW$ zW8Y5IWA+MJJys>`OR-D%<6G&@M{GKrXEj}@FXqMqw!}S zq-z`|Yy3@S_rtRo6ic9b4cON?^_YE?Q;*r#IQ4kBgxS|O_4(|pn|jQ@w)wRE?CR;k zbUd1Yz$6~V0RC=5YlC%0x$aNomR!rp8*E{DSHbiwyI=|V^>LpOui1VydblyzI?u@D zk^EF3X#B5;JC54kqAEA{IK$h&Ax^stsdYdx7`$$e7s_Yu4RQS5W8E?Qbewqe3*_Cs zACBEE&AQk-A*OQ|w~2P4%fg+6*xzffPVchQzU0N;%^<_RrKX)DaW^lbQ=ky^s)oZG>Rb3fpOo(%i|{iJ%O_n%GqZjpXcXUCpR zqX{}>%^AdpKGv#4%Zcsl>?BsjX@4f;&W%4gq(%!dILGS-n>gl`{ z&%`z&JC$!IvQxo-zH8~Tm{2>FZyS25$QXRSo#+qcp4*OLyUV>0+lFqLu_b=%nuUwI z*MF`ABow{4zxrPaFeu+n^wiAf+lBrX+lKyD9+%5Du*U9ZW{mt4@A{`VFSw-WPrPw6 zx1UL~e^y<(IPvWJoNNtxgH;sF#G;^+oIpp)@$rJN76X6djs6Z_;~6o ztoI$SaG@C7GQE6SeW>lKlXFj$52M=eT$m^+Y_MJio6g8^Br4|-# zo#vt@)vf0EscRwCyV5s}hHs(YrC4{pp~1yTH8vb@u-lH(D(Tm+TL|4c{CmB@11Al3 zRFm@6`*>S67j@04Kkm%km{a9v#(JWV3m=AiRcJLtb)HkEMeQ8D)oRYWbDbD6OdUN# znfv>yA!>pB&!^0u{IR-W>zS?(OOI49ugY01=eZ&3&u;mY0*{8PPh*}ei}l?nio6QEN0$t(nKZs(K=7NRvO4s;iOv=dLPPs){exTW!F%0gGqsXc?zlp&|q_NY9Xeov^peX-Hz?-AWonsnUs(8-vlB-Fl} z-)&@?vhq@De1m3b%F`*nzgPY#O{u@%YwYJE)9Y_tAUbd&)x+AZFBT8Mgqi);W47?n zV>SZoF-98OT?^~!FX<_w>3>ZhuJ!-Lzt$P?N?pRwg*-#E|KI(}Eg@4!l&qFM{a(8| zV-vG7h;w?%Y96$~=Y0?kyL>yfEmOhjBeq$)uHPsv*EUIbvxK)ucq_$zGvB+XUyUWR zx9HGhqIoQ3dIo-H(v)(9VAaINHI_{5l)0`7CE?#IQE4oVgDsKEz#C9GwR=3 zqP4f2%w(iLh0}15-6L*4YYSX;oXR!*7j?Or9ySk*KXP?m$91iTqpHF9p2)B~U1l#^$&*I6z_UL`5d zH{kppw~$vgV;+5Z z<$n5P8%qXbHU5K$0ku~#i%yZaoAN{Dt=Dd15L&A7rC5~q=taPcE zyqA*ngi8n)MGvCM zLDdM6pOfL-m;YKSZ$4fe0_?1E1~85%f-j1bm#{e_Ykxow)GA}}LFIg4(3i@o;0P+P z#qv_|*Flr04|a4F7!-CyKAxWf#`yq6dC<(3FDvAuKIkVu2Za20UH*S90&dfnd%VpK zvK_`w{ygRLVL|Hs1sD7ctLv!iJ~HMjaG|;C6g_^X3(e)d24=9##kqQE&V2G1WibN=+DDdj1kqC$(OX z@AEzdnyJ?go?X*?b34^DefT&ha@Y{wFdsG^t^HUlJSqrFO`2PmDlV?E-UdC%AxEIqM*l??@N3a1kYq1}7Nyc3rHT1ee@flK zL>)q;uP7`|JOFoT5TEfNm$kxd;?0viAr)zHkt$MXqW_5$L}76v{EUz&4SFJXoCZq* zIvNNmNlF2JO>>J;HAOl|VNrs(6LF){+oA-KCgMY3F@n$%38t_-LHLLa)N;L0)okrp zx>b-T_EHN^f`Ad(sO2I=L{3uZK@b%pX&OWzh~%ch`=SJ)A5xY=FM{|EX{P1+Ku|1& z9t1%gGERet+mK}%M5u=Bp|Ciqj$@N4EJkVox@BijJkUl+RSll#UT`1^WTe{p&$j65 zop{l_Hi!uG*>F4R zJDxP9Wj%;OJZBKLg}exHX=Xo*?4V`6h;gRJT*5uZWByOe#fi9-xA{ilI0w#G+eEm< z7mN(qEKb~niru!5A+&}Pq-g012ZUKmyN@{_nwEWC?c;<%KAT>EsjWyz-)my;0Otil<+7Cvpd%Mfzdv(0?SRU|&-suCY^TS1{LM|KZSppw&*`Uu7_>eoAv(0?SWs?AwP#$vG zO-h!)hg>ljnfZ`AgBLgRA$O+4hum2bA95E-e8^oa@gaAK#E0Ca5+8DxNqoq~m?F1# zom3uj*GqiJ-5~KHcca9I++7kM?cFW$A$O0&hupmqA9D9ee8@c}@geuP#E0Ay5+8C; zN_@yYCGjCwuu1)#_+_Iu3)n1Rvw+P4HVfD+V6%YD0&ijg{l2KA%lrsLO0#w+7n_l@ z=mNryY$mS9YzD5!Z04=UY{p&rPCET*8xES8|hEQX-REOv10 z`|j=aExIWJI>9uLCC4w=+w$y5UV4krm5^oufdS0|8v6%?G-}Yi1*uLA!y;SZVF_Wq zA`%F3r+_UD!y}_3VtPf!#t;XV&x2GZVrMWI_;2v&%|{V}fsOHL3}D1&n@a`Kq-uX@3!<2;tZV_ZAxogpC_879d0&1#1lo1KH{6r6ZF`1nc3>8GTR+I9 za{C<<*Gqo2rIvCEK0bt(PhgZUAh1DQ`87-67ZiM4DBx$dem7B%ov7bVU`Iip34PzT ze(0B1)bA|nbr3W^r=cc+4{~`)ZtA3zniErM77_gPb`h4$<*)PdHBIV)cn|u2Nl_kq zU`udpE}zNzLp}NNN%>T!&BPQBjuXY+(nP2KtbXi?u31w3K&!tfua5{jne`@8Li;3k z`rc}h;i}3dIJU$`={Kz(aoO^ zkGE9s+xA=U=(y0%s^&4yf{qqc6?-+4KH>RlD5r@c-&@e-f+mXcqysM(D#|Yu`N4u7 z6f{YcS4H_mQC?H#)qnJ}lf5nmt&oY6EPlBpo1L?YF@t1Lc|Be)VHSng=WnFA)U5u& zdX0^)T-A4@@i9pQ7FXA6bV9-`j;_x?MRA$ud-hXePMA@GQLs6Farz6JIP20e+*}j9 zpF?z+^ru63Mm!rg83nsg%Z&N>xEi>a(E{TWiokeBIq?0S_{AyS5)H2Wu1LHc0{M83 z6&O1qfH9*6#(i{P+(!q-P6%M^gaF2>2^gy;V9bw!U$@%|kA$K;9s>o&bDzL?(i8aL z$SZxVa}}M)c8Jo%u2Yu&6|3kIYf|p>=7_)V`h3^x)8@FgA3U4isv2=(7wX93T)KO` z9sEbZ_|plvtC~SQ8~E2*;9d7P&f7NMw9X!7(9is2{;PVFTcIgGiwy{f7GVqm#ux;Q zF$fqQ0gOH=oA7;%LCAk$YS)&ZVe5NSK1*1mRhTi*7`{*ZyeI?8`8R_%c z<{~|2n~HYEj*O#(5O_#zhuJkWhV75d0{@W(5UF|{#R}b`FqP-X^_Y(n`9x6gWy))L z`7g@J$IJ3$EZ2m5)CUY|F0YDTf}%Vq^e~qfKadZKa%L=*gWT6q`FL3_4gO75Ss3pC@x1=}tmh=YSl3ehirzMJVmMHYJM4_i83Oy}ROLD<~M|S|;`cA4% zAKug0>|nD4vw+PHY<^(#1KV|A>j!N8fUO^}rhY)b4`Zy@nd{A+lsIm;qc6P4%6+Y6 zM9rY*-{TD2Bk`-A5nrvD5znfX5n~?)BI?+Of%>oy0~q@-fUyq)82d1Qu@3_n`!Im9 z4+9wcFrW|iVIUv-Fo3ZS0~q@-fUyq)82d1Qu@3_n`!Im94+HvO9|rQV4+HpJ-vh&2 zZy*o*Fo3ZS0~q@-fUyq)IN@qWe=s`p-4DoLap~e(R>JJj-`LR?U}dF8m*xtX-@#?` z=!|E?=Fu7Lky($f-L7k`jNcca&j6f>%znyX#JxlFE9vRG?r}jFjH2fx>A7M3aCr!y zwUPfN3!qQ?Pxm2Ybss$X!rqnoS;Uomybs1+f~>s}zl)wS{f+hA5@v5@(C6=!@IH#Y zx_-A@ud$Dq_n2rNOOF58D8Vg7Kk6msJ;0dv0At<*jCl_*<~_ie_W)zw1B`hOFy=kL znD+o<-UEzz4>0CEz?k;{W8MRdc@Ho~31G~7fHCg@#=Hj@^B!Q#dw?+WnX9Nx z!|2a{!OkK=M)dn=qSh>&kN@<1o}Xjigxkq_A|%?gT=S=>nqC*?Xf*M{c>avkhdnvp z6ve3;%lYrHyrzfMMUva{_7g$zK5k%8@KsR`G^xF?p9Ov*>J^lTXHG2ft<64achNpU zNhHTfodiELBMueiKux0jWl=sEz0Aoc0j1>tmRwTe!kpN>w0*n<2 zFjgec2P+cfV?_dt6$vm_B*0ja0Aoc0j1>tmRwTe!kw71;NRW>e2{2Y9z*vz0V?_dt z6$vm_B#?&{3G%Tb0mcj#7%LKBtVn>dA_2yV1Q?+v;KNH&qnT@%J^CM8V?zU*NApV) z_ON+$Mti)@qgRwTAEfW2KNK$@cHq@$cgxCW>AeNqB?HC11tHKJH3CK-4U9e-7=1J_ z`evF#2d<^wGfRqk+*!1EY@yhDQUVj|N5`4U9e- z7=1J_`evF#2d<^wGfRqk+*!1H+?%(MJQLj|P6) z_WRHE9{y5X^dGw)jJOK=Y+!geFgzR>9u5o-2Zo0O!^45$;lS{4V0bt%JRBGv4h#g!-3)9!0>QjcsMXT92g!B3=ao}hXcdIf#Ko6@Ni&wI50dM7#9u5o-2k!Ud+WmSD?=DCF|1l58ybK-=3=ao}hXcdIaW@Daj(m7H zFgzR>9u5o-2Zo0O!^45$;lS{4=mQT&K0F*49u5o-2Zo0O!^45$;lS{4V0bw6frld> y9u5o-2Zo0O!^45$;lS{4$b*L?A07@24+n;a1H;3C;o-pWaA0^i@Y6`g&i@a(rJ%O} literal 0 HcmV?d00001 diff --git a/tests/unit/data/verify_golden_stations_tolerant.nc b/tests/unit/data/verify_golden_stations_tolerant.nc new file mode 100644 index 0000000000000000000000000000000000000000..8242a88de2e9524d998ee845d1615d45b6bc000f GIT binary patch literal 57675 zcmeHQ2|!gv+nx&-Lj}c@Z~;X8G{p^Z&-5Ho$X!BlO#$Jm7y=hCHFG!ja<^PFtt`td z)2!6WT+^)7)UxGM)Rf#(Gu!*;Ju}a}hf&h&2X^acO6ETA%$fH+bDr7f%$eG@O|)-` ziX}X}y&2q}(sej4B;hh(#_x+myR>Z+>*2HBBh-+2Hj_c!gT;?C(l1q6;1V9pLcjb2 znU#B<7KX%C-zSRVd9qT>o1(^dZXJ`Eo?*+(PO)XNvOI`1dhmab*fO(H($iW*)K0Ty zCnXPP(Xei_x(#X%N*>UtMg6+<8`f>i%J68{yH9$p?q!(o+uLF(N8uh8i-nRT2j8ZL_m&A40j!t|6RADD?C{_=EMK- zyEbmAVe9nt%;c1`#B6@=eq0UZ9b+Me!K98HeERKEi?AmvCNi&e$5vgt8kt76k}+-C zb&BfTwOyCa@m*WD>KMg*s*v21d5H@;x9SvS#J@|%&q^PfnPfB4^)w!ek9a7;yf!IA zQj+r;#;+4Dt>wD$B9{Y9x$G%J9lj6QPq!vAQD6?h?z*0`$O}79d8bG2A(3f`rc9&9L=L{|nO2ftElO|kZ;Yha9WA|7AY9M1 z8aL)yDV~_%SzzY*K{=N{)z z<5z%!%PY8of@Rz3kXDkvU{7{hd{X+*wCq4BNB(7}XD6nLAR{D5xo&{GPWJa4_ctIz zTqFtSkY`J|EcOoh!_r{hfiWEHx+rmNMII*!Sx1o%mk+eAB3~nkgN`cFhu-^?_|V&2 zkwfqOiX3`dD01j+tH`0Zog#P^JBb>3fu9x## znz3zCj^~B@+OEiPUpo{z?rWza$9?Tm?ohZQ;W9#`bh zo1@5~_k<#c-jj+PdQT~G=>1-iL+>R;4!xHZIrLspigp2~_G_Jk>NlHcd z*`UZ_&zp)I_H0$;(EEWRhu&?99D27aa_HTm$f0+ql%u@vQt)mCf2iO+3f?QRRVHY& z$2X#%A6+kRLh^mr%j=W;ejn3XL-nQ*JpXE-*}dX<_|phzohiE?3A>DwzQ-Xu1)cXZ1gNWGC_4ZTDp^~*Fte6=$F65>#Er9-MtYRB1JW~XNr_q6@tL*(woIPGw)pInAvPWo zU2H{#=&oI(SVPQT@oUadu33!50{&xf_m6t9Bg?Lp!30_jcJG&5&K>4)$>3XF#1M~U z?9vZ5B0W}* zmbyfFVrb9$wZ>9{E?=1XQ<{7nLNjBn&O|&F9osgB5ji}J6;v67CEpxb@1ayXCV4QJ zmcG$IA=$t1dD}?!h~Y%iA-sB+Bd0ViETVGyw~^i2w&AOtJYGp2mRzI7YlsgmZw-F! zuo0^m#VT#YvdQTe3q!0LW4)uIyYjRc?m}j{V+LFr!e_o?*l-rYp|Hs(L4;sN6Mp*8 zAyp!~b&O_SQ9eIf*m$lygvF(uxr_Pp8w+9&wOiuLPV;N4u}>B~@5c^uX&~#l`*ImJ zo$Id5?#UT=H(TXJ(lh(p__I_)x;W=tMHb7YLF}mgkuc`ZA9FB!dS3HzR-db=#J(Lf zvI*P9BUEN(*EMOusD~KH=1xD+ighwlc)xe^NaoMu2C^^aeILm>8kbecxYCBD8`|FO z{8Tg>$;}C7dw-nLkqtB=*za51i}~=xstVx7?wB~&1na4-byDj)KPc#h|zeF=sg~} zreUfPK@i2UTp>6XfNKSJd~L(hAt~}rIbg2K=XVhS2{PQ}Dmm^%_sycO5n<|U@Iq!d z_^x;BX=o$*8osims(sHw?Q76Vpv!#?AHLpnqa7K3JTYEdv6ybXVxml|w~P&nO4)zI z?JcB#YPpX4qT9B~+drAn!1II!*~*HEtuq&SY;qPWim{b;e_&W?v7$H^6FVF*uC!QD zuzj<#O{@4AUs|kK!7P<`;Ta1X7AwkrN%NS6I@^Xhs0sq|>2_nFY_ZCjM0R(#F5Ti| zVxw9+^9i$$7Hd_JF7s>CCcb{SX(1dF02V9nejBc|lG#V(wH7Nq4RXbeTZ(Mgs_W5v zg~VblBYrS%Bdk-`sL;j;LPaHmT1D+>iXc>4BE{BSVq=FlM;L0bdyEuKnG~O$*gw_w zdrPS?ACJl^6@waV{{a!@wgwu+KHjDeZTT(|3}aEBD+$(TywBpvit`>4CsqbAg6RMqS z_P_{-4>|7*7!DZWq6#S!qR9@N&Kj0v4gciP&6~NzwLgtP_dJxEGwu?Mzs{|J*8B&q zMm^)wCOO_eF&uoftzQ@`kmXvG*=IB{j8sV(L&EyU{JtcnvoW-%!DR z8sjsAQuXA`tRJQ1+-%peUDs~BCFk`f|3{4>^=v66*ahE6QFeQIl_k225>xpxhf~MV zjA<#{j46LzZZyvA>ioL76*|w#k&{U3LBpVxyv?Vxu~@j*2&)L1X^V=@}#^W+!IZvYihjna?z& z(SB0qP2oa&$PG7I!Xj|#*X^f@lSMx6PS2-HE#7JvYAo`*RE}s2M-D#jot_GD4rPhL zoAW+gL+!ftC-dbz7 z1S43Gmqf@a9|GK&<58c$oDTIZ%;D6=8!Dt8PUt+CDt~i0wS;C&(`9B%yL8N$cJr7q z%_SNg5<^I10;h{xvZY$bKKr$}S{L8`Jmaty-TK|;2a0{Bv16@#&b)d=%g#-FV&(Kh z8tZ4RbPz*Y(_cv1YKe4rd~KOWe2l#lvh znq5Cq#}++gpH@Dr?@GFq-Le->bg7#r$%5THA6PHgxgRvzyL}JnT_MXF>Bz_6=eB zpt$=R?0=(-zIlD~jG4W?^rQYWQ_6hjrh71y0dBxFK748y@zqgU#`YZmw(Y3~%a3510z?kX)UcLU*CHJgVuTd?*oiU@q zeY}-bRIAEbwFl1>Ghurwm}Xea;k^~yhhsiLxP9vsK>+49%_^p6uPVz=H>ZGTkV{ip zD9VrnyL#R2O`)%lHP&)z-p!wpPvK4dS$Vb; z;8gl6^vGP9O2^bBd~@-$bKA@~Q^8pZCLWmM4^{9Z3Z`{Xm#Nt->hgPLt~XLjoXZLR z9DCg0I>?Cj2N>-SFp2^&CbsEGl6xv=X@q;WPrKyAWKsLX9~0^(GhxSbxfxGT@I(ba zuHZ==Gq!2l2MgV;eG;uqc~Pmm&_NN`dCKLiyscBGe#Y&xMZ#mwZ8H!s_R5Uyyb-rt zn|6ZAtVk_aG~;ZJjON4FqHeUyfTroGMYw(e-L+anRBfLv3F5p*Aw{7D3(7|dV-UW% z^J53|ZdKrxwD*0kuf8c?L$)1 zK9xd;`u-WCJ>vOlrKJaQK8|p{#y>)gK6$V#@aEuF zQ5D*VYb;ht2j)zj*7^@qw~RZBHCX&$6&>1`*s@td+$NQ(}}irzJ@&)({!) zG?nMisQyA~u~xY*nCFpwFoHv`4^AB3(1zgu$l#K4BL(GB(;*10CPPceP@YHW)HwaM z6J^-_wg9YXB_q#%7L6=5BGs#H@eM`tMObL(}^C%j=qLU zCR>YJn0OJ_ouSP^+zk|-vwOgv@AOPSKJU)(+Jg3ov%syjwk2#@@{cO>jzW)IcXUEz z4lD?7n)3Y%mTtKZ*!8ntt@`q-(RMv<&-e!_#n|=p9l!fIIMc4L>$`A7lYVx+kAHOH zAIWz8nZ2E#Il{vm|Bx`ZX1ZOk99K24+=F)gzky>%J=@=|uRLGAN_mw?w+4ZH}eycfty)X!!1Xd9B36n`mckGA(k+uEof zH)wHj{SPFDebe%jO|~kJ9htA$zFK{0NAuTQ-cqpgWLWaK-UqcVbxY+o%Q~o~+m0q? zr0mufp9q}Ot;5@zmb)f8xzaZ6y^6DHSO00R%jsWgytMNhxkt3X1HY8^<9^({d*7W= z`+#<8-HSgR^#4$s`g-W1_dnaGJ#n(`bH|SCcRBsySI5=zjQUPXk7_^ebL~6re9+m6 z505^lwI6)8^}JW8>Eiy%BcAk;Pvs;=K0&Tt7ouE;@ETGMA9~!7A9{Sp&U4y-d^c?C z?RQ!`GwCT?Lez2X=T|R&x}n-hZDM=-sQUUzZGpbSYv(s7vi$$c4iAo0fB7^Pf2uea6+cNum2K`kYx`pS;+}LqC1kyZgQ{5B-m= zFZ`C4;Gqv%HtU1>gMZP!_un$^gPMOj{2om!Rt{+UhXy}-{r>NMEdRR(-&ib>D<*V6 z#gxX;pREoouh;5v?di*p2k6_rJusledw1(G|5+Z?<09Xx@Wz?;mNLHj)i$F>&K%{V zi+gsp95{2zJqH29S_un_1_3-D(3v*6agr))x!jCjQ}F8wUZ>zUIA&}|=LKgRgMe-% z`TO5}5S{gu+nl`AWG6rGwkY|XuFeStOwl8UvX)$^U+_6}oZp8zhjJ8rLcu2$ ze2U|V&yU(5sy}|JyIVB~BOl`Cako%-w=nBdNOQA%$!}$B{d}g2_zKM_nl7@ewXg@N zjBE}q_Sol6&b^~FYd!_KgkOR`elr6~Dws|hFo*jpm=4-AhtmOjW-JcYo9|Wn3+El} zPZAFV%~HX>-$zH-nJu8B>&%#rtTSUes?LnX5p_qrLNlBehb9-M#cdS}=+HZ}MeP;b zLBSo3$m?dFI4)mv;>KRAIQSk>D%-kQseEIA$V+2SAOHFHtxmbxc!eyrjweoFFW3Ay zaVH7FbA3su@zU9o%?V!pUlR8L|9 z7DpUze7HXm0nOar2J2S|+9tUfZMP(f%tR;Z@1t!S%Ey9uubQ!bFK%DuH~*?5`v>)y zY`k*|hBs0PSzx*gw$2YzxK&HRwG~`P!E_wz+m_s5XW1#4nD=oPT%Vb-SI>5yc9YJ) zlnK_DJUccI89gh;5`BKJ;1N^E@#m&zVo}=X`)N@g&cB8d^*S`P5FSyZ;r7}D*5B^U*&BL{( zLsYmQzYxUtMUkTv$yIIH-QPWX;%nss>R@5011Cs9;skD$=Bt* zF+U%DWv#Jeh!jU<4qE4)eSj=E4ER^|0UizRd)C?R{hCgXz`2#CV2i_6g)DSjF#q2{ zHaaHQjOln_GZx1J*ZO_(GIN1Aq812OJNCmZ$1j-^Y%+>$V(gfq`!@Kz&i9hPIl*CL zM1l#zbG?AD@B2UG=%rz6B0A0IvqZAM;~&k7lj<%iD+<(f%k965A@f2!53T^M4(3&9 z4uc%zYW&;TxAS|ZNb`X3Zl&EQyf-n7oLNlf#*3EdDf~OmSxsNn^(`oT))#3u3!+KV+V1}E`du$!?sd_5!XAn$*ojHyp7BE0n|_iVJJ{s-@jy8G(o3E#SBqtrom?EXbnB4eK?&f#q>-YP{ z*g89`$%&;;Hqg(GDC;?ReFJ^+k1bz%&99+;ZE?#c1JWDlO}d_LTFSG5KCXFmeB^`( zm(vez-FN-($rrVa{oZZVV$TI_NLon73-*iJbG1B!r@W#50HtpB8oTHBpr(`HY) z-9y>DrCs9g`5;kZ?3DFGIDACR$H#RD28PdMIu4A8(hrg0>dsr(ft|R(HDrA3h`kxU z_2(-L`sc=?cpia(;vBn;GFj;(8;|!u*5!rm%b$`I(R;?pu!+{#c2}O1xa+_a%0L|1YoRZrL2& z@bAkM{$F5`>n_^QdI`JX>k`8sD=~bbchc$I zdd%HgtTkk4dAV44-JA6oI$5lkM?^yOUVY{WalOSFrr5yWst=fmz}g5DUk{+y;``CD zE8JPEaaQI=Zh_3lSJD0|P|hb+?5Ppv780 zZ!&aOSKdT25i(HP#&2r>dW1E|8^+HGuswyahDuoC@U@cf7HcCJ=F_@WOjLZQxVofr z!6sZPxAn4;2$j06!{DRrBH+`u{jV9-adBCv>=9@_Kgww_Wb9~y; zVa8M*%(w@~oBw?2viXjV$ib`5!^jo>xrAT=l{j)qIj0Xg*A{Rp1ZD|v3afFV3bKHH_gIgUTLq#KU#Bv!IyzEGXY+hk}!fzf~hqX7lJ8AAq*Il04K z>y#s7#a$%Go4JK&Ono}no-s%YAN?dEDG`>)TX4^|U>X(FAvCffJArBVV#YLJF=M`+ zCl97E%P%z>MRAY$c=oLx*{xk?Rt$%v$?s%NszSm1OQW{=;z!8c^Be=4(3go>}nRt5C2&?j5}1Gt(e70k_)Za0TP46_CD)3o;@G=eud!TgJO z(#3FCi+>zcI`N5EjPm98UL8!Y$MuFX4D;EW{MsOXPtu96-C``yF7nVohT%Ia$!nJ& z$G3S|4GzN?M&w#XnZTlwaoX@d}=x;E5bpth3+@ z_nh4~bRBp}HXJ&a!0228ule+?itfswCPfUo@fQxzb)YK4%owUO%($6??^AGdj%%$v z_p+a0fV;hPc5NEBkGplk`}?*M6LM+5r_oix`)&Rcj4U*@GbFLIDy6wtB526B9BnC8I9L5@!i(_9!i z$nl9`^iv>*UUXO>hhFqtActO>$sz|i^r9~VIrP>L(H-%j7ySsthhB6jActP`EFgzo zbTJ@@UbO0vLob?l$f1{tEIG)bmx?Yq$f1{tFgeJf*BDcorB~-uDza^uDjip?9w$hu)7AIrQ#Rvm6y<>1 z_NVa$6T@L+BEqp55f~c~fwAcj7`qICb6b3sW2|d>MXqU*!Nh9HB)Jk55F`gX<(o$~ z$&nTmdH>EF03!;Ukq0NHKklw|WDhwBlpr5*Q60%M!d%DAR4`3ynZs!^%ZzCf%ZzFA zYESn|V$K?jL!IWq1Ig{!bT!3Y8S>OPVPk^u{tYW71-t2nrhCm#oThorn5K8ln5K2j zn5J_p)Z2KKd(O31w1QyF&g%Bw?#hs6I}FoVg7AtMLiu2R!u=JTq+nWjFvq6_2Q#LH zhCk9?`^jAiFl+)hkh=uq@85VXL%}Vk7HkzkGxz2PN;CFmOf&UnOf&Rm{65E9Qdf&T zKj;{vqYaFXHZZz1!02cLql)1>gbas93E}8y1EZr2jE*+&mzx&4EkidJb#wlDy9ceM znHfiGXlA@z!Otpq1;^em)HuXF?mUIZKrv9H?xa3BGUy-fdF$c^^79=h9msf0`jF>$ z9pq|!Tn9Pg<>#yMaUJCOT?aYh!^d^-5idU<@o+tS$n(1na>Rp=>*1G|`I(=ObPyld z?9W|sK`T6;!3weImL5_I&`G}XFpWk(myZQ!rQG8~&qI@XI z2g-rJn+GNIIKvg!#q%g3PJUQj!hvrb?s1cF)O#Y;8~ruE!5Ujf#E=`pn9nO7d^3JQ zo|;xongIW(-1_e!<)Lz?Jn#;=c@B8E+yn=VQ`gYLDIX~kLO4F{4UDh+0w0o_*nrRU z%6mP&PK0oL&k6W&&%E(*a2xD7-7jyryT8wkK@9E>gBV~8Vt_G-0Y(oK_m4pg!ZC;e z#vle5gBW0}769WE0^lFs%A3D0Pw2NS?+F$Bm5QmUR5gUK2h)6qv-kMe zt|oh=PrWgWMtZ^Q%kSP8$Nc%v3T)`t<0dol7M1_Ut!A(XdE7vzvCsvK4$ZH`%EWp; z&1ecEj9H%wT+GCqUV+O#d69`1sJr;L2t9s zJcTMOs9Nn^@^z?c0aJI$SCWdoasPfMUMTDS>e<8c^{Fetmku-WLf7EsLyyYWrb1HL z&+@?ZN%FX$Au#Ls2A8jl!+E`$; zvA}3!F)@WU7U5`PfzieSqm2bd8w-p!78q?TFxpsPw6U-UZ7jmk#scHa6kwd00*o_L zfbot7Fy7Gs#yc9oct->F;2jNw;~foPw6VZwV}Vg)0i%rtMjH!#Xk!tMHWnCdEHK(w zV6?HoXk&rVNCNjCou6qk+L_U@#gOj0Ogyfx&2CFd7() z1_q;n!DwJG8W@ZQ2BU$&XkaiJ7>ou6qk+L_U@#gOj0OgyfziWaa$f?UVBMS8DdMzN=Z*+-qb%4)`l<-i^WpSsA=iI z0pU;wBgGRe#dxxlD%fl=oIqs|3JoePXQ7Z`OeFzQ@j)VaW@bAeIk0)x@O zsB?i)=K`b71xB3buKXKTwv6>z^HS9QRf1q&ILxD3yeA!7Re#dxxipF zFzQ@j)VaV-{~PRSX83ipt#_>7jJXjo92g7-2E&2DaMW>NIKshjU@#mQ3rS`42A=P;lN-xFc=ODh697)z~%OpKWb+9b#onejN#~Kf#JYlI4~Fv42END z4Td8e3 members[:, sel]).sum(axis=0) # no ties + np.testing.assert_allclose( + _sel(result, "T_2M.rank_hist", **at).values, + np.bincount(ranks, minlength=n_members + 1), + ) + + +def test_crps_matches_scores_library(): + import scores + + rng = np.random.default_rng(2) + members = rng.gamma(2.0, size=(4, 30)) + truth = rng.gamma(2.0, size=30) + fcst, obs = _datasets(members, truth) + + result = _run(fcst, obs) + + expected = scores.probability.crps_for_ensemble( + fcst["T_2M"], obs["T_2M"], "member", method="fair" + ) + assert _sel(result, "T_2M.CRPS").item() == pytest.approx(expected.item()) + + +def test_identical_members_crps_equals_mae(): + """With identical members the fair CRPS is the MAE and the spread is 0.""" + rng = np.random.default_rng(3) + single = rng.normal(size=25) + truth = rng.normal(size=25) + fcst, obs = _datasets(np.stack([single] * 3), truth) + + result = _run(fcst, obs) + det = verify(fcst.mean("member"), obs, "ens", "truth", [_ALL], num_workers=1) + + assert _sel(result, "T_2M.CRPS").item() == pytest.approx( + _sel(det, "T_2M.MAE").item() + ) + assert _sel(result, "T_2M.ENS_VAR").item() == pytest.approx(0.0, abs=1e-12) + assert _sel(result, "T_2M.MEMBER_MSE").item() == pytest.approx( + _sel(det, "T_2M.MSE").item() + ) + + +# --------------------------------------------------------------------------- +# Spread-skill ratio and rank histogram shape +# --------------------------------------------------------------------------- + + +def _calibration_case(spread_scale, n=20000, n_members=10, seed=4): + """Truth ~ N(mu, 1); members ~ N(mu, spread_scale**2), mu varies per point.""" + rng = np.random.default_rng(seed) + mu = rng.normal(scale=3.0, size=n) + truth = mu + rng.normal(size=n) + members = mu + spread_scale * rng.normal(size=(n_members, n)) + fcst, obs = _datasets(members, truth) + ens = _run(fcst, obs) + det = verify(fcst.mean("member"), obs, "ens", "truth", [_ALL], num_workers=1) + # The ensemble scores merge into verify()'s output of the ensemble mean + merged = xr.merge([det, ens], compat="no_conflicts", join="outer") + n_members = ens.attrs["n_members"] + ssr = np.sqrt( + _sel(merged, "T_2M.ENS_VAR").item() / _sel(merged, "T_2M.MSE").item() + ) * np.sqrt((n_members + 1) / n_members) + hist = _sel(merged, "T_2M.rank_hist").values + return ssr, hist / hist.sum() + + +def test_perfect_ensemble_ssr_one_flat_rank_histogram(): + ssr, hist = _calibration_case(1.0) + assert ssr == pytest.approx(1.0, abs=0.03) + np.testing.assert_allclose(hist, 1 / 11, atol=0.01) + + +def test_underdispersive_ensemble_ssr_below_one_u_shaped_histogram(): + ssr, hist = _calibration_case(0.5) + assert ssr < 0.7 + assert hist[0] > 2 / 11 and hist[-1] > 2 / 11 + assert hist[5] < 1 / 11 + + +def test_overdispersive_ensemble_ssr_above_one_dome_histogram(): + ssr, hist = _calibration_case(2.0) + assert ssr > 1.5 + assert hist[0] < 0.5 / 11 and hist[-1] < 0.5 / 11 + assert hist[5] > 1 / 11 + + +# --------------------------------------------------------------------------- +# Ties, missing values +# --------------------------------------------------------------------------- + + +def test_ties_are_split_evenly_between_ranks(): + """Truth 0 mm with members (0, 0, 1, 2): ranks 0, 1 and 2 share each point.""" + n = 9 + members = np.tile(np.array([0.0, 0.0, 1.0, 2.0])[:, None], (1, n)) + fcst, obs = _datasets(members, np.zeros(n), param="TOT_PREC") + hist = _sel(_run(fcst, obs), "TOT_PREC.rank_hist").values + np.testing.assert_allclose(hist, [3, 3, 3, 0, 0], atol=1e-12) + + +def test_all_tied_spreads_over_all_ranks(): + """Dry everywhere: the truth ties every member, so all N+1 ranks get 1/(N+1).""" + n = 10 + fcst, obs = _datasets(np.zeros((4, n)), np.zeros(n), param="TOT_PREC") + result = _run(fcst, obs) + np.testing.assert_allclose( + _sel(result, "TOT_PREC.rank_hist").values, np.full(5, n / 5), atol=1e-12 + ) + assert _sel(result, "TOT_PREC.CRPS").item() == 0.0 + + +def test_nan_member_makes_point_missing(): + rng = np.random.default_rng(5) + n = 20 + members = rng.normal(size=(4, n)) + truth = rng.normal(size=n) + members[2, 7] = np.nan # one member missing at one point + truth[3] = np.nan # truth missing elsewhere: not a forecast gap + fcst, obs = _datasets(members, truth) + thresholds = {"T_2M": {"gt": [0.0]}} + + strict = _run(fcst, obs, threshold_dict=thresholds) + for var in ("CRPS", "ENS_VAR", "MEMBER_MSE", "rank_hist", "prob_table"): + assert np.isnan(_sel(strict, f"T_2M.{var}").values).all(), var + + tolerant = _run(fcst, obs, threshold_dict=thresholds, max_missing_fraction=0.1) + keep = np.ones(n, bool) + keep[[3, 7]] = False + assert _sel(tolerant, "T_2M.CRPS").item() == pytest.approx( + _fair_crps_brute(members[:, keep], truth[keep]).mean() + ) + assert _sel(tolerant, "T_2M.rank_hist").sum().item() == keep.sum() + assert _sel(tolerant, "T_2M.prob_table").sum().item() == keep.sum() + + +def test_truth_gaps_do_not_mask_scores(): + rng = np.random.default_rng(6) + truth = rng.normal(size=10) + truth[:5] = np.nan + fcst, obs = _datasets(rng.normal(size=(3, 10)), truth) + result = _run(fcst, obs) + assert not np.isnan(_sel(result, "T_2M.CRPS").item()) + assert _sel(result, "T_2M.rank_hist").sum().item() == 5 + + +# --------------------------------------------------------------------------- +# Probability tables +# --------------------------------------------------------------------------- + + +def test_prob_table_counts_and_threshold_dim_match_verify(): + rng = np.random.default_rng(7) + n, n_members = 50, 4 + members = rng.gamma(0.8, size=(n_members, n)) + truth = rng.gamma(0.8, size=n) + fcst, obs = _datasets(members, truth, param="TOT_PREC") + thresholds = {"TOT_PREC": {"gt": [0.1, 1.0], "lt": [0.5]}} + holdout = ["S00002", "S00020", "S00041"] + + result = _run( + fcst, + obs, + regions=[_ALL, _WEST], + threshold_dict=thresholds, + holdout_stations=holdout, + ) + det = verify( + fcst.mean("member"), + obs, + "ens", + "truth", + [_ALL, _WEST], + threshold_dict=thresholds, + holdout_stations=holdout, + num_workers=1, + ) + + table = result["TOT_PREC.prob_table"] + assert table.dims == ( + "forecast_reference_time", + "source", + "region", + "station_group", + "TOT_PREC.threshold", + "n_members_yes", + "obs", + ) + assert list(table["TOT_PREC.threshold"].values) == list( + det["TOT_PREC.threshold"].values + ) + assert list(table["n_members_yes"].values) == list(range(n_members + 1)) + assert list(table["obs"].values) == [0, 1] + # Every valid point is counted once per threshold + totals = table.sum(["n_members_yes", "obs"]) + n_valid = det["TOT_PREC.contingency_table"].sel(contingency="total_count") + np.testing.assert_array_equal( + totals.squeeze("forecast_reference_time").sel(source="ens").values, + n_valid.squeeze("forecast_reference_time") + .sel(source="ens") + .transpose(*totals.dims[2:]) + .values, + ) + # Brute force for one threshold, region "west", all stations + west = fcst["longitude"].values < 7.5 + k = (members[:, west] > 1.0).sum(axis=0) + o = (truth[west] > 1.0).astype(int) + expected = np.zeros((n_members + 1, 2)) + np.add.at(expected, (k, o), 1) + np.testing.assert_array_equal( + _sel( + result, + "TOT_PREC.prob_table", + region="west", + **{"TOT_PREC.threshold": "gt_1p0"}, + ).values, + expected, + ) + + +# --------------------------------------------------------------------------- +# Grids, special parameters, dims +# --------------------------------------------------------------------------- + + +def test_gridded_truth(): + """INCA-like (y, x) truth: works and has the single 'all' station group.""" + ny, nx, n_members = 4, 5, 3 + rng = np.random.default_rng(8) + lat, lon = np.meshgrid( + np.linspace(46.0, 47.0, ny), np.linspace(5.0, 10.0, nx), indexing="ij" + ) + coords = { + "y": np.arange(ny), + "x": np.arange(nx), + "latitude": (("y", "x"), lat), + "longitude": (("y", "x"), lon), + "forecast_reference_time": _FRT, + "step": [0, 6], + } + members = rng.normal(size=(n_members, 2, ny, nx)) + truth = rng.normal(size=(2, ny, nx)) + fcst, obs = _datasets(members, truth, coords=coords, dims=("step", "y", "x")) + + result = _run( + fcst.chunk({"step": 1, "member": 1}), + obs, + regions=[_ALL, _WEST], + threshold_dict={"T_2M": {"gt": [0.0]}}, + ) + + assert list(result["station_group"].values) == ["all"] + west = lon < 7.5 + for step in range(2): + crps = _fair_crps_brute(members[:, step], truth[step]) + assert _sel(result, "T_2M.CRPS", step=step * 6).item() == pytest.approx( + crps.mean() + ) + assert _sel( + result, "T_2M.CRPS", step=step * 6, region="west" + ).item() == pytest.approx(crps[west].mean()) + assert (_sel(result, "T_2M.rank_hist").sum("rank") == ny * nx).all() + assert ( + _sel(result, "T_2M.prob_table", region="west").sum(["n_members_yes", "obs"]) + == west.sum() + ).all() + with pytest.raises(ValueError, match="Holdout stations"): + _run(fcst, obs, holdout_stations=["CHM"]) + + +def test_wind_direction_and_precipitation_type(): + """No CRPS/spread for angles and class codes; PTYPE still gets prob tables.""" + n = 12 + rng = np.random.default_rng(9) + codes = np.array([0, 1, 5, 7]) + pt_members = codes[rng.integers(0, 4, size=(4, n))].astype(float) + pt_truth = codes[rng.integers(0, 4, size=n)].astype(float) + dd_members = rng.uniform(0, 360, size=(4, n)) + coords = _stations(n) + fcst = xr.Dataset( + { + "PTYPE": (("member", "values"), pt_members), + "DD_10M": (("member", "values"), dd_members), + }, + coords={**coords, "member": np.arange(1, 5)}, + ) + obs = xr.Dataset( + { + "PTYPE": ("values", pt_truth), + "DD_10M": ("values", rng.uniform(0, 360, size=n)), + }, + coords=coords, + ) + + result = _run(fcst, obs, threshold_dict={"PTYPE": {"eq": [0, 1, 5, 7]}}) + + for param in ("PTYPE", "DD_10M"): + for var in ("CRPS", "ENS_VAR", "MEMBER_MSE", "rank_hist"): + assert np.isnan(_sel(result, f"{param}.{var}").values).all() + assert "DD_10M.prob_table" not in result + rain = _sel(result, "PTYPE.prob_table", **{"PTYPE.threshold": "eq_1"}) + expected = np.zeros((5, 2)) + np.add.at(expected, ((pt_members == 1).sum(0), (pt_truth == 1).astype(int)), 1) + np.testing.assert_array_equal(rain.values, expected) + + +def test_output_dims_match_verify_continuous_scores(): + rng = np.random.default_rng(10) + coords = {**_stations(8), "step": [0, 6, 12]} + fcst, obs = _datasets( + rng.normal(size=(3, 3, 8)), + rng.normal(size=(3, 8)), + coords=coords, + dims=("step", "values"), + ) + result = _run(fcst.chunk({"step": 1}), obs, holdout_stations=["S00001"]) + det = verify( + fcst.mean("member"), + obs, + "ens", + "truth", + [_ALL], + holdout_stations=["S00001"], + num_workers=1, + ) + + assert result.attrs["n_members"] == 3 + for var in ("CRPS", "ENS_VAR", "MEMBER_MSE"): + assert result[f"T_2M.{var}"].dims == det["T_2M.MAE"].dims + assert result["T_2M.rank_hist"].dims == (*det["T_2M.MAE"].dims, "rank") + assert list(result["rank"].values) == [0, 1, 2, 3] + assert "member" not in result.dims + merged = xr.merge([det, result], compat="no_conflicts", join="outer") + assert merged.sizes["source"] == 2 # ens and truth, no new source + + +def test_needs_a_member_dim_with_two_members(): + fcst, obs = _datasets(np.zeros((1, 5)), np.zeros(5)) + with pytest.raises(ValueError, match="2 members"): + _run(fcst, obs) + with pytest.raises(ValueError, match="member"): + _run(fcst.isel(member=0), obs) diff --git a/tests/unit/test_verify_regression.py b/tests/unit/test_verify_regression.py new file mode 100644 index 00000000..34621d98 --- /dev/null +++ b/tests/unit/test_verify_regression.py @@ -0,0 +1,162 @@ +"""verify() output must not change when its internals are refactored. + +The golden files under ``tests/unit/data/`` were written by verify() before the +region/station-group mask construction and the missing-forecast rule were factored out +into helpers shared with verify_ensemble() (commit be3dc5c). Golden files rather than +loading the old code from git: the test does not depend on the git history, and any +later change of verify()'s output shows up here too. + +The inputs are built from closed-form expressions (no random numbers), so they do not +depend on the numpy version. To regenerate the files after an intended change of +verify()'s output: ``uv run python tests/unit/test_verify_regression.py``. +""" + +from pathlib import Path + +import numpy as np +import pytest +import xarray as xr + +from verification import verify + +DATA = Path(__file__).parent / "data" +_FRT = np.datetime64("2024-01-01T00", "ns") +_REGIONS = [ + {"type": "bbox", "name": "all", "bbox": [1.5, 16, 43, 49.5]}, + {"type": "bbox", "name": "west", "bbox": [1.5, 7.5, 43, 49.5]}, +] + + +def station_inputs() -> tuple[xr.Dataset, xr.Dataset]: + """12 stations, 3 steps; gaps in forecast and truth; circular, categorical, + truth-only and forecast-only parameters.""" + n, steps = 12, np.array([0, 6, 12]) + i = np.arange(n)[None, :] + s = np.arange(steps.size)[:, None] + coords = { + "values": [f"S{k:02d}" for k in range(n)], + "longitude": ("values", np.linspace(5.0, 10.0, n)), + "latitude": ("values", np.linspace(46.0, 47.0, n)), + "forecast_reference_time": _FRT, + "step": steps, + } + dims = ("step", "values") + t_obs = 280.0 + 3.0 * np.sin(0.7 * i + 0.3 * s) + t_fcst = t_obs + 0.8 * np.cos(1.3 * i - 0.5 * s) + 0.2 + t_fcst[1, 4] = np.nan # a forecast gap at one step + t_obs[2, 7] = np.nan # a truth gap + tp_obs = np.maximum(0.0, 2.0 * np.sin(0.9 * i + s)) + tp_fcst = np.maximum(0.0, 2.0 * np.sin(0.9 * i + s + 0.4) + 0.3) + dd_obs = np.mod(40.0 * i + 25.0 * s, 360.0) + dd_fcst = np.mod(dd_obs + 30.0 * np.sin(i + s) + 5.0, 360.0) + pt_obs = np.array([0, 1, 5, 7])[(i + s) % 4] * np.ones_like(t_obs) + pt_fcst = np.array([0, 1, 5, 7])[(i + 2 * s + 1) % 4] * np.ones_like(t_obs) + fcst = xr.Dataset( + { + "T_2M": (dims, t_fcst), + "TOT_PREC": (dims, tp_fcst), + "DD_10M": (dims, dd_fcst), + "PTYPE": (dims, pt_fcst), + "FCST_ONLY": (dims, t_fcst), + }, + coords=coords, + ) + obs = xr.Dataset( + { + "T_2M": (dims, t_obs), + "TOT_PREC": (dims, tp_obs), + "DD_10M": (dims, dd_obs), + "PTYPE": (dims, pt_obs), + "CAPE_MU": (dims, 100.0 + 10.0 * i + s * np.ones_like(t_obs)), + }, + coords=coords, + ) + return fcst, obs + + +STATION_THRESHOLDS = { + "T_2M": {"lt": [280.0], "gt": [278.5, 282.0]}, + "TOT_PREC": {"gt": [0.0, 1.0]}, + "DD_10M": {"gt": [180.0]}, + "PTYPE": {"eq": [0, 1, 5, 7]}, +} + + +def gridded_inputs() -> tuple[xr.Dataset, xr.Dataset]: + """INCA-like truth on a regular (y, x) grid, with gaps.""" + ny, nx, steps = 5, 6, np.array([0, 6]) + y = np.arange(ny)[None, :, None] + x = np.arange(nx)[None, None, :] + s = np.arange(steps.size)[:, None, None] + lat, lon = np.meshgrid( + np.linspace(46.0, 47.0, ny), np.linspace(5.0, 10.0, nx), indexing="ij" + ) + coords = { + "y": np.arange(ny), + "x": np.arange(nx), + "latitude": (("y", "x"), lat), + "longitude": (("y", "x"), lon), + "forecast_reference_time": _FRT, + "step": steps, + } + dims = ("step", "y", "x") + t_obs = 275.0 + np.sin(0.5 * y + 0.8 * x + s) + t_fcst = t_obs + 0.5 * np.cos(x - y + s) + t_fcst[1, 0, 0] = np.nan + t_obs[0, 4, 5] = np.nan + tp_obs = np.maximum(0.0, np.sin(y * x + s)) + tp_fcst = np.maximum(0.0, np.sin(y * x + s + 0.2)) + fcst = xr.Dataset( + {"T_2M": (dims, t_fcst), "TOT_PREC": (dims, tp_fcst)}, coords=coords + ) + obs = xr.Dataset({"T_2M": (dims, t_obs), "TOT_PREC": (dims, tp_obs)}, coords=coords) + return fcst, obs + + +def cases() -> dict: + """name -> verify() keyword arguments.""" + st_fcst, st_obs = station_inputs() + gr_fcst, gr_obs = gridded_inputs() + common = dict(fcst_label="fcst", obs_label="obs", regions=_REGIONS, num_workers=1) + return { + "stations": dict( + fcst=st_fcst.chunk({"step": 1}), + obs=st_obs, + threshold_dict=STATION_THRESHOLDS, + holdout_stations=["S00", "S03", "S07", "S11"], + **common, + ), + "stations_tolerant": dict( + fcst=st_fcst, + obs=st_obs, + threshold_dict=STATION_THRESHOLDS, + max_missing_fraction=0.5, + **common, + ), + "gridded": dict( + fcst=gr_fcst.chunk({"step": 1}), + obs=gr_obs, + threshold_dict={"TOT_PREC": {"gt": [0.1, 0.5]}}, + **common, + ), + } + + +@pytest.mark.parametrize("name", list(cases())) +def test_verify_output_unchanged(name): + result = verify(**cases()[name]) + expected = xr.load_dataset(DATA / f"verify_golden_{name}.nc") + + xr.testing.assert_identical(result, expected) + for var in expected.variables: + assert result[var].dims == expected[var].dims, var + if expected[var].dtype.kind in "fiub": + assert result[var].dtype == expected[var].dtype, var + # Bit for bit, NaN at the same places + np.testing.assert_array_equal(result[var].values, expected[var].values) + + +if __name__ == "__main__": + DATA.mkdir(exist_ok=True) + for name, kwargs in cases().items(): + verify(**kwargs).to_netcdf(DATA / f"verify_golden_{name}.nc") From 8701de6041ccf66041204c91d5c20b0d2d2f342b Mon Sep 17 00:00:00 2001 From: Francesco Zanetta Date: Thu, 5 Nov 2026 16:35:45 +0100 Subject: [PATCH 2/7] aggregate ensemble scores: spread, ssr, brier and its decomposition --- .../data/aggregation_golden_deterministic.nc | Bin 0 -> 25964 bytes tests/unit/test_verification_aggregation.py | 216 ++++++++++++++++++ workflow/scripts/verification_aggregation.py | 71 +++++- 3 files changed, 285 insertions(+), 2 deletions(-) create mode 100644 tests/unit/data/aggregation_golden_deterministic.nc create mode 100644 tests/unit/test_verification_aggregation.py diff --git a/tests/unit/data/aggregation_golden_deterministic.nc b/tests/unit/data/aggregation_golden_deterministic.nc new file mode 100644 index 0000000000000000000000000000000000000000..1e1db2713fb37299e6dc1fe5d2f6943b3a87bf59 GIT binary patch literal 25964 zcmeHP30PId`<@F|1i2`RD`H+1i9`{}EK?j2xB@B~kd+&vpkfFHDDGQm?z=7i&CI2z z(sC1=geraU|6V@ek?}5E2&DH6)^QP*^yr!(uwJzgp~%Q)Oi^7%DKIJNP`K!w!mw zRANAj?ImoIu!DpvGdw!Wkj`_P^Rh=f?q@vbb5dIa*T!6`sR3i%TH(+nI8+E~=^ryB zCW<(+Qu8)o(G?!~x0K6+wsOiJFghWh_ApgdOu)g#U|=K(9<2En=8Ak1uR2#KjH*4U zP?<}tn<&hQ$P*knog})xfG-HI$Wol{Zz;Y!sx@mP6K|)85K&9N`1pjVn7GI!W-le` z${LF4R-D#1F|rIS8yjJ8ew}^8!W2W}#o~rmYye??zMTVz$(75kS!*j68sgg}K#5;T zjUN*ilN2#1epmv-&R^$GvV?^N5I40W7|ev6V*&8nl()MXA4P?Q$exTlp zQFBcrmN9d%K8);-)VBDzReG9OgWjvZ^3j^tJuxYIsFLdg51G6e+*m4*CB4~K=7L41 zBxHcVf5-okhk9%jtx7blia7{=>AS)?R(P~>io2*uH?=}&W#gk9L3q(cU!1_lROVU> zPn8yk*@Q%^Qsyn$VjgG-4C+9z0O3JXTZjnz)j368Y7CjOVlJOXO>D%xmw}jAa&92~ z?^&eR}JK+^co>%*I(MoKkP&8Ab&%rN!D_v!I3zD{|B!w{vcYjM8*)T<= z?umqW&WZn>(^|rL)IVgav_olU;=v&k`|yNS{&vd~{@@GG(1= z=iFnfViZc!QmyhZ^ugx^;ML*GbJmqv#9yrX?Ag+hFCHy1v9aAFV~0hPdoj%y`6?kr z?a4A?u3&uGJ}xA)Jh>(;=0-?(q=ZLFcr?QkTLlm7=-0KgQet&}W#9PI% z%7kHj;Xnvq7ipbE3)uIL5)M=`*y^WZ^%<*(y+~C)TuRUD&z(BNlU`l|a~`;Z{Aslt zaXz?%JXx&(=Y>1S*YHMXj`PDE~;1?Q=9@PPWxm*lAL0!fbgE|lb`?;=T# z`Yw~?(EGk5hu-Cq9C}wsa_C(t$)R_xB!}L0k{o*1OLFMlAjzS3qa=sk9g-Y+Qzbd{ zrb%+>-6_eTcb6oG-tQzi^d6Mt(0fRdL+@cp4!uVtIrN^A%~YIm+2ilB1kKk{o*RfGz9oB*lkbvG^&rE9mVa#fM&DC2ZY=4gptCExb zM1~0ulc7j+eBg@e$k35I-|H`wP*nQGx9VCW!>Z(3pTAH-k@RP8w~u&+Xf>m&kUw;* zO%2Ym_u@I4Uf?;E{ia}PH#_a$tN7;X39r<5&i3p73k{m@%0B04`mVKm{70RodFJHX zvFp!My>P)wfx2~cx9`68{cTfnuh6|6YHpZ2`6|5@lX7XP(%e{DekV)Nc&BF?KhpJN39oa>{)P}e)F0O`<&3_M=$-Wg{Rwg$j_p0 zdfN@JvEvkV+0}r3PEn;gwNC4nUoU??Kw!U@@5A!`4b1mVwK$MrJfh$rrB0zZ*LZ2 z+Y!H?|MR9H?6aS4&*~CC>ihjP=LUCYX&=9l4ws(!O1GW#(vNE6*Dw8G8a+C=+VX-K z2WaNVr5Piao}t~R+ihK4{d?W=>*c?y(P*3HM~~3d9o(hUa7#4%T%;o11*#W*KVBfu z?f$>kE8RBJ$=NS#eqwSqUHlGrS@eTd++EV;M=yQzt(nVbc_RJTb~PbWAy$v zoAta*{S|jnxBPnfyKbIyw%w#`+VN#RLZz$ay+J!ZAG>n5`(EWE&}ChI^wK9!`D)v{4RUDf-IF>+kItaGy-e(LohEEH6qDX@!I`o3%Qy93-6(s8>_KD%aYp4!iD=<=hN z{zb**P^G&r%%NsFuCO2{kJ79K?Xni+>Xu(Gf4;YaagXnVYBAq0)ndM{s>OVNmhZtL zeJ1ESCy+56G_eD zq<{L2C9^^Vu_`vR*>X2G4CV~v`pS;qcL)#kXE$|NJZ>LR?jt;uK&As)j?VG@Nr`2{ zV^tBc*yRE4ies6&bm(tNUxwJ4uqdRH73D%9cqz4s}V5gT6DBtBEe;TtW4DWDk=%lTZ9o-X=AaWEWpbUQIq$w9(ax8^~5g+pLc= z_mDP%Z`6>a5Hs)5XIKgpIuN3CwfKbtfl>vpRa9dZLXtkmqE6l4Z zLWK=h0~8ISDt3sn8I%$m$$np+@M%*E}BF6T?j26C*koDp1|YBom->?PaF$x!Sn8r$Jcbj@a*y=lI4}G|Tz?yXI%^Pa-;<&}Cro}bj~@8i66jy@ z665(*ktqfvy0J~7azP5GVg@6+yIs+(bDYQ-j2Jg;#O-MujTwv>L2SegVVsQ_j2KId z{7x*6#tcUEZF|j4$S>3}%V0!TuXgu^y54FL3`UG_PDPO{zM&E9RGOVWKc}iS7%}kO z6^or+V~5UQ#2{82t6NupHLoTudDAdEHDYYLD;CRZjyfY7j2QBYV_7-`sd*JXkGwDv z-W7}Gg*{7y5d&m#tbp(^HLu@%kr#$SyYL9_E}k8Ne8YbC6Jh+ncE4;-fRLMcr(>?v zyPOM7FI=O|@Ulf@^w1eLT=2vHL~qV7H#Z)|DDbjbRP6=fJ;)_kC2HC4A(qp@d>$i4 zJ9y8RFrT)y;R_|ah+#rzcI!~lS)1zjB&of^1*i%zw8jf5G!-;U=&$OhpNQ<=3U`_0timmEs|7Kjb}g}fdH zvw@st%VJ_NY*|Q}?wj{{`Rrutd10A*+{qXT6%;B;8sai0(;DXc5?(Ih6%t;_Fd=Iu ze6vm%Kwvo=WpB4q3P=VE?;jjx*+YD(uB8uicqxh(8p8qWdWx&Pr=d^z>}<=Q5G8!v z*%%M^>TGM7T-({!NqD`4`IxSazmZ|KsJh%p`Mm-TDQ~(fViG94(l5FaD!%jK&cc=KO@UQS(U%=L@yuhD)E@gfcaq<;CrlBrjI|M&CtU4EtoTN z_L(bn$!CUnmpDj{dQxOmlwp~!Ms~2wwRVy!;WP>Fl<+Qwtxox8ieY7bIUB;Sv(j+~ zgT&+VpZuav8G*|4LbSb?A6B z*2ynwYpt_T!u(>kHvA6>-(fgrXJT*Wkv?-*xRlSyZ|c*$e2-U}3JLd`{Jvu$+Gai| z;X@KWEa4*z6Y_Kbom)PKDe9{f@IO6FeH+|$b@4EzH!K|pkK0ty?b{~JBC=AS8{L_G zNZO&Vz8U*P^-AqKR53dqUsWkhFI;eb{~H{;eyhD?&d&T)=lvZ`r>EaaCH`&O^uDz> z^@|S+2j2ShH{BCvz4VsFvp#Ly=tk=L*A|Zq*-|l0IT^3kG3}5krcy#(Pu=qC<-gOi za^mJt>$K!;KUN#_f7eqx&hlWNU#U;bsim%MYNzRi>+z$8-&ya}hJBi+dDU(C+Xi3z z_HQ-5G-2;hr}lp}eYKD5_VoOzcKaBo>GNjnwo3b$yDNpX;-nT8zx>1ZLCb&N@@@S( zJ@&P>`%+1tv%c3Ow$p?`IR!mjuG3{6w_UCrJ4^M#_3G!o_T6{Cx7fdsDf8xGM*M9WUf3j)Hs>Zucq^|tHBV%rzA5w2*c+_pv{1mlYccAi> z&=b1l*UP^{?10)Uj^@*;Yq~F6|70%pXtZL6X?-58eDGGv>bAL5FIE=4t=CsLPLD`ks!PmwoVI1|4>(=1A+_ zS#(vs?mfQRmrY;)BBjQz{7l{Q>*b#@zf$kuy?M0R?T!XN%T2oUa!h_y`+Pb#bZ-6b z>;{cqxE?<)o;%*VdM1~)-?J)s)SY~~+%IHK+?E{L>GB_QM)%XVAHDROcedPOHztqP z{q?;`{Z8Mcr_x-EhQ>GP*Y$Te%=gaIEx%s=>8V2pwpfr)_e}dPD!ocBZTN1dA53nU z^zBoVR~>J5gX)Fr@iR8JL5N%L4Ao8u>Ar01x3}x%(hWZ)9`v4~&yQaE^rj7_+-Z|X z<3{(mn0@{Ry}G2uQtQ-QI{3pbJtCgS)h)kX{-S>O|6e{A^FKXa%>Vg2;!d#^)h^~w zOE_J^XC!=bAz&M%(#u*H-yy^N^pSW>giu}XW z4h`Jw70w6hK*;UzlGwd_NgU#Nk#8Q7#FtO@Rv7f>?NKHdGJ2KvoQxj>aR*;(7yL)O zTFigE8+NtT_Tt!16;>@Od4+)GP!Me};?9n(a;lw}6a~%=M%>@Ai;GB%95OUEI?=PID2u^}%RC~>uyFqX zw2i^&q87NKl4Yp({xDVW*!xok#3aDp7^qnNKmST1ak3iXBMSYFTK0cz-;t1MbEj`3 z0 z)E2XO*CAK!)_CWrlC4y>@6O`yj?fmFzcWIM`MV;tn7<=Ji}|}D@_IC=rj!wD2~2;$ z$JZ9*m(OK(w|0VD<+z4@)gLYXV_ks1^F`a<&JwOGVg7CxZG7=gmm{Q)hcE#4+D%tZ zxfOdo9Ig}H*Duvu{KPPRyKSi(RJk!&s-rrrGvRZUywaiwcB(ObCQ3pG?vx6=DdTvQ zwttUlB;@=7f;$r%RjQWhqX_PP!pEZsijPMTU_6Qd<52_{k0QX>We3LLA23$Rz<3k^ z#&7|Q6*4fE#lTqZ0%MU2oO3?VqCArPZA@aI1CzjD5*SPZgGpd82@EEI!6Y!41O}79 zU=kQi0)t6lFbN#zy~3N@WX5;h60|mZOto1r9>|=DO$*h^)MjxL89r=Q=qtVl4h*J& z!89U>X=q1A}Q`Fbxc*fx$E|9^Zh$JTRCCc3o%FQp@;Em2tBg;qMzIFgnAw zfx$R17zYO9z+fC0j01ylU@#60#(}{&Fc=30!@yt|*mhRx(z4vx!RP~?fWZ?mcmf7b zz~BiOJOP6zVDJPCo`As@8z+e;@i~@sE zU@!^{MuEX7Fc<{}qrhMk7>oj=uLFZ&U|1sXg=@>dyw#@j(}(w6qG#TlWS1PnPHQF% zUQo#+n;tuJd-KK5FH^m6!PNavp>tOMF*W(!V|0gmx>KjmGw6n+uU@b9?s>Yv&i?rJ y NaN SSR, SPREAD is 0 + assert np.isnan(sel["T_2M.SSR"].isel(step=0)).all() + assert (sel["T_2M.SPREAD"].isel(step=0) == 0).all() + assert out.attrs["n_members"] == N + + +def test_ssr_nan_where_rmse_zero(): + ds = _ens_dataset() + ds["T_2M.MSE"] = ds["T_2M.MSE"] * 0 + ds["T_2M.BIAS"] = ds["T_2M.BIAS"] * 0 + out = aggregate_results(ds) + assert np.isnan(out["T_2M.SSR"]).all() + + +def test_ssr_n_inferred_from_rank(): + ds = _ens_dataset() + ds.attrs.pop("n_members") + out = aggregate_results(ds) + ref = aggregate_results(_ens_dataset()) + xr.testing.assert_allclose(out["T_2M.SSR"], ref["T_2M.SSR"]) + assert out.attrs["n_members"] == N + + +def test_counts_conversion_with_stratification(): + ds = _ens_dataset() + out = aggregate_results(ds) + # rank_hist: sum over inits of the per-init counts, in every stratum + total = ds["T_2M.rank_hist"].sum("forecast_reference_time") + np.testing.assert_allclose( + out.sel(season="all", init_hour=-999)["T_2M.rank_hist"], total + ) + djf = ( + ds["T_2M.rank_hist"] + .isel(forecast_reference_time=[0, 1]) + .sum("forecast_reference_time") + ) + np.testing.assert_allclose( + out.sel(season="DJF", init_hour=-999)["T_2M.rank_hist"], djf + ) + jja_12 = ds["T_2M.rank_hist"].isel(forecast_reference_time=2) + np.testing.assert_allclose( + out.sel(season="JJA", init_hour=12)["T_2M.rank_hist"], jja_12 + ) + # prob_table split per threshold, dims n_members_yes, obs + pt = out["TOT_PREC.prob_table_gt_0p1"] + assert set(pt.dims) >= {"n_members_yes", "obs"} + exp = ( + ds["TOT_PREC.prob_table"] + .isel({"TOT_PREC.threshold": 0}) + .sum("forecast_reference_time") + ) + np.testing.assert_allclose(pt.sel(season="all", init_hour=-999), exp) + assert "TOT_PREC.prob_table" not in out + + +def test_brier_and_murphy_decomposition(): + ds = _ens_dataset() + out = aggregate_results(ds).sel(season="all", init_hour=-999) + counts = ds["TOT_PREC.prob_table"].sum("forecast_reference_time") + for i, label in enumerate(["gt_0p1", "gt_1p0"]): + c = ( + counts.isel({"TOT_PREC.threshold": i}) + .transpose(..., "n_members_yes", "obs") + .values + ) + t = c.sum(axis=(-1, -2)) + f = np.arange(N + 1) / N + n1, n0 = c[..., 1], c[..., 0] + nk = n1 + n0 + bs = ((n1 * (f - 1) ** 2 + n0 * f**2).sum(-1)) / t + obar = n1.sum(-1) / t + ok = n1 / nk + rel = (nk * (f - ok) ** 2).sum(-1) / t + res = (nk * (ok - obar[..., None]) ** 2).sum(-1) / t + unc = obar * (1 - obar) + got = { + k: out[f"TOT_PREC.{k}_{label}"].transpose("source", "region", "step").values + for k in ["BS", "BS_REL", "BS_RES", "BS_UNC"] + } + for k, ref in zip(["BS", "BS_REL", "BS_RES", "BS_UNC"], [bs, rel, res, unc]): + np.testing.assert_allclose(got[k], ref[0] if ref.ndim == 4 else ref) + np.testing.assert_allclose( + got["BS"], got["BS_REL"] - got["BS_RES"] + got["BS_UNC"] + ) + + +def test_brier_hand_example(): + # one init pair, N=2; bins k=0,1,2 with (obs0, obs1) counts + ds = _ens_dataset().isel(forecast_reference_time=[0], region=[0], step=[0]) + tab = np.zeros((1, 1, 1, 1, 2, 3, 2)) + tab[..., 0, :, :] = [[6, 0], [1, 1], [0, 2]] # k=0: 6 no,0 yes; k=1: 1,1; k=2: 0,2 + tab[..., 1, :, :] = tab[..., 0, :, :] + ds = ds.isel(rank=[0, 1, 2], n_members_yes=[0, 1, 2]) + ds["TOT_PREC.prob_table"] = (ds["TOT_PREC.prob_table"].dims, tab) + ds.attrs["n_members"] = 2 + out = ( + aggregate_results(ds) + .sel(season="all", init_hour=-999) + .isel(source=0, region=0, step=0) + ) + # T=10, f=(0,.5,1): BS = (0.25*(1+1) + 0)/10 = 0.05 ; obar = 3/10 + assert np.isclose(out["TOT_PREC.BS_gt_0p1"], 0.05) + assert np.isclose(out["TOT_PREC.BS_UNC_gt_0p1"], 0.21) + # REL = (2*(0.5-0.5)^2 + 0)/10 = 0 ; RES = (6*.09 + 2*.04 + 2*.49)/10 = 0.16 + assert np.isclose(out["TOT_PREC.BS_REL_gt_0p1"], 0.0) + assert np.isclose(out["TOT_PREC.BS_RES_gt_0p1"], 0.16) diff --git a/workflow/scripts/verification_aggregation.py b/workflow/scripts/verification_aggregation.py index 60c7c492..8d7acb5e 100644 --- a/workflow/scripts/verification_aggregation.py +++ b/workflow/scripts/verification_aggregation.py @@ -34,6 +34,40 @@ def get_season(time): return "SON" +def _brier_from_prob_table(table: xr.DataArray) -> xr.Dataset: + """Brier score and Murphy decomposition from a prob_table of counts. + + ``table`` has dims ``n_members_yes`` (0..N) and ``obs`` (0, 1) besides the others. + The forecast probability of bin k is k/N. BS = REL - RES + UNC holds exactly as + the forecast probability takes one value per bin. + """ + n_bins = table.sizes["n_members_yes"] + n_members = n_bins - 1 + f = xr.DataArray( + np.arange(n_bins) / n_members, + dims="n_members_yes", + coords={"n_members_yes": table["n_members_yes"]}, + ) + n_obs1 = table.isel(obs=1, drop=True) + n_obs0 = table.isel(obs=0, drop=True) + n_k = n_obs0 + n_obs1 + total = n_k.sum("n_members_yes", skipna=False) + total = total.where(total > 0) + o_k = (n_obs1 / n_k.where(n_k > 0)).fillna(0.0) + o_bar = n_obs1.sum("n_members_yes", skipna=False) / total + bs = (n_obs1 * (f - 1) ** 2 + n_obs0 * f**2).sum("n_members_yes", skipna=False) + rel = (n_k * (f - o_k) ** 2).sum("n_members_yes", skipna=False) + res = (n_k * (o_k - o_bar) ** 2).sum("n_members_yes", skipna=False) + return xr.Dataset( + { + "BS": bs / total, + "BS_REL": rel / total, + "BS_RES": res / total, + "BS_UNC": o_bar * (1 - o_bar), + } + ) + + def aggregate_results(ds: xr.Dataset) -> xr.Dataset: """Compute mean metric values aggregated by season""" LOG.info("Aggregation results") @@ -74,8 +108,23 @@ def aggregate_results(ds: xr.Dataset) -> xr.Dataset: ds_mean.append(ds_result) out = xr.merge(ds_mean, compat="no_conflicts", join="outer") - for var in out.data_vars: - if "contingency_table" in var: + n_members_attr = ds.attrs.get("n_members") + for var in list(out.data_vars): + if "rank_hist" in var: + # ensemble rank histogram: back to per-stratum counts + out[var] = out[var] * out["n_samples"] + if n_members_attr is None: + n_members_attr = out.sizes["rank"] - 1 + elif "prob_table" in var: + # forecast-probability/observation counts: back to per-stratum counts, + # then Brier score and its decomposition per threshold + out[var] = out[var] * out["n_samples"] + brier = _brier_from_prob_table(out[var]) + for name, da in brier.data_vars.items(): + out[var.replace("prob_table", name)] = da + if n_members_attr is None: + n_members_attr = out.sizes["n_members_yes"] - 1 + elif "contingency_table" in var: # convert contingency table means back to per-stratum counts out[var] = out[var] * out["n_samples"] contingency_manager = scores.categorical.BasicContingencyManager( @@ -128,6 +177,24 @@ def aggregate_results(ds: xr.Dataset) -> xr.Dataset: out[var] = np.sqrt(out[var]) out = out.rename(var_transform) + # Ensemble spread and spread-skill ratio (after the sqrt/rename above: needs RMSE) + for var in [v for v in out.data_vars if v.endswith(".ENS_VAR")]: + prefix = var[: -len("ENS_VAR")] + spread = np.sqrt(out[var]) + out[f"{prefix}SPREAD"] = spread + rmse_var = f"{prefix}RMSE" + if rmse_var in out.data_vars and n_members_attr is not None: + n = float(n_members_attr) + rmse = out[rmse_var] + ssr = spread / rmse.where(rmse != 0) * np.sqrt((n + 1) / n) + out[f"{prefix}SSR"] = ssr.where(spread != 0) + out = out.drop_vars(var) + if n_members_attr is not None and any( + v.endswith(".SPREAD") or "rank_hist" in v or "prob_table" in v + for v in out.data_vars + ): + out.attrs["n_members"] = int(n_members_attr) + LOG.info("Computed aggregation in %.2f seconds", time.time() - start) LOG.info("Aggregated results: \n %s", out) From d1e6ed0289c3c8328204fcc3c4486e2ac74730ea Mon Sep 17 00:00:00 2001 From: Francesco Zanetta Date: Thu, 5 Nov 2026 16:35:48 +0100 Subject: [PATCH 3/7] show ensemble scores in dashboard, scorecards and metric plots --- .../report/dashboard/template.html.jinja2 | 1 + tests/unit/test_ensemble_consumers.py | 160 ++++++++++++++++++ workflow/scripts/ensemble_scores.py | 83 +++++++++ .../scripts/report_experiment_dashboard.py | 10 +- workflow/scripts/report_scorecard.py | 18 +- workflow/scripts/verification_plot_metrics.py | 100 ++++++++++- 6 files changed, 365 insertions(+), 7 deletions(-) create mode 100644 tests/unit/test_ensemble_consumers.py create mode 100644 workflow/scripts/ensemble_scores.py diff --git a/resources/report/dashboard/template.html.jinja2 b/resources/report/dashboard/template.html.jinja2 index 0da748f0..0415f5f4 100644 --- a/resources/report/dashboard/template.html.jinja2 +++ b/resources/report/dashboard/template.html.jinja2 @@ -191,6 +191,7 @@ + {% if crps_note %}

{{ crps_note | e }}.

{% endif %}
diff --git a/tests/unit/test_ensemble_consumers.py b/tests/unit/test_ensemble_consumers.py new file mode 100644 index 00000000..613646f1 --- /dev/null +++ b/tests/unit/test_ensemble_consumers.py @@ -0,0 +1,160 @@ +"""Consumers of aggregated verification (dashboard, scorecard, plots) and ensemble scores.""" + +import sys +from argparse import Namespace +from pathlib import Path + +import numpy as np +import pytest +import xarray as xr + +ROOT = Path(__file__).parents[2] +sys.path.insert(0, str(ROOT / "workflow" / "scripts")) +import report_experiment_dashboard as dash # noqa: E402 +import report_scorecard as sc # noqa: E402 +import verification_plot_metrics as vpm # noqa: E402 +from ensemble_scores import ( # noqa: E402 + add_deterministic_crps, + drop_non_score_vars, + ensemble_sources, + fill_crps_from_mae, +) + +N = 3 +STEPS = np.array([6, 12], dtype="timedelta64[h]").astype("timedelta64[ns]") + + +def _coords(source): + return { + "source": [source], + "region": ["icon"], + "season": ["all"], + "init_hour": [-999], + "station_group": ["all"], + "step": STEPS, + } + + +def _var(val, extra=None): + dims = ["source", "region", "season", "init_hour", "station_group", "step"] + shape = [1, 1, 1, 1, 1, 2] + arr = np.full(shape, val, dtype=float) + if extra: + for d, n in extra.items(): + dims.append(d) + shape.append(n) + arr = np.random.default_rng(0).uniform(1, 5, shape) + return dims, arr + + +def make_file(source, ensemble=False): + ds = xr.Dataset(coords=_coords(source)) + for name, val in [("T_2M.RMSE", 1.0), ("T_2M.MAE", 0.8), ("T_2M.BIAS", 0.1)]: + ds[name] = _var(val) + ds["n_samples"] = (("season", "init_hour"), [[4]]) + if ensemble: + for name, val in [ + ("T_2M.CRPS", 0.5), + ("T_2M.SPREAD", 0.9), + ("T_2M.SSR", 1.1), + ("T_2M.MEMBER_RMSE", 1.2), + ("T_2M.BS_gt_0p1", 0.2), + ]: + ds[name] = _var(val) + ds["T_2M.rank_hist"] = _var(0, {"rank": N + 1}) + ds["T_2M.prob_table_gt_0p1"] = _var(0, {"n_members_yes": N + 1, "obs": 2}) + return ds + + +def test_drop_non_score_vars(): + ds = drop_non_score_vars(make_file("ens", ensemble=True)) + assert "T_2M.rank_hist" not in ds and "T_2M.prob_table_gt_0p1" not in ds + assert {"T_2M.CRPS", "T_2M.BS_gt_0p1", "T_2M.SSR"} <= set(ds.data_vars) + + +def test_deterministic_crps_is_mae_for_display_only(): + ens, det = make_file("ens", True), make_file("det") + ds = vpm._concat_sources([ens, det]) + assert ensemble_sources(ds) == ["ens"] + out, filled = add_deterministic_crps(ds) + assert filled + assert out["T_2M.CRPS"].sel(source="det").isel(step=0).item() == 0.8 + assert out["T_2M.CRPS"].sel(source="ens").isel(step=0).item() == 0.5 + assert "T_2M.CRPS" not in det # input file untouched + + +def test_deterministic_only_is_unchanged(): + ds = vpm._concat_sources([make_file("a"), make_file("b")]) + out, filled = add_deterministic_crps(ds) + assert not filled and out is ds + single, filled = fill_crps_from_mae(make_file("a")) + assert filled and single["T_2M.CRPS"].equals(single["T_2M.MAE"]) + ens, filled = fill_crps_from_mae(make_file("e", True)) + assert not filled and float(ens["T_2M.CRPS"].max()) == 0.5 + + +def _run_dashboard(tmp_path, files): + paths = [] + for i, f in enumerate(files): + p = tmp_path / f"f{i}.nc" + f.to_netcdf(p) + paths.append(p) + out = tmp_path / "out" + dash.main( + Namespace( + verif_files=paths, + template=ROOT / "resources/report/dashboard/template.html.jinja2", + script=ROOT / "resources/report/dashboard/script.js", + header_text="", + stratification=[], + station_holdout=False, + configfile=tmp_path / "none.yaml", + output=out, + label_map={}, + ) + ) + return (out / "dashboard.html").read_text() + + +def test_dashboard(tmp_path): + html = _run_dashboard(tmp_path, [make_file("ens", True), make_file("det")]) + for m in ["CRPS", "SPREAD", "SSR", "MEMBER_RMSE", "BS > 0.1"]: + assert f'data-metric="{m}"' in html + assert "rank_hist" not in html and "prob_table" not in html + assert "CRPS = MAE for deterministic forecasts" in html + + +def test_dashboard_deterministic_only_has_no_ensemble_content(tmp_path): + html = _run_dashboard(tmp_path, [make_file("a"), make_file("b")]) + assert "CRPS" not in html + assert "deterministic forecasts" not in html + + +def test_scorecard_crps(tmp_path): + ens, det = make_file("ens", True), make_file("det") + pe, pd_ = tmp_path / "e.nc", tmp_path / "d.nc" + ens.to_netcdf(pe) + det.to_netcdf(pd_) + cfg = { + "model": {"path": str(pe), "source": "ens", "label": "ens"}, + "baseline": {"path": str(pd_), "source": "det", "label": "det"}, + "stratification": "region", + "lead_times": "6/12/6", + "all_metrics": sc.ALL_METRICS if hasattr(sc, "ALL_METRICS") else ["RMSE"], + "variables": {"T_2M": ["CRPS"]}, + } + diff = sc._load_relative_diff(cfg) + assert cfg["crps_filled"] + assert "T_2M.rank_hist" not in diff + # CRPS 0.5 vs MAE 0.8 -> -37.5 % + assert float(diff["T_2M.CRPS"].isel(step=0)) == pytest.approx(-37.5) + assert sc._is_model_better(-37.5, "CRPS", ["CRPS"], []) is True + + +def test_plot_ensemble_diagnostics(tmp_path): + ds = vpm._concat_sources([make_file("ens", True), make_file("det")]) + ds = ds.sel(station_group="all", drop=True) + args = Namespace(output_dir=tmp_path, label_map={}) + vpm.plot_ensemble_diagnostics(ds, args) + names = {p.name for p in tmp_path.iterdir()} + assert names == {"rank_hist_T_2M.png", "reliability_T_2M_gt_0p1.png"} diff --git a/workflow/scripts/ensemble_scores.py b/workflow/scripts/ensemble_scores.py new file mode 100644 index 00000000..362461f2 --- /dev/null +++ b/workflow/scripts/ensemble_scores.py @@ -0,0 +1,83 @@ +"""Helpers for consumers of aggregated verification files (ensemble scores). + +Shared by the dashboard, the scorecards and the metric plots. Importing this module +has no side effects (no plotting style, no logging setup). +""" + +import numpy as np +import xarray as xr + +# Variables carrying one of these dims are counts/tables kept for plots, not scores. +NON_SCORE_DIMS = frozenset({"rank", "n_members_yes", "obs"}) + +CRPS_NOTE = "CRPS = MAE for deterministic forecasts" + + +def is_score_var(da: xr.DataArray) -> bool: + """True if the variable is a score that can be tabulated or drawn as a curve.""" + return not (set(da.dims) & NON_SCORE_DIMS) + + +def drop_non_score_vars(ds: xr.Dataset) -> xr.Dataset: + """Drop the rank histogram and probability tables (not scores).""" + drop = [v for v in ds.data_vars if not is_score_var(ds[v])] + return ds.drop_vars(drop) if drop else ds + + +def _spread_vars(ds: xr.Dataset) -> list[str]: + return [v for v in ds.data_vars if v.endswith(".SPREAD")] + + +def ensemble_sources(ds: xr.Dataset) -> list: + """Sources of *ds* that carry ensemble scores (any finite SPREAD).""" + found: set = set() + for v in _spread_vars(ds): + if "source" not in ds[v].dims: + continue + other = [d for d in ds[v].dims if d != "source"] + has = ds[v].notnull().any(other) if other else ds[v].notnull() + found.update(ds["source"].values[has.values].tolist()) + return [s for s in ds["source"].values.tolist() if s in found] + + +def add_deterministic_crps(ds: xr.Dataset) -> tuple[xr.Dataset, bool]: + """Multi-source dataset: give deterministic sources ``

.CRPS =

.MAE``. + + Only when the dataset contains an ensemble source; otherwise returned unchanged. + The result is for display only. Returns ``(dataset, filled)``. + """ + ens = ensemble_sources(ds) + if not ens: + return ds, False + is_det = xr.DataArray( + ~np.isin(ds["source"].values, ens), + dims="source", + coords={"source": ds["source"]}, + ) + ds = ds.copy() + filled = False + for mae in [v for v in ds.data_vars if v.endswith(".MAE")]: + if "source" not in ds[mae].dims: + continue + name = mae[: -len("MAE")] + "CRPS" + base = ds[name] if name in ds else xr.full_like(ds[mae], np.nan) + out = xr.where(is_det, ds[mae], base) + ds[name] = out.transpose(*ds[mae].dims) + ds[name].attrs = dict(ds[mae].attrs) + filled = True + return ds, filled + + +def fill_crps_from_mae(ds: xr.Dataset) -> tuple[xr.Dataset, bool]: + """Single-participant dataset: if it is deterministic, set ``

.CRPS =

.MAE``.""" + if any(bool(ds[v].notnull().any()) for v in _spread_vars(ds)): + return ds, False + ds = ds.copy() + filled = False + for mae in [v for v in ds.data_vars if v.endswith(".MAE")]: + name = mae[: -len("MAE")] + "CRPS" + if name in ds and bool(ds[name].notnull().any()): + continue + ds[name] = ds[mae] + filled = True + return ds, filled diff --git a/workflow/scripts/report_experiment_dashboard.py b/workflow/scripts/report_experiment_dashboard.py index 75a13efb..c45d13a1 100644 --- a/workflow/scripts/report_experiment_dashboard.py +++ b/workflow/scripts/report_experiment_dashboard.py @@ -14,6 +14,7 @@ _ensure_unique_lead_time, _select_best_sources, ) +from ensemble_scores import CRPS_NOTE, add_deterministic_crps, is_score_var from verification import decode_metric LOG = logging.getLogger(__name__) @@ -69,10 +70,16 @@ def main(args): dfs = [_ensure_unique_lead_time(d) for d in dfs] dfs = _select_best_sources(dfs) ds = _concat_sources(dfs) + ds, crps_filled = add_deterministic_crps(ds) LOG.info("Loaded verification netcdf: \n%s", ds) # extract only non-spatial variables to pd.DataFrame - nonspatial_vars = [d for d in ds.data_vars if "spatial" not in d and "." in d] + # (rank histograms and probability tables are counts for plots, not scores) + nonspatial_vars = [ + d + for d in ds.data_vars + if "spatial" not in d and "." in d and is_score_var(ds[d]) + ] df = ds[nonspatial_vars].to_array("stack").to_dataframe(name="value").reset_index() df[["param", "metric"]] = df["stack"].str.split(".", n=1, expand=True) df["metric"] = df.metric.apply(decode_metric) @@ -175,6 +182,7 @@ def _sanitize(v): station_holdout=station_holdout, station_groups=station_groups, header_text=args.header_text, + crps_note=CRPS_NOTE if crps_filled else "", configfile_content=open(args.configfile, "r").read() if args.configfile.is_file() else "", diff --git a/workflow/scripts/report_scorecard.py b/workflow/scripts/report_scorecard.py index 5d8b4662..40e0ba04 100644 --- a/workflow/scripts/report_scorecard.py +++ b/workflow/scripts/report_scorecard.py @@ -6,7 +6,7 @@ diff = (model − baseline) / |baseline| × 100 [%] Positive diff is better for higher-is-better metrics (R², ETS, POD) and worse -for lower-is-better metrics (RMSE, MAE, STDE, FAR). +for lower-is-better metrics (RMSE, MAE, STDE, FAR, CRPS). """ import logging @@ -20,6 +20,7 @@ import matplotlib.pyplot as plt from matplotlib.transforms import ScaledTranslation +from ensemble_scores import CRPS_NOTE, drop_non_score_vars, fill_crps_from_mae from verification import decode_metric LOG = logging.getLogger(__name__) @@ -202,7 +203,7 @@ def _build_config(args) -> dict: # All recognised metrics — every entry must also appear in metric_directions. "all_metrics": ["RMSE", "MAE", "STDE", "R2", "ETS", "POD", "FAR"], "metric_directions": { - "lower_is_better": ["RMSE", "MAE", "STDE", "FAR"], + "lower_is_better": ["RMSE", "MAE", "STDE", "FAR", "CRPS"], "higher_is_better": ["R2", "ETS", "POD"], }, # None value → use all_metrics for that variable; list → restrict to that subset. @@ -259,10 +260,14 @@ def _load_relative_diff(cfg: dict) -> xr.Dataset: f"Fix: delete '{fewer}' and rerun the pipeline." ) - model_ds = model_ds.sel(**sel_coords).squeeze(drop=True) - baseline_ds = baseline_ds.sel(**{**sel_coords, "source": baseline_source}).squeeze( - drop=True + model_ds = drop_non_score_vars(model_ds.sel(**sel_coords).squeeze(drop=True)) + baseline_ds = drop_non_score_vars( + baseline_ds.sel(**{**sel_coords, "source": baseline_source}).squeeze(drop=True) ) + # Deterministic participants have no CRPS: show their MAE (display only). + model_ds, model_filled = fill_crps_from_mae(model_ds) + baseline_ds, baseline_filled = fill_crps_from_mae(baseline_ds) + cfg["crps_filled"] = model_filled or baseline_filled common_vars = [v for v in model_ds.data_vars if v in baseline_ds.data_vars] common_leads = sorted(set(model_ds.step.values) & set(baseline_ds.step.values)) @@ -793,6 +798,9 @@ def _render_scorecard(diff: xr.Dataset, cfg: dict, outfn: Path): cfg, ) + if cfg.get("crps_filled") and any(m == "CRPS" for _, m in rows): + fig.text(0.01, 0.0, CRPS_NOTE, fontsize=small_fs, va="top", ha="left") + outfn.parent.mkdir(parents=True, exist_ok=True) plt.savefig(outfn, bbox_inches="tight") LOG.info("Saved %s", outfn) diff --git a/workflow/scripts/verification_plot_metrics.py b/workflow/scripts/verification_plot_metrics.py index dd0b6e54..b051a6ad 100644 --- a/workflow/scripts/verification_plot_metrics.py +++ b/workflow/scripts/verification_plot_metrics.py @@ -7,6 +7,12 @@ import matplotlib.pyplot as plt import pandas as pd import xarray as xr +from ensemble_scores import ( + add_deterministic_crps, + ensemble_sources, + is_score_var, + CRPS_NOTE, +) from verification import decode_metric LOG = logging.getLogger(__name__) @@ -111,9 +117,13 @@ def main(args: Namespace) -> None: dfs = _select_best_sources(dfs) # 3) Concatenate by source; outer join to keep the union of lead_times ds = _concat_sources(dfs) + ds, crps_filled = add_deterministic_crps(ds) # extract only non-spatial variables to pd.DataFrame - nonspatial_vars = [d for d in ds.data_vars if "spatial" not in d] + # (rank histograms and probability tables are not scores: drawn separately) + nonspatial_vars = [ + d for d in ds.data_vars if "spatial" not in d and is_score_var(ds[d]) + ] all_df = ( ds[nonspatial_vars].to_array("stack").to_dataframe(name="value").reset_index() ) @@ -156,6 +166,10 @@ def _subset_df(df): title = f"{metric} - {param} - {region}" title += f"- {season} - {init_hour}" if args.stratify else "" + if crps_filled and metric == "CRPS": + ax.annotate( + CRPS_NOTE, (0.0, -0.3), xycoords="axes fraction", fontsize=9, va="top" + ) for source, df in sub_df.groupby("source"): display_label = args.label_map.get(source, source) df.plot( @@ -181,6 +195,90 @@ def _subset_df(df): plt.savefig(args.output_dir / fn, bbox_inches="tight") plt.close(fig) + plot_ensemble_diagnostics(ds, args) + + +def _counts_for_plot(da: xr.DataArray, source, region, keep: list[str]) -> xr.DataArray: + """Counts of one source at the default stratum, summed over everything but *keep*.""" + da = da.sel(source=source) + for dim, val in (("region", region), ("season", "all"), ("init_hour", -999)): + if dim in da.dims: + da = da.sel({dim: val}, drop=True) + return da.sum([d for d in da.dims if d not in keep], skipna=True) + + +def plot_ensemble_diagnostics(ds: xr.Dataset, args: Namespace) -> None: + """Rank histograms and reliability diagrams, for ensemble participants only. + + Counts are summed over lead times, at the default region, all seasons and init hours. + """ + ens_sources = ensemble_sources(ds) + if not ens_sources: + return + region = ds["region"].values[0] + args.output_dir.mkdir(parents=True, exist_ok=True) + label = lambda s: args.label_map.get(s, s) # noqa: E731 + + for name in [v for v in ds.data_vars if v.endswith(".rank_hist")]: + param = name.rsplit(".", 1)[0] + fig, ax = plt.subplots(figsize=(10, 6)) + n_plotted = 0 + for source in ens_sources: + counts = _counts_for_plot(ds[name], source, region, ["rank"]) + total = float(counts.sum()) + if not total > 0: + continue + ax.plot( + counts["rank"].values, + counts.values / total, + marker="o", + label=label(source), + ) + n_plotted += 1 + if n_plotted: + n_ranks = ds.sizes["rank"] + ax.axhline(1 / n_ranks, color="grey", ls="--", lw=1) + ax.set( + title=f"Rank histogram - {param} - {region}", + xlabel="Rank of observation among members", + ylabel="Relative frequency", + ) + ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.15), ncol=n_plotted) + fig.savefig(args.output_dir / f"rank_hist_{param}.png", bbox_inches="tight") + plt.close(fig) + + for name in [v for v in ds.data_vars if ".prob_table_" in v]: + param, rest = name.split(".prob_table_", 1) + fig, (ax, axc) = plt.subplots( + 2, 1, figsize=(10, 8), sharex=True, height_ratios=[3, 1] + ) + n_plotted = 0 + for source in ens_sources: + tab = _counts_for_plot(ds[name], source, region, ["n_members_yes", "obs"]) + n_yes = tab["n_members_yes"].values + n_mem = int(n_yes.max()) + n_bin = tab.sum("obs") + freq = tab.sel(obs=1) / n_bin.where(n_bin > 0) + if not bool(n_bin.sum() > 0): + continue + prob = n_yes / n_mem + ax.plot(prob, freq.values, marker="o", label=label(source)) + axc.plot(prob, n_bin.values, marker="o") + n_plotted += 1 + if n_plotted: + ax.plot([0, 1], [0, 1], color="grey", ls="--", lw=1) + ax.set( + title=f"Reliability - {param} {decode_metric(rest)} - {region}", + ylabel="Observed frequency", + ) + axc.set(xlabel="Forecast probability (k/N)", ylabel="Count", yscale="log") + ax.legend(loc="upper left") + fig.savefig( + args.output_dir / f"reliability_{param}_{rest}.png", + bbox_inches="tight", + ) + plt.close(fig) + if __name__ == "__main__": parser = ArgumentParser(description="Verify results from KENDA-1 data.") From 1541c50d3df44cb1be6c612df14300a156a0c67f Mon Sep 17 00:00:00 2001 From: Francesco Zanetta Date: Thu, 5 Nov 2026 16:35:50 +0100 Subject: [PATCH 4/7] verify ensemble runs on all members, ensemble mean as deterministic forecast --- README.md | 22 +- tests/unit/test_ensemble_wiring.py | 288 ++++++++++++++++++++ workflow/rules/common.smk | 45 ++- workflow/rules/verification.smk | 24 +- workflow/scripts/verification_metrics.py | 182 ++++++++++++- workflow/scripts/verification_valid_mask.py | 13 + 6 files changed, 551 insertions(+), 23 deletions(-) create mode 100644 tests/unit/test_ensemble_wiring.py diff --git a/README.md b/README.md index 3838cc92..ed0d199c 100644 --- a/README.md +++ b/README.md @@ -382,11 +382,31 @@ runs: - A member that finished keeps a `.ok` file under `logs/inference_execute/{run_id}-{init_time}/`, so a rerun after a failure only runs the missing members. -- Verification, plots and MEC use member 001 only, with a warning. +- Verification scores the ensemble (see below); plots, SAL, score maps and MEC use + member 001 only, with a warning. - Not allowed on the inline `forecaster` of a `temporal_downscaler`. - Changing `members` or `seed` gives the run a new `run_id`; the inference environment is reused. +#### Ensemble verification + +An ensemble run is verified from all its members, in the usual `verif_*.nc` files of +the run: + +- Its deterministic scores (`BIAS`, `RMSE`, `MAE`, contingency scores, ...) are those + of the ensemble mean. There is no separate score for a single member: to verify one, + run it as a deterministic forecaster. +- Ensemble scores are added next to them: the fair `CRPS`, the ensemble spread + `SPREAD` and spread-skill ratio `SSR`, the mean score of the single members + `MEMBER_RMSE`, Brier scores per threshold, and the counts of the rank histogram and + reliability diagram. +- The forecast statistics (`mean`, `var`, `min`, `max`) are averaged over the members. +- A point where any member is missing counts as missing for the ensemble (and, with + `score_on_common_points`, is left out for every participant). +- Each member is loaded like a deterministic run (derived and period parameters, + lapse-rate correction with the orography of member 001). Memory and runtime of the + verification jobs grow with the number of members. + ## Installation diff --git a/tests/unit/test_ensemble_wiring.py b/tests/unit/test_ensemble_wiring.py new file mode 100644 index 00000000..1ade73d0 --- /dev/null +++ b/tests/unit/test_ensemble_wiring.py @@ -0,0 +1,288 @@ +"""Wiring of ensemble verification: the common.smk helpers of the verification rules, +member loading and scoring in verification_metrics.py, and the valid mask of an +ensemble in verification_valid_mask.py.""" + +import ast +import sys +from argparse import Namespace +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest +import xarray as xr + +ROOT = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(ROOT / "workflow" / "scripts")) + +import verification_metrics # noqa: E402 +import verification_valid_mask # noqa: E402 +from verification import verify # noqa: E402 + +COMMON = ROOT / "workflow/rules/common.smk" +FUNCTIONS = { + "run_grib_dir", + "ensemble_members", + "verification_forecast_arg", + "ensemble_scaled", +} + + +def _common_namespace(run_configs, out_root): + tree = ast.parse(COMMON.read_text()) + body = [ + node + for node in tree.body + if isinstance(node, ast.FunctionDef) and node.name in FUNCTIONS + ] + namespace = {"Path": Path, "RUN_CONFIGS": run_configs, "OUT_ROOT": out_root} + exec( + compile(ast.Module(body=body, type_ignores=[]), str(COMMON), "exec"), namespace + ) + return namespace + + +RUN_CONFIGS = { + "forecaster-aaaa-bbbb/cccc": {"steps": "0/24/6"}, + "forecaster-aaaa-bbbb/dddd": { + "steps": "0/24/6", + "ensemble": {"members": 4, "seed": 0}, + }, +} +DET, ENS = RUN_CONFIGS + + +def test_verification_forecast_arg_deterministic_unchanged(tmp_path): + """Deterministic runs: the same Path as before the ensemble wiring. Snakemake + does not record Path params, so the param never triggers a rerun.""" + ns = _common_namespace(RUN_CONFIGS, tmp_path) + arg = ns["verification_forecast_arg"](DET, "202406150000") + assert isinstance(arg, Path) + assert arg == (tmp_path / f"data/runs/{DET}/202406150000/grib").resolve() + assert arg == ns["run_grib_dir"](DET, "202406150000").resolve() + + +def test_verification_forecast_arg_ensemble(tmp_path): + ns = _common_namespace(RUN_CONFIGS, tmp_path) + grib = (tmp_path / f"data/runs/{ENS}/202406150000/grib").resolve() + assert ns["verification_forecast_arg"](ENS, "202406150000") == ( + f"{grib} --members 4" + ) + # plots, SAL, score maps and MEC keep member 001 + assert ns["run_grib_dir"](ENS, "202406150000") == ( + tmp_path / f"data/runs/{ENS}/202406150000/grib/001" + ) + + +def test_ensemble_scaled_resources(tmp_path): + ns = _common_namespace(RUN_CONFIGS, tmp_path) + mem = ns["ensemble_scaled"](80_000, 6_000) + runtime = ns["ensemble_scaled"](60, 6, unit="m") + assert mem(Namespace(run_id=DET)) == 80_000 + assert runtime(Namespace(run_id=DET)) == "60m" + assert mem(Namespace(run_id=ENS)) == 104_000 + assert runtime(Namespace(run_id=ENS)) == "84m" + + +# --- member loading ----------------------------------------------------------- + + +def _member_dataset(member, n_values=3): + values = np.full((1, 2, n_values), float(member)) + return xr.Dataset( + {"T_2M": (("forecast_reference_time", "step", "values"), values)}, + coords={ + "forecast_reference_time": [np.datetime64("2024-06-15T00:00", "ns")], + "step": pd.to_timedelta([0, 6], unit="h"), + "latitude": ("values", np.linspace(46.0, 47.0, n_values)), + "longitude": ("values", np.linspace(6.0, 9.0, n_values)), + "valid_time": ( + ("forecast_reference_time", "step"), + np.array( + [["2024-06-15T00:00", "2024-06-15T06:00"]], dtype="datetime64[ns]" + ), + ), + # differs per member on purpose: member 001's must be kept + "elevation": ("values", np.full(n_values, 100.0 * member)), + }, + ) + + +@pytest.fixture +def member_dirs(tmp_path, monkeypatch): + grib = tmp_path / "grib" + for m in (1, 2, 3): + (grib / f"{m:03d}").mkdir(parents=True) + (grib / f"{m:03d}" / "x_0.grib").touch() + calls = [] + + def fake_load(root, reftime, steps, params, member="000"): + calls.append((Path(root).name, params)) + return _member_dataset(int(Path(root).name)) + + monkeypatch.setattr(verification_metrics, "load_forecast_data", fake_load) + return grib, calls + + +def test_load_ensemble_forecast_concats_members(member_dirs): + grib, calls = member_dirs + fcst = verification_metrics.load_ensemble_forecast( + grib, None, [0, 6], ["T_2M"], members=3 + ) + assert [c[0] for c in calls] == ["001", "002", "003"] + assert fcst["member"].values.tolist() == [1, 2, 3] + assert fcst["T_2M"].sel(member=2).values.tolist() == [[[2.0] * 3] * 2] + assert "member" not in fcst["elevation"].dims + np.testing.assert_array_equal(fcst["elevation"].values, 100.0) + + +def test_load_ensemble_forecast_fails_on_missing_member(member_dirs): + grib, _ = member_dirs + with pytest.raises(FileNotFoundError, match="member 004 of 4"): + verification_metrics.load_ensemble_forecast( + grib, None, [0, 6], ["T_2M"], members=4 + ) + + +def test_load_aligned_data_routes_members(member_dirs, monkeypatch): + grib, _ = member_dirs + truth = xr.Dataset( + {"T_2M": (("time", "values"), np.zeros((2, 3)))}, + coords={ + "time": np.array( + ["2024-06-15T00:00", "2024-06-15T06:00"], dtype="datetime64[ns]" + ) + }, + ) + monkeypatch.setattr(verification_metrics, "load_truth_data", lambda *a: truth) + monkeypatch.setattr(verification_metrics, "map_forecast_to_truth", lambda f, t: f) + args = Namespace( + forecast=grib, reftime=None, steps=[0, 6], truth=None, members=3, member="000" + ) + fcst, _ = verification_metrics.load_aligned_data(args, ["T_2M"]) + assert fcst.sizes["member"] == 3 + + +# --- scoring -------------------------------------------------------------------- + +_REGION = [{"type": "bbox", "name": "all", "bbox": [1.5, 16, 43, 49.5]}] + + +def _stations(values): + """(step, values) station field inside _REGION; values (..., n).""" + values = np.asarray(values, dtype=float) + n = values.shape[-1] + dims = ("member", "step", "values") if values.ndim == 3 else ("step", "values") + coords = { + "step": pd.to_timedelta([6], unit="h"), + "values": np.arange(n), + "longitude": ("values", np.linspace(5.0, 10.0, n)), + "latitude": ("values", np.linspace(46.0, 47.0, n)), + "forecast_reference_time": np.datetime64("2024-06-15T00:00", "ns"), + } + if values.ndim == 3: + coords["member"] = np.arange(1, values.shape[0] + 1) + return xr.Dataset({"T_2M": (dims, values)}, coords=coords) + + +TRUTH = _stations([[1.0, 2.0, 3.0, 4.0]]) +MEMBERS = _stations( + [ + [[0.0, 2.0, 2.0, 6.0]], + [[2.0, 4.0, 3.0, 2.0]], + [[4.0, 0.0, 7.0, 4.0]], + ] +) + + +def _run(fcst, truth=TRUTH, **kwargs): + return verification_metrics.verify_ensemble_run( + fcst, + truth, + "ens", + "truth", + regions=_REGION, + threshold_dict={"T_2M": {"gt": [2.5]}}, + holdout_stations=None, + **kwargs, + ) + + +def test_verify_ensemble_run_scores_the_ensemble_mean(): + out = _run(MEMBERS) + ref = verify( + MEMBERS.mean("member"), + TRUTH, + "ens", + "truth", + regions=_REGION, + threshold_dict={"T_2M": {"gt": [2.5]}}, + num_workers=1, + ) + for name in ("T_2M.BIAS", "T_2M.MSE", "T_2M.MAE", "T_2M.contingency_table"): + xr.testing.assert_allclose(out[name], ref[name]) + # truth statistics are those of verify() + for name in ("T_2M.mean", "T_2M.var", "T_2M.min", "T_2M.max"): + xr.testing.assert_allclose( + out[name].sel(source="truth"), ref[name].sel(source="truth") + ) + for name in ("T_2M.CRPS", "T_2M.ENS_VAR", "T_2M.MEMBER_MSE", "T_2M.rank_hist"): + assert name in out + assert "T_2M.prob_table" in out + assert out.attrs["n_members"] == 3 + assert out["source"].values.tolist() == ["ens", "truth"] + + +def test_verify_ensemble_run_member_mean_statistics(): + out = _run(MEMBERS) + sel = dict(source="ens", region="all", station_group="all") + values = MEMBERS["T_2M"].values[:, 0] # member, values + assert out["T_2M.mean"].sel(**sel).item() == pytest.approx(values.mean()) + assert out["T_2M.var"].sel(**sel).item() == pytest.approx(values.var(axis=1).mean()) + assert out["T_2M.min"].sel(**sel).item() == pytest.approx(values.min(axis=1).mean()) + assert out["T_2M.max"].sel(**sel).item() == pytest.approx(values.max(axis=1).mean()) + + +def test_verify_ensemble_run_missing_member_value(): + """A point where one member is missing makes the scores NaN, as in verify().""" + members = MEMBERS.copy(deep=True) + members["T_2M"][1, 0, 2] = np.nan + out = _run(members) + sel = dict(source="ens", region="all", station_group="all") + for name in ("T_2M.MSE", "T_2M.CRPS", "T_2M.mean", "T_2M.var"): + assert np.isnan(out[name].sel(**sel).item()), name + assert np.isfinite(out["T_2M.mean"].sel(source="truth", region="all").item()) + + +def test_verify_ensemble_run_circular_ensemble_mean(): + """The ensemble mean of wind direction is circular: 350 and 10 deg give 0 deg.""" + members = _stations([[[350.0, 90.0, 180.0, 270.0]], [[10.0, 90.0, 180.0, 270.0]]]) + members = members.rename({"T_2M": "DD_10M"}) + truth = _stations([[0.0, 90.0, 180.0, 270.0]]).rename({"T_2M": "DD_10M"}) + out = verification_metrics.verify_ensemble_run( + members, + truth, + "ens", + "truth", + regions=_REGION, + threshold_dict=None, + holdout_stations=None, + ) + sel = dict(source="ens", region="all", station_group="all") + assert out["DD_10M.MAE"].sel(**sel).item() == pytest.approx(0.0, abs=1e-9) + + +# --- valid mask ----------------------------------------------------------------- + + +def test_valid_mask_of_an_ensemble_is_the_and_over_members(tmp_path, monkeypatch): + members = MEMBERS.copy(deep=True) + members["T_2M"][1, 0, 2] = np.nan # member 2 missing at point 2 + monkeypatch.setattr( + verification_valid_mask, "load_aligned_data", lambda args, p: (members, None) + ) + args = Namespace(params=["T_2M"], output=tmp_path / "mask.nc", members=3) + verification_valid_mask.main(args) + mask = xr.open_dataset(args.output) + assert "member" not in mask.dims + assert mask["T_2M"].values.ravel().tolist() == [True, True, False, True] diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index b979677a..ac0fc7b7 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -687,8 +687,9 @@ def scoremap_runtime(wc): def run_grib_dir(run_id, init_time): """GRIB directory of a run's inference output, as read by downstream rules. - Deterministic runs write ``grib/``; ensemble runs write ``grib/{NNN}/`` per member - and, until ensemble verification lands, downstream rules use member 001 only. + Deterministic runs write ``grib/``; ensemble runs write ``grib/{NNN}/`` per member. + Plots, SAL, score maps and MEC read member 001 only; verification reads all + members (see ``verification_forecast_arg``). """ base = Path(OUT_ROOT) / f"data/runs/{run_id}/{init_time}/grib" if RUN_CONFIGS[run_id].get("ensemble") is not None: @@ -696,11 +697,47 @@ def run_grib_dir(run_id, init_time): return base +def ensemble_members(run_id) -> int | None: + """Number of members of an ensemble run, None for a deterministic run.""" + ensemble = RUN_CONFIGS[run_id].get("ensemble") + return None if ensemble is None else ensemble["members"] + + +def verification_forecast_arg(run_id, init_time): + """``--forecast`` argument of verification_metrics and verification_valid_mask. + + Deterministic runs: the ``grib/`` directory, as a Path. Snakemake does not + record Path params, so this param never triggers a rerun. + Ensemble runs: ``grib/`` followed by ``--members N``; the scripts read every + ``grib/{NNN}/`` and fail on a missing member. Passing the flag in this param + keeps the rules' shell text, and so the deterministic outputs, unchanged. + """ + base = (Path(OUT_ROOT) / f"data/runs/{run_id}/{init_time}/grib").resolve() + members = ensemble_members(run_id) + if members is None: + return base + return f"{base} --members {members}" + + +def ensemble_scaled(base, per_member, unit=None): + """Resource callable: ``base`` for deterministic runs, ``base + N * per_member`` + for an ensemble of N members (``{run_id}`` wildcard); as a string with ``unit`` + appended (e.g. ``"m"`` for a runtime) when ``unit`` is given.""" + + def resource(wc): + members = ensemble_members(wc.run_id) + value = base if members is None else base + members * per_member + return value if unit is None else f"{value}{unit}" + + return resource + + for _run_id, _rc in RUN_CONFIGS.items(): if _rc.get("ensemble") is not None: logging.getLogger("snakemake").warning( - "Run %s is an ensemble of %d members: verification, plots and MEC use " - "member 001 only until ensemble verification is implemented.", + "Run %s is an ensemble of %d members: verification scores the ensemble " + "mean and the ensemble scores over all members; plots, SAL, score maps " + "and MEC use member 001 only.", _run_id, _rc["ensemble"]["members"], ) diff --git a/workflow/rules/verification.smk b/workflow/rules/verification.smk index 97014091..e4ac997f 100644 --- a/workflow/rules/verification.smk +++ b/workflow/rules/verification.smk @@ -63,14 +63,18 @@ rule verification_valid_mask: OUT_ROOT / f"data/runs/{{run_id}}/{{init_time}}/valid_mask_{VALID_MASK_HASH}.nc", log: OUT_ROOT / "logs/verification_valid_mask/{run_id}-{init_time}.log", + # Ensembles: every member is loaded, so memory and runtime grow with N. resources: - cpus_per_task=balfrin_cpus(40_000, actual_cpus=8), - mem_mb=40_000, - runtime="30m", + cpus_per_task=lambda wc: balfrin_cpus( + ensemble_scaled(40_000, 4_000)(wc), actual_cpus=8 + ), + mem_mb=ensemble_scaled(40_000, 4_000), + runtime=ensemble_scaled(30, 3, unit="m"), params: fcst_steps=lambda wc: RUN_CONFIGS[wc.run_id]["steps"], truth=config["truth"]["root"], - grib_out_dir=lambda wc: run_grib_dir(wc.run_id, wc.init_time).resolve(), + # grib/ (+ --members N for an ensemble: valid where every member has a value) + grib_out_dir=lambda wc: verification_forecast_arg(wc.run_id, wc.init_time), mask_params=",".join(COMMON_POINTS_PARAMS), shell: """ @@ -162,10 +166,13 @@ rule verification_metrics: OUT_ROOT / f"data/runs/{{run_id}}/{{init_time}}/verif_{VERIF_HASH}.nc", log: OUT_ROOT / "logs/verification_metrics/{run_id}-{init_time}.log", + # Ensembles: every member is loaded and scored, so memory and runtime grow with N. resources: - cpus_per_task=balfrin_cpus(80_000, actual_cpus=24), - mem_mb=80_000, - runtime="60m", + cpus_per_task=lambda wc: balfrin_cpus( + ensemble_scaled(80_000, 6_000)(wc), actual_cpus=24 + ), + mem_mb=ensemble_scaled(80_000, 6_000), + runtime=ensemble_scaled(60, 6, unit="m"), # wildcard_constraints: # run_id="^" # to avoid ambiguitiy with run_baseline_verif # TODO: implement logic to use experiment name instead of run_id as wildcard @@ -174,7 +181,8 @@ rule verification_metrics: truth=config["truth"]["root"], truth_source_id=f"truth-{TRUTH_HASH}", regions=REGIONS, - grib_out_dir=lambda wc: run_grib_dir(wc.run_id, wc.init_time).resolve(), + # grib/ (+ --members N for an ensemble: ensemble mean and ensemble scores) + grib_out_dir=lambda wc: verification_forecast_arg(wc.run_id, wc.init_time), experiment_params=",".join(EXPERIMENT_PARAMS), threshold_dict=EXPERIMENT_CONFIG.get("thresholds", {}), holdout_flag=lambda wc, input: ( diff --git a/workflow/scripts/verification_metrics.py b/workflow/scripts/verification_metrics.py index 0b6335a1..bc7bda06 100644 --- a/workflow/scripts/verification_metrics.py +++ b/workflow/scripts/verification_metrics.py @@ -5,10 +5,20 @@ from datetime import datetime from pathlib import Path +import numpy as np import pandas as pd import xarray as xr -from verification import verify, apply_lapse_rate_correction_inplace # noqa: E402 +from verification import ( # noqa: E402 + CIRCULAR_PARAMS, + _apply_masks, + _compute_statistics, + _region_and_station_masks, + _too_many_missing, + apply_lapse_rate_correction_inplace, + verify, + verify_ensemble, +) from verification.spatial import map_forecast_to_truth # noqa: E402 from verification.valid_mask import apply_valid_mask, combine_valid_masks # noqa: E402 from data_input import ( @@ -88,16 +98,58 @@ def program_summary_log(args): LOG.info("=" * 80) +def load_ensemble_forecast( + root: Path, reftime: datetime, steps: list[int], params: list[str], members: int +) -> xr.Dataset: + """Load the members ``root/001 .. root/{members:03d}`` of an ensemble run. + + Each member is loaded with ``load_forecast_data`` (same derived and period + parameters as a deterministic run) and the members are concatenated on a + ``member`` dimension (1..N). A missing member is an error. All members get + member 001's ``elevation`` (model orography, the same for every member), so + the lapse-rate correction uses one elevation for the whole ensemble. + """ + datasets = [] + for m in range(1, members + 1): + member_dir = Path(root) / f"{m:03d}" + if not any(member_dir.glob("*.grib")): + raise FileNotFoundError( + f"Ensemble member {m:03d} of {members}: no GRIB files in {member_dir}" + ) + datasets.append(load_forecast_data(member_dir, reftime, steps, params)) + elevation = datasets[0].coords.get("elevation") + datasets = [ds.drop_vars("elevation", errors="ignore") for ds in datasets] + fcst = xr.concat( + datasets, + dim=pd.Index(range(1, members + 1), name="member"), + coords="minimal", + compat="override", + join="exact", + ) + if elevation is not None: + fcst = fcst.assign_coords(elevation=elevation) + return fcst + + def load_aligned_data(args, params: list[str]) -> tuple[xr.Dataset, xr.Dataset]: """Load forecast and truth for ``params``, with the forecast mapped to the truth - points and the truth selected at the forecast valid times.""" + points and the truth selected at the forecast valid times. + + With ``args.members`` set, ``args.forecast`` is the ``grib/`` directory of an + ensemble run and the forecast has a ``member`` dimension (load_ensemble_forecast). + """ # get baseline forecast data now = datetime.now() - fcst = load_forecast_data( - args.forecast, args.reftime, args.steps, params, member=args.member - ) + if getattr(args, "members", None): + fcst = load_ensemble_forecast( + args.forecast, args.reftime, args.steps, params, args.members + ) + else: + fcst = load_forecast_data( + args.forecast, args.reftime, args.steps, params, member=args.member + ) LOG.info( "Loaded forecast data in %s seconds: \n%s", @@ -130,6 +182,101 @@ def load_aligned_data(args, params: list[str]) -> tuple[xr.Dataset, xr.Dataset]: return fcst, truth +def member_mean_statistics( + fcst: xr.Dataset, + truth: xr.Dataset, + regions: list[dict], + holdout_stations: list[str] | None = None, + params: list[str] | None = None, +) -> xr.Dataset: + """Forecast statistics ``{p}.{mean,var,min,max}`` of an ensemble, averaged over + the members (without a ``source`` dim). + + Each member is masked as ``verify`` masks the forecast (regions, station groups, + NaN where too many forecast values are missing; a point is missing if any member + is). The ensemble mean is smoother than any member, so its own variance and + extremes would not describe the forecast fields. + """ + dim = ["x", "y"] if "x" in fcst.dims and "y" in fcst.dims else ["values"] + fcst, truth = xr.align(fcst, truth, join="inner", copy=False) + masks, station_masks = _region_and_station_masks(truth, regions, holdout_stations) + stats = [] + for param in params or list(fcst.data_vars): + if param not in truth.data_vars: + continue + members = fcst[param].where(fcst[param].notnull().all("member")) + too_many_missing = _too_many_missing( + members.isnull().any("member"), + _apply_masks(truth[param], masks, station_masks), + dim, + 0.0, + ) + member_stats = _compute_statistics( + _apply_masks(members, masks, station_masks), + prefix=param + ".", + dim=dim, + circular=param in CIRCULAR_PARAMS, + ) + stats.append( + member_stats.mean("member") + .squeeze("source", drop=True) + .where(~too_many_missing) + ) + return xr.merge(stats).compute() + + +def verify_ensemble_run( + fcst: xr.Dataset, truth: xr.Dataset, source_id: str, truth_source_id: str, **kwargs +) -> xr.Dataset: + """Scores of an ensemble run (``fcst`` with a ``member`` dimension). + + The deterministic scores are those of the ensemble mean (``verify``), under the + run's own ``source``; the ensemble scores (``verify_ensemble``: CRPS, ENS_VAR, + MEMBER_MSE, rank_hist, prob_table) are merged next to them. A point where any + member is missing is missing in the ensemble mean; the ensemble mean of a circular + parameter (wind direction) is the direction of the mean unit vector. The run's forecast statistics + are averaged over the members (``member_mean_statistics``), except for circular + parameters (wind direction), which keep those of the ensemble mean. + """ + n_members = fcst.sizes["member"] + LOG.info("Verifying an ensemble of %d members", n_members) + ens_mean = fcst.mean("member", skipna=False) + for param in CIRCULAR_PARAMS.intersection(fcst.data_vars): + # wind direction: direction of the mean unit vector (350 and 10 deg -> 0 deg) + rad = np.deg2rad(fcst[param]) + ens_mean[param] = np.mod( + np.degrees( + np.arctan2( + np.sin(rad).mean("member", skipna=False), + np.cos(rad).mean("member", skipna=False), + ) + ), + 360.0, + ).assign_attrs(fcst[param].attrs) + results = verify(ens_mean, truth, source_id, truth_source_id, **kwargs) + + stat_params = [p for p in fcst.data_vars if p not in CIRCULAR_PARAMS] + member_stats = member_mean_statistics( + fcst, + truth, + kwargs["regions"], + kwargs.get("holdout_stations"), + params=stat_params, + ) + for name, value in member_stats.data_vars.items(): + current = results[name] + results[name] = ( + xr.where(current["source"] == source_id, value, current) + .transpose(*current.dims) + .assign_attrs(current.attrs) + ) + + ens_scores = verify_ensemble(fcst, truth, source_id, truth_source_id, **kwargs) + out = xr.merge([results, ens_scores], join="outer", compat="no_conflicts") + out.attrs = {**results.attrs, "n_members": n_members} + return out + + def main(args: ScriptConfig): """Main function to verify baseline forecast data.""" @@ -154,15 +301,19 @@ def main(args: ScriptConfig): # compute metrics and statistics now = datetime.now() - results = verify( - fcst, - truth, - args.source_id, - args.truth_source_id, + verify_kwargs = dict( regions=args.regions, threshold_dict=args.threshold_dict, holdout_stations=holdout_stations, ) + if "member" in fcst.dims: + results = verify_ensemble_run( + fcst, truth, args.source_id, args.truth_source_id, **verify_kwargs + ) + else: + results = verify( + fcst, truth, args.source_id, args.truth_source_id, **verify_kwargs + ) LOG.info( "Computed verification metrics in %s seconds", (datetime.now() - now).total_seconds(), @@ -268,6 +419,17 @@ def main(args: ScriptConfig): default="000", help="Ensemble member to load: '000' for control, 'median' for the pre-computed median, 'mean' to average all members, or any 3-digit member ID.", ) + parser.add_argument( + "--members", + type=int, + default=None, + help=( + "Number of members of an ensemble run. --forecast is then the run's grib/ " + "directory holding one sub-directory per member (001, 002, ...); the " + "ensemble mean is scored as the run's deterministic forecast and the " + "ensemble scores are added." + ), + ) parser.add_argument( "--valid_masks", type=Path, diff --git a/workflow/scripts/verification_valid_mask.py b/workflow/scripts/verification_valid_mask.py index ced13ece..ebd32c40 100644 --- a/workflow/scripts/verification_valid_mask.py +++ b/workflow/scripts/verification_valid_mask.py @@ -26,6 +26,10 @@ def main(args): # Loaded exactly as for verification, so the mask matches the scored points. fcst, _ = load_aligned_data(args, args.params) + if "member" in fcst.dims: + # An ensemble has a value where every member has one: the ensemble mean + # (scored by verification_metrics.py) is missing where any member is. + fcst = fcst.mean("member", skipna=False) mask = forecast_valid_mask(fcst, args.params) LOG.info("Valid mask for %s: \n%s", args.params, mask) args.output.parent.mkdir(parents=True, exist_ok=True) @@ -66,5 +70,14 @@ def main(args): default="000", help="Ensemble member to load (see verification_metrics.py).", ) + parser.add_argument( + "--members", + type=int, + default=None, + help=( + "Number of members of an ensemble run (--forecast is then its grib/ " + "directory): the mask is where every member has a value." + ), + ) parser.add_argument("--output", type=Path, required=True, help="Output NetCDF.") main(parser.parse_args()) From 9dd34eac02f1665f55e43ce513e13154f3b48d32 Mon Sep 17 00:00:00 2001 From: Francesco Zanetta Date: Thu, 5 Nov 2026 16:55:55 +0100 Subject: [PATCH 5/7] ensemble review fixes: nan ptype scores, per-ensemble n in plots, plain members flag, lighter member stats, report helpers in verification --- README.md | 4 + src/verification/__init__.py | 18 ++- src/verification/ensemble_report.py | 112 +++++++++++++ tests/unit/test_ensemble_consumers.py | 84 +++++++++- tests/unit/test_ensemble_wiring.py | 150 ++++++++++++++++-- workflow/rules/common.smk | 28 ++-- workflow/rules/report.smk | 1 + workflow/rules/verification.smk | 21 ++- workflow/scripts/ensemble_scores.py | 83 ---------- .../scripts/report_experiment_dashboard.py | 2 +- workflow/scripts/report_scorecard.py | 6 +- workflow/scripts/verification_metrics.py | 53 ++++++- workflow/scripts/verification_plot_metrics.py | 43 +++-- 13 files changed, 452 insertions(+), 153 deletions(-) create mode 100644 src/verification/ensemble_report.py delete mode 100644 workflow/scripts/ensemble_scores.py diff --git a/README.md b/README.md index ed0d199c..3d11bd68 100644 --- a/README.md +++ b/README.md @@ -400,6 +400,10 @@ the run: `SPREAD` and spread-skill ratio `SSR`, the mean score of the single members `MEMBER_RMSE`, Brier scores per threshold, and the counts of the rank histogram and reliability diagram. +- Categorical parameters (`PTYPE`) have no ensemble mean (rain, 1, and snow, 5, + would average to freezing rain, 3), so their contingency scores (`ETS`, `POD`, `FAR`, + ...) are NaN for an ensemble run. The Brier scores and reliability diagrams of their + `eq` thresholds score them instead. - The forecast statistics (`mean`, `var`, `min`, `max`) are averaged over the members. - A point where any member is missing counts as missing for the ensemble (and, with `score_on_common_points`, is left out for every participant). diff --git a/src/verification/__init__.py b/src/verification/__init__.py index cd781630..ef3cf377 100644 --- a/src/verification/__init__.py +++ b/src/verification/__init__.py @@ -351,6 +351,14 @@ def _compute_statistics( return stats +def default_num_workers() -> int: + """Threads for the dask computations: the CPUs this process may run on.""" + try: + return len(os.sched_getaffinity(0)) + except AttributeError: + return max((os.cpu_count() or 6) - 2, 1) + + def _merge_metrics(ds: xr.Dataset, num_workers: int = 4) -> xr.Dataset: out = xr.merge(ds, compat="no_conflicts") if "forecast_reference_time" not in out.dims: @@ -560,10 +568,7 @@ def verify( start = time.time() if num_workers is None: - try: - num_workers = len(os.sched_getaffinity(0)) - except AttributeError: - num_workers = max((os.cpu_count() or 6) - 2, 1) + num_workers = default_num_workers() if dim is None: if "x" in fcst.dims and "y" in fcst.dims: @@ -854,10 +859,7 @@ def verify_ensemble( raise ValueError(f"Ensemble verification needs >= 2 members, got {n_members}.") if num_workers is None: - try: - num_workers = len(os.sched_getaffinity(0)) - except AttributeError: - num_workers = max((os.cpu_count() or 6) - 2, 1) + num_workers = default_num_workers() if dim is None: dim = ["x", "y"] if "x" in fcst.dims and "y" in fcst.dims else ["values"] diff --git a/src/verification/ensemble_report.py b/src/verification/ensemble_report.py new file mode 100644 index 00000000..05dc84bf --- /dev/null +++ b/src/verification/ensemble_report.py @@ -0,0 +1,112 @@ +"""Ensemble scores in the reports: helpers for consumers of aggregated verification +files (dashboard, scorecards, metric plots). + +Importing this module has no side effects (no plotting style, no logging setup). +""" + +import numpy as np +import xarray as xr + +# Variables carrying one of these dims are counts/tables kept for plots, not scores. +NON_SCORE_DIMS = frozenset({"rank", "n_members_yes", "obs"}) + +CRPS_NOTE = "CRPS = MAE for deterministic forecasts" + +# Source id prefix of the truth (``truth-``), which only has statistics. +TRUTH_SOURCE_PREFIX = "truth-" + + +def is_score_var(da: xr.DataArray) -> bool: + """True if the variable is a score that can be tabulated or drawn as a curve.""" + return not (set(da.dims) & NON_SCORE_DIMS) + + +def drop_non_score_vars(ds: xr.Dataset) -> xr.Dataset: + """Drop the rank histogram and probability tables (not scores).""" + drop = [v for v in ds.data_vars if not is_score_var(ds[v])] + return ds.drop_vars(drop) if drop else ds + + +def _ensemble_vars(ds: xr.Dataset) -> list[str]: + """Rank histograms and probability tables: only an ensemble has them.""" + return [v for v in ds.data_vars if ".rank_hist" in v or ".prob_table" in v] + + +def _sources_with_values(ds: xr.Dataset, names: list[str]) -> set: + """Sources with a value (not NaN) in any of the variables *names*.""" + found: set = set() + for name in names: + da = ds[name] + if "source" not in da.dims: + continue + other = [d for d in da.dims if d != "source"] + has = da.notnull().any(other) if other else da.notnull() + found.update(ds["source"].values[has.values].tolist()) + return found + + +def ensemble_sources(ds: xr.Dataset) -> list: + """Sources of *ds* that are ensembles: those with a rank histogram or probability + table (deterministic sources have NaN there after concatenation).""" + found = _sources_with_values(ds, _ensemble_vars(ds)) + return [s for s in ds["source"].values.tolist() if s in found] + + +def is_ensemble(ds: xr.Dataset) -> bool: + """True if a single-participant dataset holds ensemble scores.""" + return any(bool(ds[v].notnull().any()) for v in _ensemble_vars(ds)) + + +def _mae_vars(ds: xr.Dataset) -> list[str]: + return [v for v in ds.data_vars if v.endswith(".MAE")] + + +def add_deterministic_crps(ds: xr.Dataset) -> tuple[xr.Dataset, bool]: + """Multi-source dataset: give deterministic forecasts ``

.CRPS =

.MAE``. + + Only when the dataset contains an ensemble source; otherwise returned unchanged. + Deterministic forecasts are the sources (candidate runs, baselines) with an MAE + and no ensemble scores; the truth source is never filled. The result is for + display only. Returns ``(dataset, filled)``, *filled* True if a source was filled. + """ + ens = ensemble_sources(ds) + if not ens: + return ds, False + with_mae = _sources_with_values(ds, _mae_vars(ds)) + deterministic = [ + s + for s in ds["source"].values.tolist() + if s in with_mae and s not in ens and not str(s).startswith(TRUTH_SOURCE_PREFIX) + ] + if not deterministic: + return ds, False + is_det = xr.DataArray( + np.isin(ds["source"].values, deterministic), + dims="source", + coords={"source": ds["source"]}, + ) + ds = ds.copy() + for mae in _mae_vars(ds): + if "source" not in ds[mae].dims: + continue + name = mae[: -len("MAE")] + "CRPS" + base = ds[name] if name in ds else xr.full_like(ds[mae], np.nan) + out = xr.where(is_det, ds[mae], base) + ds[name] = out.transpose(*ds[mae].dims) + ds[name].attrs = dict(ds[mae].attrs) + return ds, True + + +def fill_crps_from_mae(ds: xr.Dataset) -> tuple[xr.Dataset, bool]: + """Single-participant dataset: if it is deterministic, set ``

.CRPS =

.MAE``.""" + if is_ensemble(ds): + return ds, False + ds = ds.copy() + filled = False + for mae in _mae_vars(ds): + name = mae[: -len("MAE")] + "CRPS" + if name in ds and bool(ds[name].notnull().any()): + continue + ds[name] = ds[mae] + filled = True + return ds, filled diff --git a/tests/unit/test_ensemble_consumers.py b/tests/unit/test_ensemble_consumers.py index 613646f1..4c816e2a 100644 --- a/tests/unit/test_ensemble_consumers.py +++ b/tests/unit/test_ensemble_consumers.py @@ -13,7 +13,7 @@ import report_experiment_dashboard as dash # noqa: E402 import report_scorecard as sc # noqa: E402 import verification_plot_metrics as vpm # noqa: E402 -from ensemble_scores import ( # noqa: E402 +from verification.ensemble_report import ( # noqa: E402 add_deterministic_crps, drop_non_score_vars, ensemble_sources, @@ -47,7 +47,7 @@ def _var(val, extra=None): return dims, arr -def make_file(source, ensemble=False): +def make_file(source, ensemble=False, n=N): ds = xr.Dataset(coords=_coords(source)) for name, val in [("T_2M.RMSE", 1.0), ("T_2M.MAE", 0.8), ("T_2M.BIAS", 0.1)]: ds[name] = _var(val) @@ -61,8 +61,10 @@ def make_file(source, ensemble=False): ("T_2M.BS_gt_0p1", 0.2), ]: ds[name] = _var(val) - ds["T_2M.rank_hist"] = _var(0, {"rank": N + 1}) - ds["T_2M.prob_table_gt_0p1"] = _var(0, {"n_members_yes": N + 1, "obs": 2}) + ds["T_2M.rank_hist"] = _var(0, {"rank": n + 1}) + ds["T_2M.prob_table_gt_0p1"] = _var(0, {"n_members_yes": n + 1, "obs": 2}) + ds = ds.assign_coords(rank=np.arange(n + 1), n_members_yes=np.arange(n + 1)) + ds.attrs["n_members"] = n return ds @@ -158,3 +160,77 @@ def test_plot_ensemble_diagnostics(tmp_path): vpm.plot_ensemble_diagnostics(ds, args) names = {p.name for p in tmp_path.iterdir()} assert names == {"rank_hist_T_2M.png", "reliability_T_2M_gt_0p1.png"} + + +def test_plot_ensemble_diagnostics_per_source_members(tmp_path, monkeypatch): + """Two ensembles of 3 and 5 members: each keeps its own N + 1 ranks, 1/(N+1) + reference line and probabilities k/N (the outer join pads the smaller one).""" + ds = vpm._concat_sources( + [make_file("small", True, n=3), make_file("big", True, n=5)] + ) + ds = ds.sel(station_group="all", drop=True) + assert ds.sizes["rank"] == 6 + figures, close = [], vpm.plt.close + monkeypatch.setattr(vpm.plt, "close", figures.append) + vpm.plot_ensemble_diagnostics(ds, Namespace(output_dir=tmp_path, label_map={})) + rank_fig, rel_fig = figures + lines = rank_fig.axes[0].get_lines() + curves = {ln.get_label(): ln for ln in lines if not ln.get_label().startswith("_")} + assert set(curves) == {"small (N=3)", "big (N=5)"} + assert curves["small (N=3)"].get_xdata().tolist() == [0, 1, 2, 3] + assert curves["big (N=5)"].get_xdata().tolist() == [0, 1, 2, 3, 4, 5] + assert sum(curves["small (N=3)"].get_ydata()) == pytest.approx(1.0) + references = sorted( + ln.get_ydata()[0] for ln in lines if ln.get_label().startswith("_") + ) + assert references == pytest.approx([1 / 6, 1 / 4]) + rel = { + ln.get_label(): ln.get_xdata().tolist() + for ln in rel_fig.axes[0].get_lines() + if not ln.get_label().startswith("_") + } + assert rel["small (N=3)"] == pytest.approx([0, 1 / 3, 2 / 3, 1]) + assert rel["big (N=5)"] == pytest.approx([0, 0.2, 0.4, 0.6, 0.8, 1]) + for fig in figures: + close(fig) + + +def make_truth(source="truth-c305"): + """The truth source of an aggregated file: statistics only, no scores.""" + ds = xr.Dataset(coords=_coords(source)) + ds["T_2M.mean"] = _var(285.0) + return ds + + +def test_ensemble_detected_by_ensemble_variables_not_spread(): + """An ensemble without a finite SPREAD (e.g. only PTYPE, or step 0) is an ensemble.""" + ens = make_file("ens", True) + ens["T_2M.SPREAD"][:] = np.nan + ds = vpm._concat_sources([ens, make_file("det")]) + assert ensemble_sources(ds) == ["ens"] + single, filled = fill_crps_from_mae(ens) + assert not filled and float(single["T_2M.CRPS"].max()) == 0.5 + + +def test_truth_never_gets_crps(): + """Ensemble-only experiment: the truth source is not a deterministic forecast.""" + ds = vpm._concat_sources([make_file("ens", True), make_truth()]) + out, filled = add_deterministic_crps(ds) + assert not filled + assert np.isnan(out["T_2M.CRPS"].sel(source="truth-c305")).all() + # with a deterministic participant, only that one is filled + ds = vpm._concat_sources([make_file("ens", True), make_file("det"), make_truth()]) + out, filled = add_deterministic_crps(ds) + assert filled + assert out["T_2M.CRPS"].sel(source="det").isel(step=0).item() == 0.8 + assert np.isnan(out["T_2M.CRPS"].sel(source="truth-c305")).all() + + +def test_dashboard_ensemble_only_has_no_crps_note(tmp_path): + ens = make_file("ens", True) + truth = make_truth() + truth["n_samples"] = ens["n_samples"] + both = xr.concat([ens, truth], "source", join="outer", data_vars="minimal") + html = _run_dashboard(tmp_path, [both]) + assert 'data-metric="CRPS"' in html + assert "deterministic forecasts" not in html diff --git a/tests/unit/test_ensemble_wiring.py b/tests/unit/test_ensemble_wiring.py index 1ade73d0..1b603027 100644 --- a/tests/unit/test_ensemble_wiring.py +++ b/tests/unit/test_ensemble_wiring.py @@ -21,9 +21,10 @@ COMMON = ROOT / "workflow/rules/common.smk" FUNCTIONS = { + "run_grib_root", "run_grib_dir", "ensemble_members", - "verification_forecast_arg", + "members_flag", "ensemble_scaled", } @@ -52,22 +53,21 @@ def _common_namespace(run_configs, out_root): DET, ENS = RUN_CONFIGS -def test_verification_forecast_arg_deterministic_unchanged(tmp_path): - """Deterministic runs: the same Path as before the ensemble wiring. Snakemake - does not record Path params, so the param never triggers a rerun.""" +def test_verification_reads_grib_root_deterministic_unchanged(tmp_path): + """Deterministic runs: the same grib/ directory as before, and no members flag.""" ns = _common_namespace(RUN_CONFIGS, tmp_path) - arg = ns["verification_forecast_arg"](DET, "202406150000") - assert isinstance(arg, Path) - assert arg == (tmp_path / f"data/runs/{DET}/202406150000/grib").resolve() - assert arg == ns["run_grib_dir"](DET, "202406150000").resolve() + grib = ns["run_grib_root"](DET, "202406150000") + assert grib == tmp_path / f"data/runs/{DET}/202406150000/grib" + assert grib == ns["run_grib_dir"](DET, "202406150000") + assert ns["members_flag"](Namespace(run_id=DET)) == "" -def test_verification_forecast_arg_ensemble(tmp_path): +def test_verification_reads_all_members_of_an_ensemble(tmp_path): ns = _common_namespace(RUN_CONFIGS, tmp_path) - grib = (tmp_path / f"data/runs/{ENS}/202406150000/grib").resolve() - assert ns["verification_forecast_arg"](ENS, "202406150000") == ( - f"{grib} --members 4" + assert ns["run_grib_root"](ENS, "202406150000") == ( + tmp_path / f"data/runs/{ENS}/202406150000/grib" ) + assert ns["members_flag"](Namespace(run_id=ENS)) == "--members 4" # plots, SAL, score maps and MEC keep member 001 assert ns["run_grib_dir"](ENS, "202406150000") == ( tmp_path / f"data/runs/{ENS}/202406150000/grib/001" @@ -76,11 +76,11 @@ def test_verification_forecast_arg_ensemble(tmp_path): def test_ensemble_scaled_resources(tmp_path): ns = _common_namespace(RUN_CONFIGS, tmp_path) - mem = ns["ensemble_scaled"](80_000, 6_000) + mem = ns["ensemble_scaled"](80_000, 2_000) runtime = ns["ensemble_scaled"](60, 6, unit="m") assert mem(Namespace(run_id=DET)) == 80_000 assert runtime(Namespace(run_id=DET)) == "60m" - assert mem(Namespace(run_id=ENS)) == 104_000 + assert mem(Namespace(run_id=ENS)) == 88_000 assert runtime(Namespace(run_id=ENS)) == "84m" @@ -286,3 +286,125 @@ def test_valid_mask_of_an_ensemble_is_the_and_over_members(tmp_path, monkeypatch mask = xr.open_dataset(args.output) assert "member" not in mask.dims assert mask["T_2M"].values.ravel().tolist() == [True, True, False, True] + + +def test_verify_ensemble_run_categorical_has_no_ensemble_mean_scores(): + """PTYPE: the mean of class codes is no class (rain 1 + snow 5 -> 3, freezing + rain), so the run's deterministic scores are NaN; prob_table scores the ensemble.""" + members = _stations([[[1.0, 5.0, 0.0, 1.0]], [[5.0, 1.0, 0.0, 1.0]]]) + members = members.rename({"T_2M": "PTYPE"}) + truth = _stations([[3.0, 1.0, 0.0, 1.0]]).rename({"T_2M": "PTYPE"}) + out = verification_metrics.verify_ensemble_run( + members, + truth, + "ens", + "truth", + regions=_REGION, + threshold_dict={"PTYPE": {"eq": [1, 3, 5]}}, + holdout_stations=None, + ) + sel = dict(source="ens", region="all", station_group="all") + for name in ("PTYPE.BIAS", "PTYPE.MSE", "PTYPE.MAE", "PTYPE.CORR"): + assert np.isnan(out[name].sel(**sel)).all(), name + table = out["PTYPE.contingency_table"] + assert "source" in table.dims + assert np.isnan(table.sel(**sel)).all() + # the ensemble's categorical score: 2 members, eq_3 never forecast but observed once + prob = out["PTYPE.prob_table"].sel(**sel) + assert prob.sel({"PTYPE.threshold": "eq_3", "n_members_yes": 0, "obs": 1}) == 1 + assert int(prob.sum()) == 3 * 4 + # forecast statistics (averaged over members) and truth statistics are kept + assert out["PTYPE.mean"].sel(**sel).item() == pytest.approx(14 / 8) + assert np.isfinite(out["PTYPE.mean"].sel(source="truth", region="all").item()) + + +def test_verify_ensemble_run_keeps_continuous_contingency_tables(): + out = _run(MEMBERS) + sel = dict(source="ens", region="all", station_group="all") + assert np.isfinite(out["T_2M.contingency_table"].sel(**sel)).all() + + +def _reference_member_mean_statistics(fcst, truth, regions, holdout_stations, params): + """member_mean_statistics before the memory fix (all params in one compute).""" + from verification import ( + CIRCULAR_PARAMS, + _apply_masks, + _compute_statistics, + _region_and_station_masks, + _too_many_missing, + ) + + dim = ["values"] + fcst, truth = xr.align(fcst, truth, join="inner", copy=False) + masks, station_masks = _region_and_station_masks(truth, regions, holdout_stations) + stats = [] + for param in params: + members = fcst[param].where(fcst[param].notnull().all("member")) + too_many_missing = _too_many_missing( + members.isnull().any("member"), + _apply_masks(truth[param], masks, station_masks), + dim, + 0.0, + ) + member_stats = _compute_statistics( + _apply_masks(members, masks, station_masks), + prefix=param + ".", + dim=dim, + circular=param in CIRCULAR_PARAMS, + ) + stats.append( + member_stats.mean("member") + .squeeze("source", drop=True) + .where(~too_many_missing) + ) + return xr.merge(stats).compute() + + +def test_member_mean_statistics_unchanged_by_chunking(): + """Members chunked one by one, one param at a time: same statistics as before.""" + rng = np.random.default_rng(1) + n_members, n_steps, n = 4, 3, 50 + values = rng.normal(280.0, 5.0, (n_members, n_steps, n)) + values[2, 1, 7] = np.nan # one member missing at one point and step + coords = { + "member": np.arange(1, n_members + 1), + "step": pd.to_timedelta(np.arange(n_steps) * 6, unit="h"), + "values": [f"S{i:02d}" for i in range(n)], + "longitude": ("values", np.linspace(5.0, 10.0, n)), + "latitude": ("values", np.linspace(45.9, 47.5, n)), + } + fcst = xr.Dataset( + { + "T_2M": (("member", "step", "values"), values), + "TD_2M": (("member", "step", "values"), values - 5.0), + }, + coords=coords, + ).chunk({"step": 1}) + truth = fcst.isel(member=0, drop=True).compute() + 0.5 + regions = _REGION + [{"type": "bbox", "name": "west", "bbox": [5.0, 7.5, 43, 49.5]}] + holdout = ["S00", "S03", "S10"] + params = ["T_2M", "TD_2M"] + new = verification_metrics.member_mean_statistics( + fcst, truth, regions, holdout, params=params, num_workers=2 + ) + ref = _reference_member_mean_statistics(fcst, truth, regions, holdout, params) + xr.testing.assert_allclose(new, ref, rtol=1e-12) + # S07 is missing in one member at step 1: NaN where it counts, not in holdout + assert np.isnan(new["T_2M.mean"].isel(step=1).sel(station_group="all")).all() + assert np.isfinite(new["T_2M.mean"].isel(step=1).sel(station_group="holdout")).all() + assert np.isfinite(new["T_2M.mean"].isel(step=0)).all() + + +def test_verification_rules_pass_members_flag(): + """verification_metrics and verification_valid_mask read grib/ and pass the members + flag as a plain param (empty for deterministic runs).""" + text = (ROOT / "workflow/rules/verification.smk").read_text() + for rule in ("verification_valid_mask", "verification_metrics"): + body = text.split(f"rule {rule}:", 1)[1].split("\nrule ", 1)[0] + assert ( + "grib_out_dir=lambda wc: run_grib_root(wc.run_id, wc.init_time).resolve()" + in body + ) + assert "members_flag=members_flag" in body + assert "--forecast {params.grib_out_dir}" in body + assert "{params.members_flag}" in body diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index ac0fc7b7..0dd832b1 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -684,14 +684,20 @@ def scoremap_runtime(wc): return max(60, math.ceil(1.5 * len(REFTIMES) * seconds / 60)) +def run_grib_root(run_id, init_time): + """``grib/`` directory of a run's inference output (one sub-directory per member + for an ensemble). Verification reads it, with ``members_flag`` for an ensemble.""" + return Path(OUT_ROOT) / f"data/runs/{run_id}/{init_time}/grib" + + def run_grib_dir(run_id, init_time): """GRIB directory of a run's inference output, as read by downstream rules. Deterministic runs write ``grib/``; ensemble runs write ``grib/{NNN}/`` per member. Plots, SAL, score maps and MEC read member 001 only; verification reads all - members (see ``verification_forecast_arg``). + members (``run_grib_root`` and ``members_flag``). """ - base = Path(OUT_ROOT) / f"data/runs/{run_id}/{init_time}/grib" + base = run_grib_root(run_id, init_time) if RUN_CONFIGS[run_id].get("ensemble") is not None: return base / "001" return base @@ -703,20 +709,12 @@ def ensemble_members(run_id) -> int | None: return None if ensemble is None else ensemble["members"] -def verification_forecast_arg(run_id, init_time): - """``--forecast`` argument of verification_metrics and verification_valid_mask. - - Deterministic runs: the ``grib/`` directory, as a Path. Snakemake does not - record Path params, so this param never triggers a rerun. - Ensemble runs: ``grib/`` followed by ``--members N``; the scripts read every - ``grib/{NNN}/`` and fail on a missing member. Passing the flag in this param - keeps the rules' shell text, and so the deterministic outputs, unchanged. +def members_flag(wildcards): + """``--members N`` for an ensemble run (``{run_id}`` wildcard), else empty: the + verification scripts then read every ``grib/{NNN}/`` and fail on a missing member. """ - base = (Path(OUT_ROOT) / f"data/runs/{run_id}/{init_time}/grib").resolve() - members = ensemble_members(run_id) - if members is None: - return base - return f"{base} --members {members}" + members = ensemble_members(wildcards.run_id) + return "" if members is None else f"--members {members}" def ensemble_scaled(base, per_member, unit=None): diff --git a/workflow/rules/report.smk b/workflow/rules/report.smk index 7dea3a93..5a4635dd 100644 --- a/workflow/rules/report.smk +++ b/workflow/rules/report.smk @@ -22,6 +22,7 @@ def make_header_text(): rule report_experiment_dashboard: input: "src/verification/__init__.py", + "src/verification/ensemble_report.py", script="workflow/scripts/report_experiment_dashboard.py", verif=EXPERIMENT_PARTICIPANTS.values(), template="resources/report/dashboard/template.html.jinja2", diff --git a/workflow/rules/verification.smk b/workflow/rules/verification.smk index e4ac997f..7fb5f21e 100644 --- a/workflow/rules/verification.smk +++ b/workflow/rules/verification.smk @@ -66,21 +66,23 @@ rule verification_valid_mask: # Ensembles: every member is loaded, so memory and runtime grow with N. resources: cpus_per_task=lambda wc: balfrin_cpus( - ensemble_scaled(40_000, 4_000)(wc), actual_cpus=8 + ensemble_scaled(40_000, 1_000)(wc), actual_cpus=8 ), - mem_mb=ensemble_scaled(40_000, 4_000), + mem_mb=ensemble_scaled(40_000, 1_000), runtime=ensemble_scaled(30, 3, unit="m"), params: fcst_steps=lambda wc: RUN_CONFIGS[wc.run_id]["steps"], truth=config["truth"]["root"], - # grib/ (+ --members N for an ensemble: valid where every member has a value) - grib_out_dir=lambda wc: verification_forecast_arg(wc.run_id, wc.init_time), + grib_out_dir=lambda wc: run_grib_root(wc.run_id, wc.init_time).resolve(), + # ensembles: all members (grib/{NNN}/) are read + members_flag=members_flag, mask_params=",".join(COMMON_POINTS_PARAMS), shell: """ export ECCODES_DEFINITION_PATH=$(realpath .venv/share/eccodes-cosmo-resources/definitions) uv run {input.script} \ --forecast {params.grib_out_dir} \ + {params.members_flag} \ --truth "{params.truth}" \ --reftime {wildcards.init_time} \ --steps "{params.fcst_steps}" \ @@ -169,9 +171,9 @@ rule verification_metrics: # Ensembles: every member is loaded and scored, so memory and runtime grow with N. resources: cpus_per_task=lambda wc: balfrin_cpus( - ensemble_scaled(80_000, 6_000)(wc), actual_cpus=24 + ensemble_scaled(80_000, 2_000)(wc), actual_cpus=24 ), - mem_mb=ensemble_scaled(80_000, 6_000), + mem_mb=ensemble_scaled(80_000, 2_000), runtime=ensemble_scaled(60, 6, unit="m"), # wildcard_constraints: # run_id="^" # to avoid ambiguitiy with run_baseline_verif @@ -181,8 +183,9 @@ rule verification_metrics: truth=config["truth"]["root"], truth_source_id=f"truth-{TRUTH_HASH}", regions=REGIONS, - # grib/ (+ --members N for an ensemble: ensemble mean and ensemble scores) - grib_out_dir=lambda wc: verification_forecast_arg(wc.run_id, wc.init_time), + grib_out_dir=lambda wc: run_grib_root(wc.run_id, wc.init_time).resolve(), + # ensembles: all members (grib/{NNN}/) are read + members_flag=members_flag, experiment_params=",".join(EXPERIMENT_PARAMS), threshold_dict=EXPERIMENT_CONFIG.get("thresholds", {}), holdout_flag=lambda wc, input: ( @@ -199,6 +202,7 @@ rule verification_metrics: export ECCODES_DEFINITION_PATH=$(realpath .venv/share/eccodes-cosmo-resources/definitions) uv run {input.script} \ --forecast {params.grib_out_dir} \ + {params.members_flag} \ --truth "{params.truth}" \ --reftime {wildcards.init_time} \ --steps "{params.fcst_steps}" \ @@ -260,6 +264,7 @@ use rule verification_metrics_aggregation as verification_metrics_aggregation_ba rule verification_metrics_plot: input: "src/verification/__init__.py", + "src/verification/ensemble_report.py", script="workflow/scripts/verification_plot_metrics.py", verif=list(EXPERIMENT_PARTICIPANTS.values()), output: diff --git a/workflow/scripts/ensemble_scores.py b/workflow/scripts/ensemble_scores.py deleted file mode 100644 index 362461f2..00000000 --- a/workflow/scripts/ensemble_scores.py +++ /dev/null @@ -1,83 +0,0 @@ -"""Helpers for consumers of aggregated verification files (ensemble scores). - -Shared by the dashboard, the scorecards and the metric plots. Importing this module -has no side effects (no plotting style, no logging setup). -""" - -import numpy as np -import xarray as xr - -# Variables carrying one of these dims are counts/tables kept for plots, not scores. -NON_SCORE_DIMS = frozenset({"rank", "n_members_yes", "obs"}) - -CRPS_NOTE = "CRPS = MAE for deterministic forecasts" - - -def is_score_var(da: xr.DataArray) -> bool: - """True if the variable is a score that can be tabulated or drawn as a curve.""" - return not (set(da.dims) & NON_SCORE_DIMS) - - -def drop_non_score_vars(ds: xr.Dataset) -> xr.Dataset: - """Drop the rank histogram and probability tables (not scores).""" - drop = [v for v in ds.data_vars if not is_score_var(ds[v])] - return ds.drop_vars(drop) if drop else ds - - -def _spread_vars(ds: xr.Dataset) -> list[str]: - return [v for v in ds.data_vars if v.endswith(".SPREAD")] - - -def ensemble_sources(ds: xr.Dataset) -> list: - """Sources of *ds* that carry ensemble scores (any finite SPREAD).""" - found: set = set() - for v in _spread_vars(ds): - if "source" not in ds[v].dims: - continue - other = [d for d in ds[v].dims if d != "source"] - has = ds[v].notnull().any(other) if other else ds[v].notnull() - found.update(ds["source"].values[has.values].tolist()) - return [s for s in ds["source"].values.tolist() if s in found] - - -def add_deterministic_crps(ds: xr.Dataset) -> tuple[xr.Dataset, bool]: - """Multi-source dataset: give deterministic sources ``

.CRPS =

.MAE``. - - Only when the dataset contains an ensemble source; otherwise returned unchanged. - The result is for display only. Returns ``(dataset, filled)``. - """ - ens = ensemble_sources(ds) - if not ens: - return ds, False - is_det = xr.DataArray( - ~np.isin(ds["source"].values, ens), - dims="source", - coords={"source": ds["source"]}, - ) - ds = ds.copy() - filled = False - for mae in [v for v in ds.data_vars if v.endswith(".MAE")]: - if "source" not in ds[mae].dims: - continue - name = mae[: -len("MAE")] + "CRPS" - base = ds[name] if name in ds else xr.full_like(ds[mae], np.nan) - out = xr.where(is_det, ds[mae], base) - ds[name] = out.transpose(*ds[mae].dims) - ds[name].attrs = dict(ds[mae].attrs) - filled = True - return ds, filled - - -def fill_crps_from_mae(ds: xr.Dataset) -> tuple[xr.Dataset, bool]: - """Single-participant dataset: if it is deterministic, set ``

.CRPS =

.MAE``.""" - if any(bool(ds[v].notnull().any()) for v in _spread_vars(ds)): - return ds, False - ds = ds.copy() - filled = False - for mae in [v for v in ds.data_vars if v.endswith(".MAE")]: - name = mae[: -len("MAE")] + "CRPS" - if name in ds and bool(ds[name].notnull().any()): - continue - ds[name] = ds[mae] - filled = True - return ds, filled diff --git a/workflow/scripts/report_experiment_dashboard.py b/workflow/scripts/report_experiment_dashboard.py index c45d13a1..532028a4 100644 --- a/workflow/scripts/report_experiment_dashboard.py +++ b/workflow/scripts/report_experiment_dashboard.py @@ -14,7 +14,7 @@ _ensure_unique_lead_time, _select_best_sources, ) -from ensemble_scores import CRPS_NOTE, add_deterministic_crps, is_score_var +from verification.ensemble_report import CRPS_NOTE, add_deterministic_crps, is_score_var from verification import decode_metric LOG = logging.getLogger(__name__) diff --git a/workflow/scripts/report_scorecard.py b/workflow/scripts/report_scorecard.py index 40e0ba04..a5637906 100644 --- a/workflow/scripts/report_scorecard.py +++ b/workflow/scripts/report_scorecard.py @@ -20,7 +20,11 @@ import matplotlib.pyplot as plt from matplotlib.transforms import ScaledTranslation -from ensemble_scores import CRPS_NOTE, drop_non_score_vars, fill_crps_from_mae +from verification.ensemble_report import ( + CRPS_NOTE, + drop_non_score_vars, + fill_crps_from_mae, +) from verification import decode_metric LOG = logging.getLogger(__name__) diff --git a/workflow/scripts/verification_metrics.py b/workflow/scripts/verification_metrics.py index bc7bda06..29cb3480 100644 --- a/workflow/scripts/verification_metrics.py +++ b/workflow/scripts/verification_metrics.py @@ -10,12 +10,14 @@ import xarray as xr from verification import ( # noqa: E402 + CATEGORICAL_PARAMS, CIRCULAR_PARAMS, _apply_masks, _compute_statistics, _region_and_station_masks, _too_many_missing, apply_lapse_rate_correction_inplace, + default_num_workers, verify, verify_ensemble, ) @@ -188,6 +190,7 @@ def member_mean_statistics( regions: list[dict], holdout_stations: list[str] | None = None, params: list[str] | None = None, + num_workers: int | None = None, ) -> xr.Dataset: """Forecast statistics ``{p}.{mean,var,min,max}`` of an ensemble, averaged over the members (without a ``source`` dim). @@ -196,7 +199,12 @@ def member_mean_statistics( NaN where too many forecast values are missing; a point is missing if any member is). The ensemble mean is smoother than any member, so its own variance and extremes would not describe the forecast fields. + + Memory: the masked members have region x station_group x points values per member + and step. Each member is its own dask chunk and the parameters are computed one + at a time with *num_workers* threads, so a chunk is no larger than in ``verify``. """ + num_workers = num_workers or default_num_workers() dim = ["x", "y"] if "x" in fcst.dims and "y" in fcst.dims else ["values"] fcst, truth = xr.align(fcst, truth, join="inner", copy=False) masks, station_masks = _region_and_station_masks(truth, regions, holdout_stations) @@ -204,25 +212,59 @@ def member_mean_statistics( for param in params or list(fcst.data_vars): if param not in truth.data_vars: continue - members = fcst[param].where(fcst[param].notnull().all("member")) + members = fcst[param] + members = members.where(members.notnull().all("member")) too_many_missing = _too_many_missing( members.isnull().any("member"), _apply_masks(truth[param], masks, station_masks), dim, 0.0, ) + if members.chunks: + members = members.chunk({"member": 1}) member_stats = _compute_statistics( _apply_masks(members, masks, station_masks), prefix=param + ".", dim=dim, circular=param in CIRCULAR_PARAMS, ) + member_stats = member_stats.assign(too_many_missing=too_many_missing).compute( + num_workers=num_workers, scheduler="threads" + ) + too_many_missing = member_stats["too_many_missing"] stats.append( - member_stats.mean("member") + member_stats.drop_vars("too_many_missing") + .mean("member") .squeeze("source", drop=True) .where(~too_many_missing) ) - return xr.merge(stats).compute() + return xr.merge(stats) + + +def _nan_categorical_scores(results: xr.Dataset, params, source_id: str) -> xr.Dataset: + """NaN for the scores of the categorical params (class codes) under *source_id*. + + The ensemble mean of class codes is no class (rain 1 and snow 5 average to 3, + freezing rain), so the contingency tables of the ensemble mean are meaningless. + The forecast statistics (averaged over the members) and the truth are kept. + """ + statistics = ("mean", "var", "min", "max") + is_source = results["source"] == source_id + for param in CATEGORICAL_PARAMS.intersection(params): + for name in list(results.data_vars): + if not name.startswith(param + "."): + continue + if name.split(".", 1)[1] in statistics: + continue + current = results[name] + if "source" not in current.dims: + continue + results[name] = ( + current.where(~is_source) + .transpose(*current.dims) + .assign_attrs(current.attrs) + ) + return results def verify_ensemble_run( @@ -237,6 +279,9 @@ def verify_ensemble_run( parameter (wind direction) is the direction of the mean unit vector. The run's forecast statistics are averaged over the members (``member_mean_statistics``), except for circular parameters (wind direction), which keep those of the ensemble mean. + Categorical parameters (precipitation type) have no ensemble mean: their + deterministic scores (contingency tables) are NaN; their ``prob_table`` is the + ensemble's categorical score. """ n_members = fcst.sizes["member"] LOG.info("Verifying an ensemble of %d members", n_members) @@ -254,6 +299,7 @@ def verify_ensemble_run( 360.0, ).assign_attrs(fcst[param].attrs) results = verify(ens_mean, truth, source_id, truth_source_id, **kwargs) + results = _nan_categorical_scores(results, fcst.data_vars, source_id) stat_params = [p for p in fcst.data_vars if p not in CIRCULAR_PARAMS] member_stats = member_mean_statistics( @@ -262,6 +308,7 @@ def verify_ensemble_run( kwargs["regions"], kwargs.get("holdout_stations"), params=stat_params, + num_workers=kwargs.get("num_workers"), ) for name, value in member_stats.data_vars.items(): current = results[name] diff --git a/workflow/scripts/verification_plot_metrics.py b/workflow/scripts/verification_plot_metrics.py index b051a6ad..de86a7ba 100644 --- a/workflow/scripts/verification_plot_metrics.py +++ b/workflow/scripts/verification_plot_metrics.py @@ -7,7 +7,7 @@ import matplotlib.pyplot as plt import pandas as pd import xarray as xr -from ensemble_scores import ( +from verification.ensemble_report import ( add_deterministic_crps, ensemble_sources, is_score_var, @@ -198,19 +198,28 @@ def _subset_df(df): plot_ensemble_diagnostics(ds, args) -def _counts_for_plot(da: xr.DataArray, source, region, keep: list[str]) -> xr.DataArray: - """Counts of one source at the default stratum, summed over everything but *keep*.""" +def _counts_for_plot( + da: xr.DataArray, source, region, keep: list[str], bins: str +) -> xr.DataArray: + """Counts of one source at the default stratum, summed over everything but *keep*. + + The participants are concatenated with an outer join, so an ensemble of N members + has NaN in the *bins* (``rank`` or ``n_members_yes``) above N when another ensemble + has more members: those bins are dropped, and the counts keep the source's N + 1 bins. + """ da = da.sel(source=source) for dim, val in (("region", region), ("season", "all"), ("init_hour", -999)): if dim in da.dims: da = da.sel({dim: val}, drop=True) - return da.sum([d for d in da.dims if d not in keep], skipna=True) + da = da.sum([d for d in da.dims if d not in keep], skipna=True, min_count=1) + return da.dropna(bins, how="all") def plot_ensemble_diagnostics(ds: xr.Dataset, args: Namespace) -> None: """Rank histograms and reliability diagrams, for ensemble participants only. Counts are summed over lead times, at the default region, all seasons and init hours. + Each ensemble has its own number of members N (N + 1 ranks, probabilities k/N). """ ens_sources = ensemble_sources(ds) if not ens_sources: @@ -224,20 +233,21 @@ def plot_ensemble_diagnostics(ds: xr.Dataset, args: Namespace) -> None: fig, ax = plt.subplots(figsize=(10, 6)) n_plotted = 0 for source in ens_sources: - counts = _counts_for_plot(ds[name], source, region, ["rank"]) + counts = _counts_for_plot(ds[name], source, region, ["rank"], "rank") total = float(counts.sum()) if not total > 0: continue - ax.plot( + n_ranks = counts.sizes["rank"] + (line,) = ax.plot( counts["rank"].values, - counts.values / total, + counts.fillna(0).values / total, marker="o", - label=label(source), + label=f"{label(source)} (N={n_ranks - 1})", ) + # flat histogram of this ensemble: 1 / (N + 1) + ax.axhline(1 / n_ranks, color=line.get_color(), ls="--", lw=1) n_plotted += 1 if n_plotted: - n_ranks = ds.sizes["rank"] - ax.axhline(1 / n_ranks, color="grey", ls="--", lw=1) ax.set( title=f"Rank histogram - {param} - {region}", xlabel="Rank of observation among members", @@ -254,15 +264,16 @@ def plot_ensemble_diagnostics(ds: xr.Dataset, args: Namespace) -> None: ) n_plotted = 0 for source in ens_sources: - tab = _counts_for_plot(ds[name], source, region, ["n_members_yes", "obs"]) - n_yes = tab["n_members_yes"].values - n_mem = int(n_yes.max()) + tab = _counts_for_plot( + ds[name], source, region, ["n_members_yes", "obs"], "n_members_yes" + ).fillna(0) n_bin = tab.sum("obs") - freq = tab.sel(obs=1) / n_bin.where(n_bin > 0) if not bool(n_bin.sum() > 0): continue - prob = n_yes / n_mem - ax.plot(prob, freq.values, marker="o", label=label(source)) + n_mem = tab.sizes["n_members_yes"] - 1 + freq = tab.sel(obs=1) / n_bin.where(n_bin > 0) + prob = tab["n_members_yes"].values / n_mem + ax.plot(prob, freq.values, marker="o", label=f"{label(source)} (N={n_mem})") axc.plot(prob, n_bin.values, marker="o") n_plotted += 1 if n_plotted: From b45eba0c84a51310336a16226115aec13bf671d8 Mon Sep 17 00:00:00 2001 From: Francesco Zanetta Date: Thu, 5 Nov 2026 16:56:49 +0100 Subject: [PATCH 6/7] readme: crps in scorecard metrics --- README.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 3d11bd68..e546b06c 100644 --- a/README.md +++ b/README.md @@ -114,7 +114,8 @@ experiment: # Stratification dimension to use as scorecard columns (e.g. region, season). stratification: region # Variables and metrics as scorecard rows. Format: VAR:METRIC1,METRIC2,... - # Supported metrics: RMSE, R2, ETS, POD, FAR (categorical metrics require thresholds). + # Supported metrics: RMSE, R2, ETS, POD, FAR (categorical metrics require thresholds), + # CRPS (ensemble runs; a deterministic participant uses its MAE). variables: - "T_2M:RMSE,R2" - "TOT_PREC:RMSE,ETS" From 8d71a0792544cb0ea8006adc8251bdc4cc467f65 Mon Sep 17 00:00:00 2001 From: Francesco Zanetta Date: Thu, 5 Nov 2026 17:24:58 +0100 Subject: [PATCH 7/7] compare verify golden files up to last-bit float differences --- tests/unit/test_verify_regression.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/tests/unit/test_verify_regression.py b/tests/unit/test_verify_regression.py index 34621d98..89259ec4 100644 --- a/tests/unit/test_verify_regression.py +++ b/tests/unit/test_verify_regression.py @@ -1,4 +1,5 @@ -"""verify() output must not change when its internals are refactored. +"""verify() output must not change when its internals are refactored (up to last-bit +floating-point differences between machines). The golden files under ``tests/unit/data/`` were written by verify() before the region/station-group mask construction and the missing-forecast rule were factored out @@ -147,13 +148,20 @@ def test_verify_output_unchanged(name): result = verify(**cases()[name]) expected = xr.load_dataset(DATA / f"verify_golden_{name}.nc") - xr.testing.assert_identical(result, expected) + # Same structure, attributes and NaN positions; values equal up to the last bits, + # which differ between machines for sin/cos (e.g. the circular variance of DD_10M). + xr.testing.assert_allclose(result, expected, rtol=1e-12, atol=0) + assert result.attrs == expected.attrs + assert set(result.variables) == set(expected.variables) for var in expected.variables: assert result[var].dims == expected[var].dims, var + assert result[var].attrs == expected[var].attrs, var if expected[var].dtype.kind in "fiub": assert result[var].dtype == expected[var].dtype, var - # Bit for bit, NaN at the same places - np.testing.assert_array_equal(result[var].values, expected[var].values) + if expected[var].dtype.kind == "f": + np.testing.assert_array_equal( + np.isnan(result[var].values), np.isnan(expected[var].values), var + ) if __name__ == "__main__":