Fix grouped probability and durational exceedance statistics

[rafa-guedes]

Motivation

A monthly-grouped range_probability was needed for the NZ Parent hindcast stats pipeline (config/hindcast/prax/stats/stats_nzpar_cawthron.yml). Adding group: month to the config appeared to work but silently produced the wrong numbers, which turned up three further defects in the neighbouring exceedance code.

What's here

Three commits, split by kind of change. Each passes the suite on its own, so the history bisects cleanly.

1. range_probability honours group (feature)

group is a first-class field on CallConfig and the pipeline always forwards it, but range_probability absorbed it into **kwargs and ignored it. The call still succeeded, and pipeline._apply still appended the group to the output variable name — so a group: month config wrote the ungrouped all-record probability under a name like hs_gt_1p5_range_probability_month. Wrong result, convincing name, no error.

NaN semantics are preserved by taking the denominator from a grouped count() rather than a boolean mean(), so all-NaN land cells stay NaN instead of becoming a real 0.0. That matters here: in the target store 14% of the domain is NaN land while 126 wet cells legitimately hold exactly 0.0, and a boolean mean would merge the two.

The denominator depends only on the variable, so it is now computed once and shared across every range defined on it — a config with three ranges on hs builds one grouped count instead of three.

2. Durational exceedance is actually grouped (bug fix)

With both duration and group set, the run-length reduction happened inside apply_ufunc before group was consulted, and the ungrouped whole-record value was then expand_dims'd across the group keys. Every month received the identical annual number.

Run detection now returns a per-timestep mask computed on the continuous record, with grouping as a plain groupby().mean() afterwards. Detecting before grouping is deliberate: grouping first would both truncate an event at the period boundary (a 72h storm split across two months fails a 72h threshold in both and disappears) and splice non-adjacent periods, since a month grouper places every January end to end and would invent multi-year storms at each year boundary. There are regression tests for both.

The attribution rule — each qualifying sample counts towards the group it actually falls in — is now documented in the exceedance docstring, along with its consequences. It matches DHI's published persistence methodology for fraction-of-time metrics and is the only convention under which the group values recompose to the ungrouped result (verified exact to 0.00e+00 over 20 years of synthetic 3-hourly data at durations from 0h to 48h).

Two further defects fixed in the same machinery:

• Runs touching either end of the record were dropped. Detection used scipy.signal.find_peaks with plateau_size, which only reports plateaus bounded by lower values on both sides — so a record wholly above the threshold returned 0.0 instead of 1.0. Checked against the old implementation directly: for a run at the start, old 0.000 vs new 0.667; only mid-record runs agreed. Now based on padded-diff run boundaries, which is exact at the edges and drops the scipy import from this module. Existing durational exceedance values may increase slightly.
• The timestep assumed nanosecond time. float(np.diff(time)[0]) / 3.6e12 raises TypeError under NumPy 2 for any other datetime64 resolution, and would silently mis-scale if it did not. Now goes through pd.Timedelta.

3. Fractional thresholds named 1p5, not 1.5 (breaking)

{var}_{threshold:g} produced hs_1.5, and with the pipeline suffixes hs_1.5_exceedance_month — a dot makes the name unreachable by attribute access and awkward in path-like contexts.

BREAKING: configs with fractional thresholds will write new variable names. Integral thresholds are unaffected (:g already drops trailing zeros, so 2.0 has always produced hs_2), and no config in this repo changes behaviour — the Douglas climatology uses integer thresholds. Note that when appending to an existing Zarr store the old dotted variables are not removed automatically, so a re-run would leave both hs_1.5_exceedance and hs_1p5_exceedance in place.

Negative thresholds still emit a minus sign (hs_-1), which has the same identifier problem; that convention is left open deliberately.

Testing

• Full suite: 207 passed (200 before this branch; 7 tests added).
• Per commit: 201 → 206 → 207, each green in isolation.
• Grouped range_probability verified against manual per-month computation, NaN-aware, numpy- and dask-backed paths identical.
• End-to-end append into a Zarr store mirroring the production one: pre-existing variables byte-identical, new month dimension added, consolidated metadata picks up everything.

Notes for the reviewer

• CHANGELOG.md gets an Unreleased section rather than a version bump — the changelog is already behind (__version__ is 2.5.0, top entry is 2.3.0), so the version call is yours.
• The conventions in commit 2 were chosen after reviewing how xclim (resample_before_rl), icclim, CDO, climdex and ECA&D handle spells crossing a period boundary, and the metocean persistence literature (DNV-RP-H103 §8.5.1.4, DHI's persistence methodology). The short version: the correct convention follows from the metric — fraction-of-time metrics time-share, event-count metrics assign the whole event to its start period. exceedance returns a fraction of time. An event-based weather-window product would want the DNV rule and should be a separate stat, not an option on this one.