Add config-driven SAL precipitation verification - #221
Conversation
SAL (Structure–Amplitude–Location) precipitation scoring uses pysteps.verification.salscores; scikit-image is its transitive requirement. uv.lock is gitignored, so no lockfile update.
Framework-agnostic helpers for Structure–Amplitude–Location scoring: build a near-isotropic regular lat–lon raster, nearest-neighbour remap native (triangular ICON / KENDA) fields onto it via verification.spatial.spherical_nearest_neighbor_indices, and a thin wrapper over pysteps' sal() that gates dry windows. S/A/L are normalized ratios, hence invariant to a constant precipitation rescaling; the raster is required because the Location term assumes square pixels.
New SalConfig (enabled/params/leadtimes, detection thresholds, raster extent+spacing) under experiment.sal, mirroring ScoreMapsConfig, plus a validate_sal_leadtimes model validator. Regenerate config.schema.json.
verification_sal.py computes per-init SAL over all init times for a given (participant, param, lead time), remapping forecast and truth onto a common raster, and writes a CSV (one row per init: S,A,L + domain means, with a commented metadata header). New verification_sal / verification_sal_baseline rules produce these CSVs for every run and baseline; experiment_all requests them when experiment.sal.enabled (compute-only — plotting lives in the publication-figures workflow). Add an example sal block to varda-single-1.0.yaml (disabled; requires gridded truth).
Cover the raster builder, nearest-neighbour remap, and the SAL wrapper (identical-field zero, dry-window NaN, amplitude sign, location sign, unit-invariance), plus SalConfig parsing/validation and the leadtime validator. 12 tests.
config.schema.json was committed in 57a6b26 without the regenerated SAL entries, so config.py (with SalConfig) and the checked-in schema were out of sync and the pydantic-schema pre-commit hook failed in CI. Regenerate from config.py to add the SalConfig definition and the experiment.sal field.
Add a SAL block for all example configs, and adjust the one for varda-single-1.0.yaml such that it works with INCA.
The SAL rules listed only sal.py + data_input as code inputs, so edits to verification/spatial.py (the remap kernel, imported via verification.sal) or verification/__init__.py would not retrigger SAL, risking stale results. Add both as inputs. Also drop the baseline rule cpus_per_task from 24 (copied from scoremaps) to 2 for the single-threaded KDTree+pysteps job.
Use type(item).model_fields instead of the deprecated instance attribute in the SAL and scoremaps leadtime validators (Pydantic 2.11). Extend the grid_extent validator to reject non-finite values and implausible lat/lon ranges, not just wrong length/ordering.
SAL needs a resolved field; remapping sparse station truth (jretrieve, ~150 points) onto the ~1 km raster yields meaningless scores. Add a point-density guard (point_density_per_km2 + MIN_TRUTH_POINT_DENSITY) that fails fast on the first init. Also guard compute_sal against all-NaN input, drop the unused SEASONS constant, and add unit tests for the density guard and grid_extent.
SAL is defined for precipitation only, but sal.params accepted any parameter and the driver would silently compute meaningless scores on e.g. T_2M. Add a lenient guard requiring params to start with TOT_PREC (period-accumulated like TOT_PREC6 and bare cumulative TOT_PREC both allowed): a SalConfig field validator that fails fast at config load, plus a matching check in verification_sal.py for direct invocations. Add unit tests. If SAL should be deployed later on for different variables as well (which is conceivable but unlikely), this guard would need to be relaxed again.
Extract the shared steps-vs-requested-leadtimes check from the scoremaps and SAL validators into a single _reject_unproducible_leadtimes() helper, parameterized by the config-block label used in the error message. Both validators become thin enabled-guards, so future features validating a leadtimes field reuse one implementation instead of copy-pasting.
They are redundant information, the init hour is given in the first column.
The leadtime-validation dedup (ab3c497) hand-wrapped the _reject_unproducible_leadtimes signature across three lines; ruff collapses it to one. Without this the ruff-format CI hook fails.
|
I created this PR to integrate core SAL functionality into evalML main. Would be great if one of you could have a look at this (probably not both needed, please decide yourself who is less swamped / more motivated ;-)). The plotting functionality for the paper that builds on this will follow later on in a PR on top of |
Review feedback on #221 flagged the SAL branch as too verbose. The non-gridded-truth guard was a heavy offender: a point_density_per_km2 helper (equirectangular area approximation, cos-lat correction, degenerate -input handling), a MIN_TRUTH_POINT_DENSITY constant with a paragraph of justification, a ~20-line guard block, and three unit tests. The two truth regimes — a gridded analysis (millions of points over the domain) and a station network (~150) — are orders of magnitude apart, so the density metric bought precision the discriminator never needed. Replace the whole thing with a raw truth_lat.size check against MIN_TRUTH_POINTS (10k) that logs a warning instead of raising, letting the user proceed at their own risk. Drop the helper, the density tests, and the now-unused imports. Net ~80 lines removed. compute_sal's all-NaN guard and the grid_extent tidy from the original commit are unaffected.
Reduce PR #221 verbosity (reviewer flagged it as too much code). - Remove the unused --steps argparse arg (never read; neither SAL Snakemake rule passes it). - Remove the try/except that re-wrapped load_forecast_data errors as a verbose RuntimeError; the loader already raises a clear error, so a missing forecast still fails hard. - Remove the runtime precip-param check, a duplicate of SalConfig.validate_params_are_precip which runs at config-validation time before the script is invoked. No change to computed SAL scores. Net ~23 lines removed.
Reduce PR #221 verbosity (reviewer flagged it as too much code). - Remove the two property-sign tests (amplitude-positive-when- overforecast, location-positive-when-displaced): they assert pysteps' internal S/A/L math rather than the evalml wrapper. - Remove test_sal_config_accepts_precip_params: redundant with the reject test plus the default-SalConfig assertion. - Remove test_sal_disabled_skips_leadtime_validation: trivial enabled=False short-circuit. - Fold the standalone NaN-fill test into the nearest-neighbour test so remap_field's NaN->0 coverage is retained. Core guards kept (grid build, remap, identical=0, dry=NaN, rescale invariance). 9 tests pass. Net ~29 lines removed.
Reduce PR #221 verbosity (reviewer flagged it as too much code). The branch split the original single scoremaps leadtime validator into a _reject_unproducible_leadtimes helper plus two near-identical @model_validator methods (scoremaps, sal). Collapse them back into one validator that loops over both blocks. Error strings and behavior are unchanged (hard-fail if any run cannot produce a requested lead time).
Reduce PR #221 verbosity (reviewer flagged it as too much code). The six tracked example configs carried a full sal: block (params + leadtimes) that only restated the SalConfig defaults. Reduce each to `sal:\n enabled: false` — an absent-or-disabled block behaves identically. No behavior change.
Reduce PR #221 verbosity (reviewer flagged it as too much code). thr_factor and thr_quantile were exposed as config fields, CLI args, and Snakemake plumbing but were never overridden — no shipped config sets them. Remove that surface and let compute_sal apply the pysteps defaults (DEFAULT_THR_FACTOR/DEFAULT_THR_QUANTILE in verification.sal) directly. - config.py: drop the two SalConfig fields. - verification_sal.py: drop the --thr-factor/--thr-quantile args, the compute_sal kwargs, the now-unused imports, and the CSV header line. - verification.smk: drop the two flags from _SAL_ARGS. - test_sal.py: drop the thr_factor default assertion + unused import. - config.schema.json: regenerated. SAL scores are byte-identical. The grid_* regrid knobs are intentionally kept for now. Net ~37 hand-written lines removed (plus schema).
Reduce PR #221 verbosity (reviewer flagged it as too much code). grid_extent, grid_step_lat and grid_step_lon were exposed as config fields, CLI args and Snakemake plumbing but were never overridden. The SAL raster is fixed to the greater-Alpine ~1.1 km near-square grid, so move the three values to module constants in verification.sal (next to build_regular_grid, the code that uses them) and drop the knobs. - sal.py: add DEFAULT_GRID_EXTENT / DEFAULT_GRID_STEP_LAT / _LON. - config.py: drop the three SalConfig fields and the validate_grid_extent validator. - verification_sal.py: use the constants directly; drop the --grid-* args and the local constants. - verification.smk: drop the now-empty _SAL_ARGS plumbing. - test_sal.py: drop the grid_extent default assertion and the grid-extent validation test. - config.schema.json: regenerated. SAL scores are byte-identical. Changing the raster now means editing the three constants. Net ~90 hand-written lines removed (plus schema).
Reduce PR #221 verbosity (reviewer flagged it as too much code). iter_init_dirs was a verbatim copy of the helper in verification_scoremaps.py. Since the Snakefile always passes --reftimes, the run path can build each grib_dir directly from the reftime (run_root/<reftime>/[grib]) instead of scanning the run directory; a missing init still fails hard in load_forecast_data. --reftimes is now required for runs as well as baselines.
Boil the SAL driver down from 329 to 236 lines with no change to the CLI contract or the output CSV. Concision, no behaviour change: - _native_1d: replace the size-1-dim name loop with squeeze(drop=True) - derive the run GRIB directory inside the reftime loop instead of precomputing (reftime, grib_dir) tuples - drop the isfinite/nan formatting conditionals in the per-init log line; %.3f already renders nan - pd.DataFrame(rows) without the redundant columns= list - np.allclose for the shared-grid test, one log line instead of three - --reftimes required=True instead of a manual parser.error, which also makes the "no inits processed" guard unreachable - condense the argparse help strings Dropped: - the truth valid-time pre-flight. Absent dates are already excluded via the blacklist: config knob, and a genuine gap still fails hard in data_input, so the check only bought an earlier, prettier error. - the INCA accumulation guard. It rejected TOT_PREC6 as unsupported, but the INCA reader sums native 10-minute slots and caps steps at 6, so TOT_PREC6 at step 6 runs; lead times beyond 6h already fail upstream with a clear message. Also fixes the RUF059 and FURB122 lint findings the file carried.
Reduce PR #221 verbosity further (reviewer flagged it as too much code). Net -112 lines across 5 files with no change to the CLI contract and no change to any computed value: sal.py 106 -> 87, verification_sal.py 236 -> 184, test_sal.py 107 -> 73. Concision, no behaviour change: - shorten both module docstrings - one log line per init instead of two, and drop the startup banner, the raster line and the remap-indices line; the log filename already encodes run/truth/param/leadtime - accumulate rows as tuples, so pd.DataFrame needs columns= again (re-adding what c9474a8 dropped, since rows are no longer dicts) - condense the argparse help strings a second time - drop the comment banner above the SAL rules; no other rule section in verification.smk has one - fix the experiment_all comment: SAL writes CSVs, not NetCDFs sal.py surface: - replace build_regular_grid(extent, step_lat, step_lon) with a no-arg sal_raster() returning (lat2d, lon2d). The raster has been hardcoded since 076a2b4, so the parameters were exercised only by the unit test; the two step values are now local and GRID_EXTENT is the one constant the driver still needs. - compute_sal: nan_to_num both fields and gate on max() > 0 instead of filtering to the finite subset. NaNs now count as dry, which is what the driver already guaranteed by filling NaN -> 0 in remap_field, so the pipeline path is bit-identical. Retires the "All-NaN slice" RuntimeWarning the public helper could emit. Tests dropped: - the unit-invariance test: that is a property of pysteps' SAL, not of this wrapper - the SalConfig defaults/extra-forbid test: it asserted pydantic does what pydantic does - the nearest-neighbour half of the remap test, already covered by test_spherical_nearest_neighbor_indices_returns_expected_points in test_spatial_mapping.py; the NaN -> 0 assertion is kept The CSV metadata header now reports grid_cells: 851x1311 in place of grid_step: (0.01, 0.0145) -- same information, read off the raster that was actually built. Header text only; all value columns are unchanged.
|
Hi @clairemerker @frazane, I boiled down this PR as much as I could / dared (because feedback last time was that the PR was too verbose). Would be great if one of you could have a look at it! |
clairemerker
left a comment
There was a problem hiding this comment.
Thanks, I think it looks fine :)
I added some comments to the changes, and I wonder if you should already include an integration test config for it?
|
Thanks for the feedback @clairemerker, I will look into it. Yes, good idea, I will include a longtest as well. |
Co-authored-by: Claire Merker <34312518+clairemerker@users.noreply.github.com>
Co-authored-by: Claire Merker <34312518+clairemerker@users.noreply.github.com>
Co-authored-by: Claire Merker <34312518+clairemerker@users.noreply.github.com>
Co-authored-by: Claire Merker <34312518+clairemerker@users.noreply.github.com>
Co-authored-by: Claire Merker <34312518+clairemerker@users.noreply.github.com>
Co-authored-by: Claire Merker <34312518+clairemerker@users.noreply.github.com>
Co-authored-by: Claire Merker <34312518+clairemerker@users.noreply.github.com>
Co-authored-by: Claire Merker <34312518+clairemerker@users.noreply.github.com>
The sparse-truth warning fired too late and too weakly: load_truth_data() also has a jretrievedwh branch returning ~150 SwissMetNet stations, and the script only used the .zarr suffix to decide whether to pre-open the store lazily, silently falling back to lazy_ds=None for anything else. Require a .zarr truth root at the top of main() instead, before any data is read. load_truth_data() supports only .zarr and jretrieve, so this rejects station observations and every other unsupported root in one check, and the lazy pre-open becomes unconditional. The MIN_TRUTH_POINTS heuristic goes away with it — its threshold depended on truth resolution and domain, and a warning could not prevent the meaningless scores it warned about.
Adds a baseline-only integration config (ICON-CH2-CTRL against the KENDA-CH1 zarr, two inits, TOT_PREC6 at +6/+12 h) and a longtest driving `evalml experiment` over it, exercising the full SAL chain end to end without a GPU: schema, SalConfig validators, experiment_all target expansion, verification_sal_baseline, verification_sal.py and verification.sal. The test needs no blessed reference files. It leans on internal redundancy instead: A is recomputed from the row's own mean columns (exact), and truth_mean is cross-checked across the two output files that reach the same accumulation window via different (reftime, leadtime) pairs -- with rtol 1e-9, since agreement is close but not bit-exact given the zarr accumulates chunk-wise -- plus component bounds and an all-NaN guard. S and L are deliberately not pinned: they pass through a discrete object segmentation, so a tight tolerance would be flaky and a loose one would catch little. The config writes to a dedicated output_root so the exactly-one-CSV assertion cannot trip over an unrelated experiment sharing output/. The two unit tests cover the script's pre-flight guards (lead time below the accumulation period, non-zarr truth); both fire before any I/O, so they run on GitHub Actions. No CI changes: the CSCS longtest pipeline runs `pytest tests/integration -m longtest` over the whole directory, so this activates when that job lands with #226.
|
Added a longtest in 74b3be2. |
Answers the review question "which domain does GRID_EXTENT correspond to?": it corresponded to nothing named. 076a2b4 had hardcoded the raster as module constants after the grid knobs were flagged as unused verbosity, which left the domain both unexplained and unadjustable. Reinstate it as configuration and give it a default that can be justified. The default extent is the bounding box of the ICON-CH1 analysis, i.e. the SAL truth, and also what the raster was before 076a2b4: (-1, 18, 42, 50.5) at 0.01 lat x 0.0145 lon, 851x1311 cells. KENDA-CH1 spans lon -0.8171..17.7106, lat 42.0279..50.5005, so the default rounds that outwards by <=0.29 deg. It is also the raster scripts/sal_per_init.py uses for the paper case studies, which keeps those numbers comparable. The icon-ch plotting domain (0, 17.5, 40.5, 53) was tried first and rejected: it is an animation display window that reaches 2.5 deg north and 1.5 deg south of the analysis, and the remap applies no distance cutoff, so 37% of its cells would merely repeat the nearest border value. The ICON-CH1 box leaves 14%, and cropping tighter buys little -- ICON-CH1 is a rotated-pole patch whose corners fall inside its lat/lon box. Across four real precipitation cases the larger raster left S and A within 0.01 but deflated L by 10-17%, through the longer domain diagonal that pysteps normalises distances by. - sal.py: sal_raster(extent, step_lat, step_lon) plus DEFAULT_GRID_EXTENT / DEFAULT_GRID_STEP_LAT / DEFAULT_GRID_STEP_LON. - config.py: SalConfig.grid_extent with its validator, and grid_step_lat / grid_step_lon constrained to > 0, since a non-positive step yields no raster. - verification_sal.py: --grid-extent / --grid-step-lat / --grid-step-lon, echoed to the log and into the CSV metadata header. - verification.smk: _SAL_ARGS plumbing, taken from the validated config. - test_sal.py: default, validation and custom-extent cases; the config defaults are pinned to the module constants so the two cannot drift apart. - config.schema.json: regenerated. Also replaces the raster's half-step upper-bound tolerance with a relative epsilon. An extent whose span is not a whole multiple of the spacing used to overshoot by up to half a cell, so the raster could reach outside the extent it was asked for; it now stops short of the bound instead. No effect on the default extent, whose raster is unchanged. Example configs are untouched: omitting the keys yields exactly these defaults.
Adds SAL (Structure–Amplitude–Location; Wernli et al. 2008) as a config-driven precipitation verification capability, mirroring how scorecards/scoremaps work: enabled via an
experiment.salblock and computed for every requested run and baseline. Compute-only — plotting lives in the publication-figures workflow.What
src/verification/sal.py— pystepssal()wrapper + near-isotropic raster builder + nearest-neighbour remap (reusesverification.spatial). Dry windows (either field everywhere ≤ 0, NaNs counting as dry) return(NaN, NaN, NaN)instead of raising. Detection thresholds are the pysteps defaults (thr_factor=0.067,thr_quantile=0.95) as module constants — not config knobs.workflow/scripts/verification_sal.py— per-init compute driver (one job per participant/param/lead time). Scores the init times given by--reftimes, remaps both forecast and truth onto a common raster, and writes a CSV (one row per init:reftime,S,A,L,fcst_mean,truth_mean+ a commented metadata header). Dry-window rows are retained with NaN S/A/L so a downstream wet-case filter can drop them. Every configured init must be present across forecast and truth — a missing one is a hard error, never a silent skip, so run and baseline are scored over an identical sample. SAL is a per-case scalar score, so a table is the natural container.SalConfigunderexperiment.sal—enabled,params,leadtimes, and the scoring raster (grid_extent,grid_step_lat,grid_step_lon). Validators reject non-precipitation params and malformedgrid_extent; schema regenerated. The raster defaults to the ICON-CH1 analysis bounding box ([-1.0, 18.0, 42.0, 50.5]) at0.01°lat ×0.0145°lon, i.e. ~1.1 km cells that are metrically square at 46.4°N (pysteps' Location term assumes square pixels).ConfigModel.validate_leadtimes_produciblemodel validator, which rejects lead times not produced by every participant for whichever of the two blocks is enabled.verification_sal/verification_sal_baselineSnakemake rules (output keyed by truth hash) + gated fan-out inexperiment_all.sal:block added to all six example configs asenabled: false; everything else comes from theSalConfigdefaults (TOT_PREC6at[6, 12, 18, 24, 30]h).pyproject.toml,uv.lock).Example output
The CSVs are the only data product: one per (participant, param, lead time), at
data/runs/<run_id>/sal/<param>_<leadtime>_<truth_hash>.csv(anddata/baselines/<baseline_id>/sal/…). Below issal/TOT_PREC6_6_2b83.csvfrom a five-init run —TOT_PREC6at+6 h, scored against the KENDA-CH1 analysis zarr:(Values rounded for display; the CSV carries full float64 precision.)
Reading it: one row per initialisation,
S/A/Ldimensionless and signed,fcst_mean/truth_meanthe domain-mean 6 h accumulation in mm over the scoring raster. NegativeAthroughout is a dry bias — this forecast under-predicts domain-total precipitation at every init here. The commented header pins everything needed to reproduce the numbers: the param and its accumulation period, the lead time, the baseline member, the raster extent and cell count, and the source archive.Testing
tests/unit/test_sal.py): raster shape/bounds/orientation and custom extents, the remap gather and its NaN→0 fill,compute_salon identical fields (all-zero triple) and on a dry window (all-NaN triple),SalConfiggrid defaults matching the module constants, the grid and non-precip-param validators, the lead-time producibility validator, and the two script guards (lead time below the accumulation period, non-zarr truth).tests/integration/test_sal_small.py+configs/sal_small.yaml): runsevalml experimentend to end on a minimal baseline-only config (ICON-CH2-CTRL, no inference — so no GPU, MLflow or DWH, only/store_newaccess) and checks the scores are consistent, not merely present — the documented metadata header, column set, raster size and one row per configured init; finite non-negative means;Aequal to the identity recomputed from the row's own mean columns; S/A/L within their Wernli et al. bounds; S/L not NaN for every init; and — cross-checked across two independent output files — that the truth column depends only on valid time, which pins thereftime + step→ accumulation-window mapping. All reference-free, so there is no blessed output to maintain. Markedlongtest, so it is skipped by default and on GitHub Actions and runs on the CSCS balfrin runner.Notes
grid_extentclose to the truth's footprint. The remap has no distance cutoff, so cells beyond it merely repeat the nearest border value, and a larger raster also deflates pysteps' L term via its longer diagonal.0/6/1) restrictssal.leadtimesfor all participants — use a config without such a baseline for longer lead times.jretrievedwh:…) included — is rejected up front with a clear error, since remapping ~150 points onto the raster would yield meaningless scores.verification_sal_baseline. The run path (verification_sal) depends oninference_execute, so exercising it needs a GPU or the inference-replay fixture; the only untested delta is therun_root/<reftime>/[grib]resolution inverification_sal.py.