From de78ecf7f9a32b797989c8eca362988217646c11 Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Fri, 9 Oct 2026 11:46:46 +0200 Subject: [PATCH 01/15] Add side-by-side map comparisons of runs and baselines to the showcase --- src/evalml/config.py | 40 ++++++++ src/plotting/styles.py | 86 +++++++++++++++++ workflow/Snakefile | 17 ++++ workflow/rules/common.smk | 24 +++-- workflow/rules/plot.smk | 83 ++++++++++++++++ workflow/scripts/plot_compare_maps.py | 121 ++++++++++++++++++++++++ workflow/scripts/plot_forecast_frame.py | 85 +---------------- workflow/tools/config.schema.json | 76 +++++++++++++++ 8 files changed, 443 insertions(+), 89 deletions(-) create mode 100644 src/plotting/styles.py create mode 100644 workflow/scripts/plot_compare_maps.py diff --git a/src/evalml/config.py b/src/evalml/config.py index 3980046a..5b9c965c 100644 --- a/src/evalml/config.py +++ b/src/evalml/config.py @@ -337,6 +337,42 @@ class AnimationsConfig(BaseModel): ) +class ComparisonsConfig(BaseModel): + """Configuration for side-by-side map comparisons of runs and baselines.""" + + enabled: bool = Field( + default=False, + description="Whether to generate side-by-side maps of each run and the selected baselines.", + ) + baselines: List[str] | None = Field( + default=None, + description=( + "Labels of the baselines to show next to each run. " + "None (default) means all baselines in the config." + ), + ) + domains: List[str | DomainConfig] = Field( + default=["alps"], + description=( + "Domains to plot. Same format as `animations.domains`: a named domain " + "defined in plotting.DOMAINS or a custom domain dict." + ), + ) + params: List[str] = Field( + default=["SP_10M"], + description="Parameters to plot.", + ) + lead_times: str | None = Field( + default=None, + description=( + "Lead times as start/stop/step (hours). None (default) means all lead " + "times produced by both the run and every selected baseline." + ), + ) + + model_config = {"extra": "forbid"} + + class ScorecardConfig(BaseModel): """Configuration for a single named scorecard.""" @@ -395,6 +431,10 @@ class ShowcaseConfig(BaseModel): default_factory=AnimationsConfig, description="Configuration for animation generation.", ) + comparisons: ComparisonsConfig = Field( + default_factory=ComparisonsConfig, + description="Configuration for side-by-side map comparisons.", + ) class Locations(BaseModel): diff --git a/src/plotting/styles.py b/src/plotting/styles.py new file mode 100644 index 00000000..d1bbb80f --- /dev/null +++ b/src/plotting/styles.py @@ -0,0 +1,86 @@ +"""Shared plotting styles and field preprocessing for map plots.""" + +from earthkit.meteo.utils.convert import kelvin_to_celsius +import earthkit.meteo.wind as ekm_wind +import earthkit.plots as ekp +from matplotlib.colors import Colormap +import numpy as np + +from plotting.colormap_defaults import CMAP_DEFAULTS + + +def get_style(param, units_override=None, accu=1): + """Get style and colormap settings for the plot.""" + lookup = f"{param}_{accu}H" if param == "TOT_PREC" else param + cfg = CMAP_DEFAULTS[lookup] + units = units_override if units_override is not None else cfg.get("units", "") + + bounds = cfg.get("bounds", cfg.get("levels", None)) + prebuilt_cmap = cfg.get("cmap", None) + + # When the config provides a pre-built matplotlib Colormap (e.g. a + # ListedColormap), we must use the earthkit Style's vmin/vmax path, which + # handles isinstance(colors, Colormap) correctly. The levels path calls + # cmap_and_norm → len(Colormap) → TypeError. + # + # earthkit's configure_style intercepts _STYLE_KWARGS (cmap/colors/levels/vmin/vmax) + # from tricontourf kwargs before matplotlib sees them. To work around this: + # - embed the Colormap directly in the Style via colors= + # - inject bounds as 'levels' into style._kwargs so they survive to matplotlib + # - pass only norm= as a kwarg (not in _STYLE_KWARGS, so not intercepted) + extend = cfg.get("extend", "both") + + if isinstance(prebuilt_cmap, Colormap) and bounds is not None: + style = ekp.styles.Style( + colors=prebuilt_cmap, + vmin=bounds[0], + vmax=bounds[-1], + extend=extend, + units=units, + ) + style._kwargs["levels"] = list(bounds) + return { + "style": style, + "norm": cfg.get("norm", None), + } + + return { + "style": ekp.styles.Style( + levels=bounds, + extend=extend, + units=units, + colors=cfg.get("colors", None), + ), + "norm": cfg.get("norm", None), + "cmap": prebuilt_cmap, + "vmin": cfg.get("vmin", None), + "vmax": cfg.get("vmax", None), + } + + +def preprocess_field(param: str, state: dict): + """ + - Temperatures: K -> °C + - Wind speed: sqrt(u^2 + v^2) + - Precipitation: m -> mm + Returns: (field_array, units_override or None) + """ + fields = state["fields"] + if param in ("T_2M", "TD_2M", "T", "TD"): + return kelvin_to_celsius(fields[param]), "°C" + if param == "SP_10M": + return ekm_wind.speed(fields["U_10M"], fields["V_10M"]), "m/s" + if param == "SP": + return ekm_wind.speed(fields["U"], fields["V"]), "m/s" + if param == "TOT_PREC": + return np.maximum(fields[param], 0), "mm" + if param in ("CLCT", "CLCL"): + # Avoid exact 0/1 plateaus breaking tricontourf on orthographic + # projections (tmp/reproduce_clct_bug.py). Pair with extend="neither". + # Any new bounded field with silent-blank or GeometryCollection-crash + # globe frames likely needs the same clip-away-from-boundary fix. + return np.clip(fields[param], 1e-6, 1 - 1e-6), None + if param == "SSRD": + # Same issue, bottom boundary only (night-side plateau). + return np.maximum(fields[param], 1e-6), None + return fields[param], None diff --git a/workflow/Snakefile b/workflow/Snakefile index 67838e97..68cc8ea6 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -224,6 +224,23 @@ rule showcase_all: if config["showcase"]["meteograms"]["enabled"] else [] ), + ( + [ + path + for run_id in CANDIDATES + for param in COMPARISONS_CONFIG["params"] + for path in expand( + rules.plot_compare_maps.output, + init_time=[t.strftime("%Y%m%d%H%M") for t in REFTIMES], + run_id=run_id, + param=param, + leadtime=get_compare_leadtimes(run_id, param), + showcase=EXPERIMENT_NAME, + ) + ] + if COMPARISONS_CONFIG["enabled"] + else [] + ), rule sandbox_all: diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index acf76746..c91c669c 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -116,9 +116,11 @@ def parse_regions(): return json.dumps(result) -def parse_showcase_regions(): +def parse_showcase_regions(domains=None): """Parse showcase domains from config. + ``domains`` defaults to ``showcase.animations.domains``. + Returns a dict mapping domain name -> {extent, projection, rotate, hours_per_revolution}. Named domains (strings) have extent=None and projection=None, meaning the plot script will fall back to the DOMAINS lookup. @@ -126,12 +128,14 @@ def parse_showcase_regions(): rotate/hours_per_revolution only take effect when extent is None (full-globe domains). """ + if domains is None: + domains = ( + config.get("showcase", {}) + .get("animations", {}) + .get("domains", ["globe", "europe", "switzerland"]) + ) result = {} - for r in ( - config.get("showcase", {}) - .get("animations", {}) - .get("domains", ["globe", "europe", "switzerland"]) - ): + for r in domains: if isinstance(r, str): result[r] = { "extent": None, @@ -436,6 +440,14 @@ SHOWCASE_CONFIG = { "params": _showcase.get("params", ["T_2M", "SP_10M"]), "fps": _showcase.get("animations", {}).get("frames_per_second", 2.0), } +_comparisons = _showcase.get("comparisons", {}) +COMPARISONS_CONFIG = { + "enabled": _comparisons.get("enabled", False), + "regions": parse_showcase_regions(_comparisons.get("domains", ["alps"])), + "params": _comparisons.get("params", ["SP_10M"]), + "baselines": _comparisons.get("baselines"), + "lead_times": _comparisons.get("lead_times"), +} EXPERIMENT_PARAMS = config.get("experiment", {}).get( "params", ["T_2M", "TD_2M", "SP_10M", "PS", "PMSL", "TOT_PREC6"] ) diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 5afdefd8..c6536a37 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -214,3 +214,86 @@ use rule plot_scoremaps as plot_scoremaps_baseline with: log: OUT_ROOT / "logs/plot_scoremaps/{experiment}/{baseline_id}-{param}-{score}-{region}-{season}-{init_hour}-{leadtime}.log", + + +def _get_compare_baselines(wc) -> list[dict[str, str]]: + """Baselines selected for the side-by-side comparison maps.""" + selected = COMPARISONS_CONFIG["baselines"] + baselines = _get_available_baselines(wc) + if selected is None: + return baselines + return [b for b in baselines if b["label"] in selected] + + +def get_compare_leadtimes(run_id, param): + """Lead times produced by the run and by every selected baseline.""" + requested = COMPARISONS_CONFIG["lead_times"] + requested = ( + "all" if requested is None else resolve_leadtimes(requested, param=param) + ) + leadtimes = set(resolve_leadtimes(RUN_CONFIGS[run_id]["steps"], requested, param)) + selected = COMPARISONS_CONFIG["baselines"] + for cfg in BASELINE_CONFIGS.values(): + if selected is None or cfg.get("label") in selected: + leadtimes &= set(resolve_leadtimes(cfg["steps"], requested, param)) + return sorted(leadtimes) + + +rule plot_compare_maps: + input: + script="workflow/scripts/plot_compare_maps.py", + inference_okfile=rules.inference_execute.output.okfile, + output: + expand( + OUT_ROOT + / "results/{{showcase}}/{{run_id}}/{{init_time}}/compare/{{init_time}}_{{leadtime}}_{{param}}_{region}.png", + region=list(COMPARISONS_CONFIG["regions"].keys()), + ), + log: + OUT_ROOT + / "logs/{showcase}/{run_id}/{init_time}/plot_compare_maps_{leadtime}_{param}.log", + wildcard_constraints: + leadtime=r"\d+", + resources: + slurm_partition="postproc", + cpus_per_task=4, + runtime="30m", + params: + fcst_grib=lambda wc: ( + Path(OUT_ROOT) / f"data/runs/{wc.run_id}/{wc.init_time}/grib" + ).resolve(), + fcst_label=lambda wc: RUN_CONFIGS[wc.run_id]["label"], + baseline_roots=lambda wc: [x["root"] for x in _get_compare_baselines(wc)], + baseline_labels=lambda wc: [x["label"] for x in _get_compare_baselines(wc)], + regions_json=json.dumps(COMPARISONS_CONFIG["regions"]), + outdir=lambda wc: str( + ( + Path(OUT_ROOT) + / f"results/{wc.showcase}/{wc.run_id}/{wc.init_time}/compare" + ).resolve() + ), + shell: + """ + set -euo pipefail + export ECCODES_DEFINITION_PATH=$(realpath .venv/share/eccodes-cosmo-resources/definitions) + + BASELINE_ROOTS=({params.baseline_roots:q}) + BASELINE_LABELS=({params.baseline_labels:q}) + + CMD_ARGS=( + --forecast {params.fcst_grib:q} + --forecast_label {params.fcst_label:q} + --date {wildcards.init_time:q} + --leadtime {wildcards.leadtime:q} + --param {wildcards.param:q} + --regions_json {params.regions_json:q} + --outdir {params.outdir:q} + ) + + for i in "${{!BASELINE_ROOTS[@]}}"; do + CMD_ARGS+=(--baseline "${{BASELINE_ROOTS[$i]}}") + CMD_ARGS+=(--baseline_label "${{BASELINE_LABELS[$i]}}") + done + + uv run python {input.script} "${{CMD_ARGS[@]}}" >{log} 2>&1 + """ diff --git a/workflow/scripts/plot_compare_maps.py b/workflow/scripts/plot_compare_maps.py new file mode 100644 index 00000000..72f30074 --- /dev/null +++ b/workflow/scripts/plot_compare_maps.py @@ -0,0 +1,121 @@ +"""Plot side-by-side maps of an ML run and baselines at a common valid time.""" + +import json +import logging +from argparse import ArgumentParser +from datetime import datetime +from datetime import timedelta +from pathlib import Path + +import earthkit.plots as ekp +from earthkit.meteo.utils.convert import kelvin_to_celsius +import numpy as np + +from data_input import load_forecast_data +from data_input import parse_aggregated_param +from plotting import DOMAINS +from plotting import get_projection +from plotting import StatePlotter +from plotting.styles import get_style + +LOG = logging.getLogger(__name__) +logging.basicConfig( + level=logging.INFO, + format="%(asctime)s - %(name)s - %(levelname)s - %(message)s", +) + + +def load_source(root, init_time, lead_time, param): + """Return (lon, lat, values) as flat arrays.""" + ds = load_forecast_data(root, init_time, [lead_time], [param]).squeeze() + values = np.asarray(ds[param].values).ravel() + if param in ("T_2M", "TD_2M"): + values = kelvin_to_celsius(values) + lon = np.asarray(ds["longitude"].values).ravel() + lat = np.asarray(ds["latitude"].values).ravel() + return lon, lat, values + + +def resolve_region(region_name, region_cfg): + """Return (projection, extent) for a showcase region.""" + if region_cfg.get("extent") is not None: + projection = get_projection(region_cfg.get("projection") or "orthographic") + return projection, region_cfg["extent"] + return DOMAINS[region_name]["projection"], DOMAINS[region_name]["extent"] + + +def main(): + parser = ArgumentParser() + parser.add_argument("--forecast", type=str, help="Directory to ML run grib data") + parser.add_argument("--forecast_label", type=str, default="Forecast") + parser.add_argument( + "--baseline", type=str, action="append", default=[], help="Baseline root" + ) + parser.add_argument( + "--baseline_label", type=str, action="append", default=[], help="Label" + ) + parser.add_argument("--date", type=str, help="reference datetime") + parser.add_argument("--leadtime", type=int, help="leadtime") + parser.add_argument("--param", type=str, help="parameter") + parser.add_argument( + "--regions_json", + type=str, + help="JSON dict mapping region name -> {extent, projection}", + ) + parser.add_argument("--outdir", type=str, help="output directory") + args = parser.parse_args() + + init_time = datetime.strptime(args.date, "%Y%m%d%H%M") + lead_time = args.leadtime + param = args.param + regions = json.loads(args.regions_json) + outdir = Path(args.outdir) + outdir.mkdir(parents=True, exist_ok=True) + valid_time = init_time + timedelta(hours=lead_time) + + sources = [(args.forecast_label, args.forecast)] + list( + zip(args.baseline_label, args.baseline) + ) + data = {} + for label, root in sources: + LOG.info("Loading '%s' from %s", label, root) + data[label] = load_source(root, init_time, lead_time, param) + + base_param, accu = parse_aggregated_param(param) + style = get_style(base_param, accu=accu or 1) + ncols = len(sources) + + for region_name, region_cfg in regions.items(): + LOG.info("Plotting region %s", region_name) + projection, extent = resolve_region(region_name, region_cfg) + fig = ekp.Figure( + crs=projection, + domain=extent, + rows=1, + columns=ncols, + size=(5 * ncols, 5), + ) + for col, (label, _) in enumerate(sources): + subplot = fig.add_map(row=0, column=col) + lon, lat, values = data[label] + plotter = StatePlotter(lon, lat, outdir) + plotter.plot_field( + subplot, + values, + title=label, + colorbar=False, + **style, + ) + # One colorbar shared by all panels, which use identical levels. + fig.legend(location="bottom") + fig.title( + f"{param} valid: {valid_time:%Y-%m-%d %H:%M} UTC " + f"(init: {init_time:%Y-%m-%d %H:%M} UTC, +{lead_time}h)" + ) + outfn = outdir / f"{args.date}_{lead_time}_{param}_{region_name}.png" + fig.save(outfn, dpi=200) + LOG.info("saved: %s", outfn) + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/plot_forecast_frame.py b/workflow/scripts/plot_forecast_frame.py index ee8e564e..880ec56d 100644 --- a/workflow/scripts/plot_forecast_frame.py +++ b/workflow/scripts/plot_forecast_frame.py @@ -4,17 +4,13 @@ from pathlib import Path import cartopy.crs as ccrs -from earthkit.meteo.utils.convert import kelvin_to_celsius -import earthkit.meteo.wind as ekm_wind -import earthkit.plots as ekp -from matplotlib.colors import Colormap -import numpy as np from plotting import DOMAINS from plotting import get_projection from plotting import StatePlotter -from plotting.colormap_defaults import CMAP_DEFAULTS from plotting.compat import load_state_from_grib +from plotting.styles import get_style +from plotting.styles import preprocess_field LOG = logging.getLogger(__name__) logging.basicConfig( @@ -23,83 +19,6 @@ ) -def get_style(param, units_override=None, accu=1): - """Get style and colormap settings for the plot.""" - lookup = f"{param}_{accu}H" if param == "TOT_PREC" else param - cfg = CMAP_DEFAULTS[lookup] - units = units_override if units_override is not None else cfg.get("units", "") - - bounds = cfg.get("bounds", cfg.get("levels", None)) - prebuilt_cmap = cfg.get("cmap", None) - - # When the config provides a pre-built matplotlib Colormap (e.g. a - # ListedColormap), we must use the earthkit Style's vmin/vmax path, which - # handles isinstance(colors, Colormap) correctly. The levels path calls - # cmap_and_norm → len(Colormap) → TypeError. - # - # earthkit's configure_style intercepts _STYLE_KWARGS (cmap/colors/levels/vmin/vmax) - # from tricontourf kwargs before matplotlib sees them. To work around this: - # - embed the Colormap directly in the Style via colors= - # - inject bounds as 'levels' into style._kwargs so they survive to matplotlib - # - pass only norm= as a kwarg (not in _STYLE_KWARGS, so not intercepted) - extend = cfg.get("extend", "both") - - if isinstance(prebuilt_cmap, Colormap) and bounds is not None: - style = ekp.styles.Style( - colors=prebuilt_cmap, - vmin=bounds[0], - vmax=bounds[-1], - extend=extend, - units=units, - ) - style._kwargs["levels"] = list(bounds) - return { - "style": style, - "norm": cfg.get("norm", None), - } - - return { - "style": ekp.styles.Style( - levels=bounds, - extend=extend, - units=units, - colors=cfg.get("colors", None), - ), - "norm": cfg.get("norm", None), - "cmap": prebuilt_cmap, - "vmin": cfg.get("vmin", None), - "vmax": cfg.get("vmax", None), - } - - -def preprocess_field(param: str, state: dict): - """ - - Temperatures: K -> °C - - Wind speed: sqrt(u^2 + v^2) - - Precipitation: m -> mm - Returns: (field_array, units_override or None) - """ - fields = state["fields"] - if param in ("T_2M", "TD_2M", "T", "TD"): - return kelvin_to_celsius(fields[param]), "°C" - if param == "SP_10M": - return ekm_wind.speed(fields["U_10M"], fields["V_10M"]), "m/s" - if param == "SP": - return ekm_wind.speed(fields["U"], fields["V"]), "m/s" - if param == "TOT_PREC": - return np.maximum(fields[param], 0), "mm" - if param in ("CLCT", "CLCL"): - # Avoid exact 0/1 plateaus breaking tricontourf on orthographic - # projections (tmp/reproduce_clct_bug.py). Pair with extend="neither". - # Any new bounded field with silent-blank or GeometryCollection-crash - # globe frames likely needs the same clip-away-from-boundary fix. - return np.clip(fields[param], 1e-6, 1 - 1e-6), None - if param == "SSRD": - # Same issue, bottom boundary only (night-side plateau). - return np.maximum(fields[param], 1e-6), None - return fields[param], None - - def main(): parser = ArgumentParser() parser.add_argument( diff --git a/workflow/tools/config.schema.json b/workflow/tools/config.schema.json index ea2b84ba..b7f9c93f 100644 --- a/workflow/tools/config.schema.json +++ b/workflow/tools/config.schema.json @@ -88,6 +88,78 @@ "title": "BaselineItem", "type": "object" }, + "ComparisonsConfig": { + "additionalProperties": false, + "description": "Configuration for side-by-side map comparisons of runs and baselines.", + "properties": { + "enabled": { + "default": false, + "description": "Whether to generate side-by-side maps of each run and the selected baselines.", + "title": "Enabled", + "type": "boolean" + }, + "baselines": { + "anyOf": [ + { + "items": { + "type": "string" + }, + "type": "array" + }, + { + "type": "null" + } + ], + "default": null, + "description": "Labels of the baselines to show next to each run. None (default) means all baselines in the config.", + "title": "Baselines" + }, + "domains": { + "default": [ + "alps" + ], + "description": "Domains to plot. Same format as `animations.domains`: a named domain defined in plotting.DOMAINS or a custom domain dict.", + "items": { + "anyOf": [ + { + "type": "string" + }, + { + "$ref": "#/$defs/DomainConfig" + } + ] + }, + "title": "Domains", + "type": "array" + }, + "params": { + "default": [ + "SP_10M" + ], + "description": "Parameters to plot.", + "items": { + "type": "string" + }, + "title": "Params", + "type": "array" + }, + "lead_times": { + "anyOf": [ + { + "type": "string" + }, + { + "type": "null" + } + ], + "default": null, + "description": "Lead times as start/stop/step (hours). None (default) means all lead times produced by both the run and every selected baseline.", + "title": "Lead Times" + } + }, + "title": "ComparisonsConfig", + "type": "object" + }, "Dashboard": { "description": "Settings for the dashboard", "properties": { @@ -860,6 +932,10 @@ }, "animations": { "$ref": "#/$defs/AnimationsConfig" + }, + "comparisons": { + "$ref": "#/$defs/ComparisonsConfig", + "description": "Configuration for side-by-side map comparisons." } }, "title": "ShowcaseConfig", From 2801b9960460eaa8f499e6cb1e40d85ff7420308 Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Fri, 9 Oct 2026 11:52:26 +0200 Subject: [PATCH 02/15] Plot list of lead times --- src/evalml/config.py | 8 ++++---- workflow/rules/common.smk | 2 +- workflow/rules/plot.smk | 5 +---- workflow/tools/config.schema.json | 18 ++++++++---------- 4 files changed, 14 insertions(+), 19 deletions(-) diff --git a/src/evalml/config.py b/src/evalml/config.py index 5b9c965c..55e5b8cc 100644 --- a/src/evalml/config.py +++ b/src/evalml/config.py @@ -362,11 +362,11 @@ class ComparisonsConfig(BaseModel): default=["SP_10M"], description="Parameters to plot.", ) - lead_times: str | None = Field( - default=None, + lead_times: List[int] = Field( + default=[24], description=( - "Lead times as start/stop/step (hours). None (default) means all lead " - "times produced by both the run and every selected baseline." + "Lead times (hours) to plot. Lead times not produced by the run or by " + "a selected baseline are skipped." ), ) diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index c91c669c..8b2be13c 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -446,7 +446,7 @@ COMPARISONS_CONFIG = { "regions": parse_showcase_regions(_comparisons.get("domains", ["alps"])), "params": _comparisons.get("params", ["SP_10M"]), "baselines": _comparisons.get("baselines"), - "lead_times": _comparisons.get("lead_times"), + "lead_times": _comparisons.get("lead_times", [24]), } EXPERIMENT_PARAMS = config.get("experiment", {}).get( "params", ["T_2M", "TD_2M", "SP_10M", "PS", "PMSL", "TOT_PREC6"] diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index c6536a37..70296532 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -226,11 +226,8 @@ def _get_compare_baselines(wc) -> list[dict[str, str]]: def get_compare_leadtimes(run_id, param): - """Lead times produced by the run and by every selected baseline.""" + """Requested lead times produced by the run and by every selected baseline.""" requested = COMPARISONS_CONFIG["lead_times"] - requested = ( - "all" if requested is None else resolve_leadtimes(requested, param=param) - ) leadtimes = set(resolve_leadtimes(RUN_CONFIGS[run_id]["steps"], requested, param)) selected = COMPARISONS_CONFIG["baselines"] for cfg in BASELINE_CONFIGS.values(): diff --git a/workflow/tools/config.schema.json b/workflow/tools/config.schema.json index b7f9c93f..899d2a57 100644 --- a/workflow/tools/config.schema.json +++ b/workflow/tools/config.schema.json @@ -144,17 +144,15 @@ "type": "array" }, "lead_times": { - "anyOf": [ - { - "type": "string" - }, - { - "type": "null" - } + "default": [ + 24 ], - "default": null, - "description": "Lead times as start/stop/step (hours). None (default) means all lead times produced by both the run and every selected baseline.", - "title": "Lead Times" + "description": "Lead times (hours) to plot. Lead times not produced by the run or by a selected baseline are skipped.", + "items": { + "type": "integer" + }, + "title": "Lead Times", + "type": "array" } }, "title": "ComparisonsConfig", From 335cedcccc929f37aa6cf2eb952a2a483e04b9a4 Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Fri, 9 Oct 2026 11:55:52 +0200 Subject: [PATCH 03/15] Add comparison to varda single config --- config/varda-single-1.0.yaml | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/config/varda-single-1.0.yaml b/config/varda-single-1.0.yaml index 278026e5..aceeb2c8 100644 --- a/config/varda-single-1.0.yaml +++ b/config/varda-single-1.0.yaml @@ -120,7 +120,7 @@ showcase: - SP_10M - TOT_PREC1 meteograms: - enabled: false + enabled: true stations: [JUN] #, COV, GOR, WFJ, SAE, SAM, DAV, ZER, ANT, VSBAS, BRT, LTB, GOS, CEV, BIA] animations: enabled: true @@ -131,6 +131,14 @@ showcase: # - alps - icon-ch - switzerland + comparisons: + enabled: true + baselines: [ICON-CH1-ctrl] + domains: + - alps + params: + - SP_10M + lead_times: [24] locations: output_root: output/ From f0788cd53d4417a7a2ccbdb42b98c306ddd105dc Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Fri, 16 Oct 2026 08:28:07 +0200 Subject: [PATCH 04/15] Some cleanup and fixes following reviews --- config/varda-single-1.0.yaml | 5 +- src/evalml/config.py | 27 ++++--- src/plotting/__init__.py | 49 ++++++++++--- src/plotting/styles.py | 23 ++++-- workflow/Snakefile | 22 +++--- workflow/rules/common.smk | 3 +- workflow/rules/plot.smk | 63 ++++------------ ...are_maps.py => compare_forecast_frames.py} | 72 ++++++++++--------- workflow/scripts/plot_forecast_frame.py | 20 +----- workflow/tools/config.schema.json | 35 +++------ 10 files changed, 147 insertions(+), 172 deletions(-) rename workflow/scripts/{plot_compare_maps.py => compare_forecast_frames.py} (62%) diff --git a/config/varda-single-1.0.yaml b/config/varda-single-1.0.yaml index aceeb2c8..51e55378 100644 --- a/config/varda-single-1.0.yaml +++ b/config/varda-single-1.0.yaml @@ -133,11 +133,10 @@ showcase: - switzerland comparisons: enabled: true - baselines: [ICON-CH1-ctrl] + baselines: + - ICON-CH1-CTRL domains: - alps - params: - - SP_10M lead_times: [24] locations: diff --git a/src/evalml/config.py b/src/evalml/config.py index 55e5b8cc..f4de05a2 100644 --- a/src/evalml/config.py +++ b/src/evalml/config.py @@ -342,13 +342,13 @@ class ComparisonsConfig(BaseModel): enabled: bool = Field( default=False, - description="Whether to generate side-by-side maps of each run and the selected baselines.", + description="Whether to generate side-by-side maps of each run and each selected baseline.", ) - baselines: List[str] | None = Field( - default=None, + baselines: List[str] = Field( + default_factory=list, description=( - "Labels of the baselines to show next to each run. " - "None (default) means all baselines in the config." + "Labels of the baselines to compare each run with, one plot per baseline. " + "Required (non-empty) when comparisons are enabled." ), ) domains: List[str | DomainConfig] = Field( @@ -358,20 +358,25 @@ class ComparisonsConfig(BaseModel): "defined in plotting.DOMAINS or a custom domain dict." ), ) - params: List[str] = Field( - default=["SP_10M"], - description="Parameters to plot.", - ) lead_times: List[int] = Field( default=[24], description=( - "Lead times (hours) to plot. Lead times not produced by the run or by " - "a selected baseline are skipped." + "Lead times (hours) to plot. A run or baseline without data at a lead " + "time is shown as an empty panel." ), ) model_config = {"extra": "forbid"} + @model_validator(mode="after") + def _enabled_requires_baselines(self): + if self.enabled and not self.baselines: + raise ValueError( + "comparisons.baselines must list the baseline labels to compare " + "when comparisons are enabled." + ) + return self + class ScorecardConfig(BaseModel): """Configuration for a single named scorecard.""" diff --git a/src/plotting/__init__.py b/src/plotting/__init__.py index 63b6be12..38b61d9d 100644 --- a/src/plotting/__init__.py +++ b/src/plotting/__init__.py @@ -60,6 +60,32 @@ def get_projection(name: str) -> "ccrs.Projection": } +def resolve_region( + region_name: str, region_cfg: dict, lead_time: int = 0 +) -> tuple["ccrs.Projection", list[float] | None]: + """Return (projection, extent) for a showcase region. + + ``region_cfg`` is one entry of the dict built by ``parse_showcase_regions``. + A custom ``extent`` takes precedence; otherwise the named domain in + ``DOMAINS`` is used, rotated by ``lead_time`` when ``rotate`` is set. + """ + if region_cfg.get("extent") is not None: + projection = get_projection(region_cfg.get("projection") or "orthographic") + return projection, region_cfg["extent"] + domain = DOMAINS[region_name] + if region_cfg.get("rotate"): + base = domain["projection"].proj4_params + central_longitude = ( + base["lon_0"] + 360.0 * lead_time / region_cfg["hours_per_revolution"] + ) % 360 + # central_longitude=central_longitude, central_latitude=0.0 for zero-centered + projection = ccrs.Orthographic( + central_longitude=central_longitude, central_latitude=base["lat_0"] + ) + return projection, domain["extent"] + return domain["projection"], domain["extent"] + + class StatePlotter: """A class to plot state fields on various DOMAINS.""" @@ -208,16 +234,19 @@ def plot_field( # Normalize style and color-related kwargs style_to_use, plot_kwargs = self._prepare_plot_kwargs(style, kwargs) - # Temporarily suppress earthkit-plots internal source-based kwargs - with self._temporary_plot_kwargs_override(subplot): - subplot.tricontourf( - x=x[finite], - y=y[finite], - z=field[finite], - style=style_to_use, - transform=proj, - **plot_kwargs, - ) # for earthkit.plots to work properly cmap and norm are needed here + # An all-NaN field (e.g. accumulations at lead 0) leaves nothing to + # contour: only the map layers below are drawn. + if finite.any(): + # Temporarily suppress earthkit-plots internal source-based kwargs + with self._temporary_plot_kwargs_override(subplot): + subplot.tricontourf( + x=x[finite], + y=y[finite], + z=field[finite], + style=style_to_use, + transform=proj, + **plot_kwargs, + ) # for earthkit.plots to work properly cmap and norm are needed here # TODO: gridlines etc would be nicer to have in the init, but I didn't get # them to overlay the plot layer diff --git a/src/plotting/styles.py b/src/plotting/styles.py index d1bbb80f..db06b895 100644 --- a/src/plotting/styles.py +++ b/src/plotting/styles.py @@ -66,21 +66,32 @@ def preprocess_field(param: str, state: dict): Returns: (field_array, units_override or None) """ fields = state["fields"] - if param in ("T_2M", "TD_2M", "T", "TD"): - return kelvin_to_celsius(fields[param]), "°C" if param == "SP_10M": return ekm_wind.speed(fields["U_10M"], fields["V_10M"]), "m/s" if param == "SP": return ekm_wind.speed(fields["U"], fields["V"]), "m/s" + return prepare_values(param, fields[param]) + + +def prepare_values(param: str, values): + """Unit conversion and clipping of an already derived ``param`` field. + + Same as :func:`preprocess_field`, for callers that load derived params + (e.g. ``SP_10M``) directly. Returns: (field_array, units_override or None) + """ + if param in ("T_2M", "TD_2M", "T", "TD"): + return kelvin_to_celsius(values), "°C" + if param in ("SP_10M", "SP"): + return values, "m/s" if param == "TOT_PREC": - return np.maximum(fields[param], 0), "mm" + return np.maximum(values, 0), "mm" if param in ("CLCT", "CLCL"): # Avoid exact 0/1 plateaus breaking tricontourf on orthographic # projections (tmp/reproduce_clct_bug.py). Pair with extend="neither". # Any new bounded field with silent-blank or GeometryCollection-crash # globe frames likely needs the same clip-away-from-boundary fix. - return np.clip(fields[param], 1e-6, 1 - 1e-6), None + return np.clip(values, 1e-6, 1 - 1e-6), None if param == "SSRD": # Same issue, bottom boundary only (night-side plateau). - return np.maximum(fields[param], 1e-6), None - return fields[param], None + return np.maximum(values, 1e-6), None + return values, None diff --git a/workflow/Snakefile b/workflow/Snakefile index 68cc8ea6..43e058d2 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -225,19 +225,15 @@ rule showcase_all: else [] ), ( - [ - path - for run_id in CANDIDATES - for param in COMPARISONS_CONFIG["params"] - for path in expand( - rules.plot_compare_maps.output, - init_time=[t.strftime("%Y%m%d%H%M") for t in REFTIMES], - run_id=run_id, - param=param, - leadtime=get_compare_leadtimes(run_id, param), - showcase=EXPERIMENT_NAME, - ) - ] + expand( + rules.compare_forecast_frames.output, + init_time=[t.strftime("%Y%m%d%H%M") for t in REFTIMES], + run_id=CANDIDATES, + baseline=COMPARISONS_CONFIG["baselines"], + param=SHOWCASE_CONFIG["params"], + leadtime=COMPARISONS_CONFIG["lead_times"], + showcase=EXPERIMENT_NAME, + ) if COMPARISONS_CONFIG["enabled"] else [] ), diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 8b2be13c..d515135d 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -444,8 +444,7 @@ _comparisons = _showcase.get("comparisons", {}) COMPARISONS_CONFIG = { "enabled": _comparisons.get("enabled", False), "regions": parse_showcase_regions(_comparisons.get("domains", ["alps"])), - "params": _comparisons.get("params", ["SP_10M"]), - "baselines": _comparisons.get("baselines"), + "baselines": _comparisons.get("baselines", []), "lead_times": _comparisons.get("lead_times", [24]), } EXPERIMENT_PARAMS = config.get("experiment", {}).get( diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 70296532..5e5e49a3 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -216,41 +216,22 @@ use rule plot_scoremaps as plot_scoremaps_baseline with: / "logs/plot_scoremaps/{experiment}/{baseline_id}-{param}-{score}-{region}-{season}-{init_hour}-{leadtime}.log", -def _get_compare_baselines(wc) -> list[dict[str, str]]: - """Baselines selected for the side-by-side comparison maps.""" - selected = COMPARISONS_CONFIG["baselines"] - baselines = _get_available_baselines(wc) - if selected is None: - return baselines - return [b for b in baselines if b["label"] in selected] - - -def get_compare_leadtimes(run_id, param): - """Requested lead times produced by the run and by every selected baseline.""" - requested = COMPARISONS_CONFIG["lead_times"] - leadtimes = set(resolve_leadtimes(RUN_CONFIGS[run_id]["steps"], requested, param)) - selected = COMPARISONS_CONFIG["baselines"] - for cfg in BASELINE_CONFIGS.values(): - if selected is None or cfg.get("label") in selected: - leadtimes &= set(resolve_leadtimes(cfg["steps"], requested, param)) - return sorted(leadtimes) - - -rule plot_compare_maps: +rule compare_forecast_frames: input: - script="workflow/scripts/plot_compare_maps.py", + script="workflow/scripts/compare_forecast_frames.py", inference_okfile=rules.inference_execute.output.okfile, output: expand( OUT_ROOT - / "results/{{showcase}}/{{run_id}}/{{init_time}}/compare/{{init_time}}_{{leadtime}}_{{param}}_{region}.png", + / "results/{{showcase}}/{{run_id}}/{{init_time}}/compare/{{baseline}}/{{init_time}}_{{leadtime}}_{{param}}_{region}.png", region=list(COMPARISONS_CONFIG["regions"].keys()), ), log: OUT_ROOT - / "logs/{showcase}/{run_id}/{init_time}/plot_compare_maps_{leadtime}_{param}.log", + / "logs/{showcase}/{run_id}/{init_time}/compare_forecast_frames_{baseline}_{leadtime}_{param}.log", wildcard_constraints: leadtime=r"\d+", + baseline=r"[^/]+", resources: slurm_partition="postproc", cpus_per_task=4, @@ -260,37 +241,23 @@ rule plot_compare_maps: Path(OUT_ROOT) / f"data/runs/{wc.run_id}/{wc.init_time}/grib" ).resolve(), fcst_label=lambda wc: RUN_CONFIGS[wc.run_id]["label"], - baseline_roots=lambda wc: [x["root"] for x in _get_compare_baselines(wc)], - baseline_labels=lambda wc: [x["label"] for x in _get_compare_baselines(wc)], + baseline_root=lambda wc: BASELINE_CONFIGS[resolve_baseline_id(wc.baseline)][ + "root" + ], regions_json=json.dumps(COMPARISONS_CONFIG["regions"]), outdir=lambda wc: str( ( Path(OUT_ROOT) - / f"results/{wc.showcase}/{wc.run_id}/{wc.init_time}/compare" + / f"results/{wc.showcase}/{wc.run_id}/{wc.init_time}/compare/{wc.baseline}" ).resolve() ), shell: """ - set -euo pipefail export ECCODES_DEFINITION_PATH=$(realpath .venv/share/eccodes-cosmo-resources/definitions) - - BASELINE_ROOTS=({params.baseline_roots:q}) - BASELINE_LABELS=({params.baseline_labels:q}) - - CMD_ARGS=( - --forecast {params.fcst_grib:q} - --forecast_label {params.fcst_label:q} - --date {wildcards.init_time:q} - --leadtime {wildcards.leadtime:q} - --param {wildcards.param:q} - --regions_json {params.regions_json:q} - --outdir {params.outdir:q} - ) - - for i in "${{!BASELINE_ROOTS[@]}}"; do - CMD_ARGS+=(--baseline "${{BASELINE_ROOTS[$i]}}") - CMD_ARGS+=(--baseline_label "${{BASELINE_LABELS[$i]}}") - done - - uv run python {input.script} "${{CMD_ARGS[@]}}" >{log} 2>&1 + uv run python {input.script} \ + --forecast {params.fcst_grib:q} --forecast_label {params.fcst_label:q} \ + --baseline {params.baseline_root:q} --baseline_label {wildcards.baseline:q} \ + --date {wildcards.init_time:q} --leadtime {wildcards.leadtime:q} \ + --param {wildcards.param:q} --regions_json {params.regions_json:q} \ + --outdir {params.outdir:q} >{log} 2>&1 """ diff --git a/workflow/scripts/plot_compare_maps.py b/workflow/scripts/compare_forecast_frames.py similarity index 62% rename from workflow/scripts/plot_compare_maps.py rename to workflow/scripts/compare_forecast_frames.py index 72f30074..3aa88539 100644 --- a/workflow/scripts/plot_compare_maps.py +++ b/workflow/scripts/compare_forecast_frames.py @@ -1,4 +1,4 @@ -"""Plot side-by-side maps of an ML run and baselines at a common valid time.""" +"""Plot side-by-side maps of an ML run and a baseline at a common valid time.""" import json import logging @@ -8,15 +8,14 @@ from pathlib import Path import earthkit.plots as ekp -from earthkit.meteo.utils.convert import kelvin_to_celsius import numpy as np from data_input import load_forecast_data from data_input import parse_aggregated_param -from plotting import DOMAINS -from plotting import get_projection +from plotting import resolve_region from plotting import StatePlotter from plotting.styles import get_style +from plotting.styles import prepare_values LOG = logging.getLogger(__name__) logging.basicConfig( @@ -28,32 +27,19 @@ def load_source(root, init_time, lead_time, param): """Return (lon, lat, values) as flat arrays.""" ds = load_forecast_data(root, init_time, [lead_time], [param]).squeeze() - values = np.asarray(ds[param].values).ravel() - if param in ("T_2M", "TD_2M"): - values = kelvin_to_celsius(values) + base_param, _ = parse_aggregated_param(param) + values, _ = prepare_values(base_param, np.asarray(ds[param].values).ravel()) lon = np.asarray(ds["longitude"].values).ravel() lat = np.asarray(ds["latitude"].values).ravel() return lon, lat, values -def resolve_region(region_name, region_cfg): - """Return (projection, extent) for a showcase region.""" - if region_cfg.get("extent") is not None: - projection = get_projection(region_cfg.get("projection") or "orthographic") - return projection, region_cfg["extent"] - return DOMAINS[region_name]["projection"], DOMAINS[region_name]["extent"] - - def main(): parser = ArgumentParser() parser.add_argument("--forecast", type=str, help="Directory to ML run grib data") parser.add_argument("--forecast_label", type=str, default="Forecast") - parser.add_argument( - "--baseline", type=str, action="append", default=[], help="Baseline root" - ) - parser.add_argument( - "--baseline_label", type=str, action="append", default=[], help="Label" - ) + parser.add_argument("--baseline", type=str, help="Baseline root") + parser.add_argument("--baseline_label", type=str, default="Baseline") parser.add_argument("--date", type=str, help="reference datetime") parser.add_argument("--leadtime", type=int, help="leadtime") parser.add_argument("--param", type=str, help="parameter") @@ -73,32 +59,49 @@ def main(): outdir.mkdir(parents=True, exist_ok=True) valid_time = init_time + timedelta(hours=lead_time) - sources = [(args.forecast_label, args.forecast)] + list( - zip(args.baseline_label, args.baseline) - ) + sources = [ + (args.forecast_label, args.forecast), + (args.baseline_label, args.baseline), + ] data = {} for label, root in sources: LOG.info("Loading '%s' from %s", label, root) - data[label] = load_source(root, init_time, lead_time, param) + try: + lon, lat, values = load_source(root, init_time, lead_time, param) + except Exception: + # Loaders fail in different ways for lead times a source does not + # provide (outside its steps, no accumulation window, ...). + LOG.warning( + "Could not load '%s' at +%sh, plotting an empty panel", + label, + lead_time, + exc_info=True, + ) + continue + # Triangulate each grid once, shared by all regions. + data[label] = (StatePlotter(lon, lat, outdir), values) base_param, accu = parse_aggregated_param(param) - style = get_style(base_param, accu=accu or 1) - ncols = len(sources) + _, units = prepare_values(base_param, np.empty(0)) + style = get_style(base_param, units, accu=accu or 1) for region_name, region_cfg in regions.items(): LOG.info("Plotting region %s", region_name) - projection, extent = resolve_region(region_name, region_cfg) + projection, extent = resolve_region(region_name, region_cfg, lead_time) fig = ekp.Figure( crs=projection, domain=extent, rows=1, - columns=ncols, - size=(5 * ncols, 5), + columns=2, + size=(10, 5), ) for col, (label, _) in enumerate(sources): subplot = fig.add_map(row=0, column=col) - lon, lat, values = data[label] - plotter = StatePlotter(lon, lat, outdir) + if label not in data: + subplot.standard_layers() + subplot.title(f"{label} (not available)") + continue + plotter, values = data[label] plotter.plot_field( subplot, values, @@ -106,8 +109,9 @@ def main(): colorbar=False, **style, ) - # One colorbar shared by all panels, which use identical levels. - fig.legend(location="bottom") + if data: + # One colorbar shared by both panels, which use identical levels. + fig.legend(location="bottom") fig.title( f"{param} valid: {valid_time:%Y-%m-%d %H:%M} UTC " f"(init: {init_time:%Y-%m-%d %H:%M} UTC, +{lead_time}h)" diff --git a/workflow/scripts/plot_forecast_frame.py b/workflow/scripts/plot_forecast_frame.py index 880ec56d..754ca99e 100644 --- a/workflow/scripts/plot_forecast_frame.py +++ b/workflow/scripts/plot_forecast_frame.py @@ -5,8 +5,7 @@ import cartopy.crs as ccrs -from plotting import DOMAINS -from plotting import get_projection +from plotting import resolve_region from plotting import StatePlotter from plotting.compat import load_state_from_grib from plotting.styles import get_style @@ -108,22 +107,7 @@ def main(): LOG.info("Plotting region %s", region_name) outfn = outdir / f"frame_{lead_time}_{param}_{region_name}.png" plotter = StatePlotter(state["longitudes"], state["latitudes"], outdir) - if region_cfg.get("extent") is not None: - projection = get_projection(region_cfg.get("projection") or "orthographic") - extent = region_cfg["extent"] - elif region_cfg.get("rotate"): - base = DOMAINS[region_name]["projection"].proj4_params - central_longitude = ( - base["lon_0"] + 360.0 * lead_time / region_cfg["hours_per_revolution"] - ) % 360 - # central_longitude=central_longitude, central_latitude=0.0 for zero-centered - projection = ccrs.Orthographic( - central_longitude=central_longitude, central_latitude=base["lat_0"] - ) - extent = DOMAINS[region_name]["extent"] - else: - projection = DOMAINS[region_name]["projection"] - extent = DOMAINS[region_name]["extent"] + projection, extent = resolve_region(region_name, region_cfg, lead_time) fig = plotter.init_geoaxes( nrows=1, ncols=1, diff --git a/workflow/tools/config.schema.json b/workflow/tools/config.schema.json index 899d2a57..2f97e3cf 100644 --- a/workflow/tools/config.schema.json +++ b/workflow/tools/config.schema.json @@ -94,25 +94,17 @@ "properties": { "enabled": { "default": false, - "description": "Whether to generate side-by-side maps of each run and the selected baselines.", + "description": "Whether to generate side-by-side maps of each run and each selected baseline.", "title": "Enabled", "type": "boolean" }, "baselines": { - "anyOf": [ - { - "items": { - "type": "string" - }, - "type": "array" - }, - { - "type": "null" - } - ], - "default": null, - "description": "Labels of the baselines to show next to each run. None (default) means all baselines in the config.", - "title": "Baselines" + "description": "Labels of the baselines to compare each run with, one plot per baseline. Required (non-empty) when comparisons are enabled.", + "items": { + "type": "string" + }, + "title": "Baselines", + "type": "array" }, "domains": { "default": [ @@ -132,22 +124,11 @@ "title": "Domains", "type": "array" }, - "params": { - "default": [ - "SP_10M" - ], - "description": "Parameters to plot.", - "items": { - "type": "string" - }, - "title": "Params", - "type": "array" - }, "lead_times": { "default": [ 24 ], - "description": "Lead times (hours) to plot. Lead times not produced by the run or by a selected baseline are skipped.", + "description": "Lead times (hours) to plot. A run or baseline without data at a lead time is shown as an empty panel.", "items": { "type": "integer" }, From f159e187f5cafcaccb8dfbb889f7df655a791f08 Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Sat, 17 Oct 2026 00:15:08 +0200 Subject: [PATCH 05/15] Support aggregated params (e.g. TOT_PREC1) in showcase animations Frame plotting now decomposes aggregated params like the comparison maps do, and resolve_leadtimes only requests frames whose aggregation window is available. Importing data_input in common.smk loads earthkit.data, which exports ECCODES_ECKIT_GEO=1; unset it before running inference. --- workflow/rules/common.smk | 10 +++++++++- workflow/rules/inference.smk | 6 ++++++ workflow/scripts/plot_forecast_frame.py | 22 +++++++++++++++------- 3 files changed, 30 insertions(+), 8 deletions(-) diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index d515135d..ca89bfb4 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -6,6 +6,8 @@ import hashlib import json from urllib.parse import urlparse +from data_input import parse_aggregated_param + CONFIG_ROOT = Path("config").resolve() OUT_ROOT = Path(config["locations"]["output_root"]) @@ -521,6 +523,12 @@ def resolve_leadtimes(steps_spec, requested="all", param=None): ) valid = wanted & supported - if param in PARAMS_WITHOUT_STEP_ZERO_VALUE: + base_param, agg_hours = parse_aggregated_param(param) if param else (None, None) + if agg_hours: + # Aggregated params (e.g. TOT_PREC6) need the window start [lead - agg, lead] + valid = { + lt for lt in valid if lt - agg_hours == 0 or lt - agg_hours in supported + } + elif base_param in PARAMS_WITHOUT_STEP_ZERO_VALUE: valid = {lt for lt in valid if lt >= step} return sorted(valid) diff --git a/workflow/rules/inference.smk b/workflow/rules/inference.smk index e899ad4c..d675c74b 100644 --- a/workflow/rules/inference.smk +++ b/workflow/rules/inference.smk @@ -301,6 +301,12 @@ rule inference_execute: local VENV=/user-environment source "$VENV/bin/activate" + # Importing earthkit.data in the Snakemake process (e.g. via data_input) + # exports ECCODES_ECKIT_GEO=1, which srun would pass on. It makes ecCodes + # derive unstructured-grid coordinates that differ slightly from the + # checkpoint/dataset ones, and anemoi-inference then refuses to combine states. + unset ECCODES_ECKIT_GEO + if [ "{params.disable_local_definitions}" = "False" ]; then export ECCODES_DEFINITION_PATH="$VENV/share/eccodes-cosmo-resources/definitions" fi diff --git a/workflow/scripts/plot_forecast_frame.py b/workflow/scripts/plot_forecast_frame.py index 754ca99e..59b69167 100644 --- a/workflow/scripts/plot_forecast_frame.py +++ b/workflow/scripts/plot_forecast_frame.py @@ -5,6 +5,7 @@ import cartopy.crs as ccrs +from data_input import parse_aggregated_param from plotting import resolve_region from plotting import StatePlotter from plotting.compat import load_state_from_grib @@ -44,7 +45,9 @@ def main(): regions = json.loads(args.regions_json) outdir = Path(args.outdir) outdir.mkdir(parents=True, exist_ok=True) - accu = args.accu + # Aggregated params encode their period in the name (TOT_PREC1 → TOT_PREC, 1h) + base_param, agg_hours = parse_aggregated_param(param) + accu = agg_hours or args.accu LOG.info( "Plotting forecast frame: param=%s, init_time=%s, lead_time=%s, regions=%s", @@ -54,12 +57,12 @@ def main(): list(regions.keys()), ) - if param == "SP_10M": + if base_param == "SP_10M": paramlist = ["U_10M", "V_10M"] - elif param == "SP": + elif base_param == "SP": paramlist = ["U", "V"] else: - paramlist = [param] + paramlist = [base_param] # Load grib once — shared across all region plots # TODO: fix file pattern & globbing @@ -82,12 +85,17 @@ def main(): state = load_state_from_grib(grib_file, paramlist=paramlist) # tp is accumulated from start of forecast; de-accumulate to get period [lt-accu, lt] - if param == "TOT_PREC": + if base_param == "TOT_PREC": prev_lt = lead_time - accu if prev_lt > 0: prev_grib_files = list(grib_dir.glob(f"2*_{prev_lt}.grib")) if not prev_grib_files: prev_grib_files = list(grib_dir.glob(f"2*_{prev_lt:03d}.grib")) + if not prev_grib_files: + raise FileNotFoundError( + f"No GRIB file for lead_time={prev_lt} in {grib_dir}, needed to" + f" de-accumulate {param} at lead_time={lead_time}." + ) prev_grib_file = Path(prev_grib_files[0]) LOG.info( "De-accumulating TOT_PREC: loading previous grib file %s", @@ -100,7 +108,7 @@ def main(): ) # Preprocess field once — shared across all region plots - field, units_override = preprocess_field(param, state) + field, units_override = preprocess_field(base_param, state) validtime = state["valid_time"].strftime("%Y%m%d%H%M") for region_name, region_cfg in regions.items(): @@ -123,7 +131,7 @@ def main(): field, title=f"{param}, time: {validtime}", gridline_labels=not region_cfg.get("rotate", False), - **get_style(param, units_override, accu=accu), + **get_style(base_param, units_override, accu=accu), ) if len(state["lam_envelope"]) > 0: subplot.ax.add_geometries( From 84021869be549ec388aea1861c16c82f2eca21ad Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Wed, 21 Oct 2026 16:06:27 +0200 Subject: [PATCH 06/15] Compare truth, all runs and baselines in one figure One comparison figure per init time, lead time, parameter and domain, with a panel for the truth, every candidate run and each selected baseline (up to 4 per row, shared colorbar, max value in each panel title), instead of one figure per run and baseline pair. - comparisons.truth (default true) adds the truth panel; ignored for station truth, which cannot be drawn as a map. Baselines become optional; truth or baselines are required. - Outputs move to results///compare/. --- src/evalml/config.py | 28 ++++++--- tests/unit/test_compare_forecast_frames.py | 17 +++++ tests/unit/test_config.py | 16 +++++ workflow/Snakefile | 2 - workflow/rules/common.smk | 3 + workflow/rules/plot.smk | 58 ++++++++++++----- workflow/scripts/compare_forecast_frames.py | 69 ++++++++++++++------- workflow/tools/config.schema.json | 12 +++- 8 files changed, 154 insertions(+), 51 deletions(-) create mode 100644 tests/unit/test_compare_forecast_frames.py diff --git a/src/evalml/config.py b/src/evalml/config.py index 48ea833a..e9702a48 100644 --- a/src/evalml/config.py +++ b/src/evalml/config.py @@ -338,19 +338,27 @@ class AnimationsConfig(BaseModel): class ComparisonsConfig(BaseModel): - """Configuration for side-by-side map comparisons of runs and baselines.""" + """Configuration for side-by-side map comparisons of truth, runs and baselines.""" enabled: bool = Field( default=False, - description="Whether to generate side-by-side maps of each run and each selected baseline.", + description=( + "Whether to generate side-by-side maps: one figure per init time, lead " + "time, parameter and domain, with a panel for the truth (if gridded), " + "each run and each selected baseline." + ), ) - baselines: List[str] = Field( - default_factory=list, + truth: bool = Field( + default=True, description=( - "Labels of the baselines to compare each run with, one plot per baseline. " - "Required (non-empty) when comparisons are enabled." + "Whether to include a truth panel. Ignored for station truth " + "(jretrievedwh), which cannot be drawn as a map." ), ) + baselines: List[str] = Field( + default_factory=list, + description="Labels of the baselines to add as panels.", + ) domains: List[str | DomainConfig] = Field( default=["alps"], description=( @@ -369,11 +377,11 @@ class ComparisonsConfig(BaseModel): model_config = {"extra": "forbid"} @model_validator(mode="after") - def _enabled_requires_baselines(self): - if self.enabled and not self.baselines: + def _enabled_requires_a_reference(self): + if self.enabled and not (self.truth or self.baselines): raise ValueError( - "comparisons.baselines must list the baseline labels to compare " - "when comparisons are enabled." + "comparisons need something to compare the runs with: enable " + "`truth` or list `baselines`." ) return self diff --git a/tests/unit/test_compare_forecast_frames.py b/tests/unit/test_compare_forecast_frames.py new file mode 100644 index 00000000..cb3e9a79 --- /dev/null +++ b/tests/unit/test_compare_forecast_frames.py @@ -0,0 +1,17 @@ +import sys +from pathlib import Path + +import pytest + +sys.path.insert(0, str(Path(__file__).parents[2] / "workflow" / "scripts")) +from compare_forecast_frames import parse_source # noqa: E402 + + +def test_parse_source(): + assert parse_source("ICON-CH1=/store/x=y") == ("ICON-CH1", "/store/x=y") + + +@pytest.mark.parametrize("text", ["no-separator", "=root", "label="]) +def test_parse_source_rejects(text): + with pytest.raises(ValueError): + parse_source(text) diff --git a/tests/unit/test_config.py b/tests/unit/test_config.py index 530cfdf0..c64bd4f7 100644 --- a/tests/unit/test_config.py +++ b/tests/unit/test_config.py @@ -37,3 +37,19 @@ def test_experiment_optional(example_config): del example_config["experiment"] cfg = ConfigModel.model_validate(example_config) assert cfg.experiment is None + + +def test_comparisons_truth_only(example_config): + example_config["showcase"]["comparisons"] = {"enabled": True, "baselines": []} + cfg = ConfigModel.model_validate(example_config) + assert cfg.showcase.comparisons.truth is True + + +def test_comparisons_need_a_reference(example_config): + example_config["showcase"]["comparisons"] = { + "enabled": True, + "truth": False, + "baselines": [], + } + with pytest.raises(ValueError, match="enable `truth` or list `baselines`"): + ConfigModel.model_validate(example_config) diff --git a/workflow/Snakefile b/workflow/Snakefile index f61258cf..d70c61fd 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -226,8 +226,6 @@ rule showcase_all: expand( rules.compare_forecast_frames.output, init_time=[t.strftime("%Y%m%d%H%M") for t in REFTIMES], - run_id=CANDIDATES, - baseline=COMPARISONS_CONFIG["baselines"], param=SHOWCASE_CONFIG["params"], leadtime=COMPARISONS_CONFIG["lead_times"], showcase=EXPERIMENT_NAME, diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 7ab73be0..f28e4e34 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -474,6 +474,9 @@ COMPARISONS_CONFIG = { "regions": parse_showcase_regions(_comparisons.get("domains", ["alps"])), "baselines": _comparisons.get("baselines", []), "lead_times": _comparisons.get("lead_times", [24]), + # station truth (jretrievedwh) cannot be drawn as a map + "truth": _comparisons.get("truth", True) + and "jretrieve" not in str(config["truth"]["root"]), } EXPERIMENT_PARAMS = EXPERIMENT_CONFIG.get( "params", ["T_2M", "TD_2M", "SP_10M", "PS", "PMSL", "TOT_PREC6"] diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 19ce2656..182793d9 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -224,47 +224,75 @@ use rule plot_scoremaps as plot_scoremaps_baseline with: / "logs/plot_scoremaps/{experiment}/{baseline_id}-{param}-{score}-{region}-{season}-{init_hour}-{leadtime}.log", +def _comparison_sources(wc): + """'label=root' for each panel after the truth: all candidate runs, then the + selected baselines.""" + sources = [ + f"{RUN_CONFIGS[run_id]['label']}=" + f"{(Path(OUT_ROOT) / f'data/runs/{run_id}/{wc.init_time}/grib').resolve()}" + for run_id in CANDIDATES + ] + sources += [ + f"{label}={BASELINE_CONFIGS[resolve_baseline_id(label)]['root']}" + for label in COMPARISONS_CONFIG["baselines"] + ] + return sources + + rule compare_forecast_frames: input: script="workflow/scripts/compare_forecast_frames.py", - inference_okfile=rules.inference_execute.output.okfile, + inference_okfiles=lambda wc: expand( + rules.inference_execute.output.okfile, + run_id=CANDIDATES, + init_time=wc.init_time, + ), + truth=lambda wc: truth_file_dep(wc) if COMPARISONS_CONFIG["truth"] else [], output: expand( OUT_ROOT - / "results/{{showcase}}/{{run_id}}/{{init_time}}/compare/{{baseline}}/{{init_time}}_{{leadtime}}_{{param}}_{region}.png", + / "results/{{showcase}}/{{init_time}}/compare/{{init_time}}_{{leadtime}}_{{param}}_{region}.png", region=list(COMPARISONS_CONFIG["regions"].keys()), ), log: OUT_ROOT - / "logs/{showcase}/{run_id}/{init_time}/compare_forecast_frames_{baseline}_{leadtime}_{param}.log", + / "logs/{showcase}/{init_time}/compare_forecast_frames_{leadtime}_{param}.log", wildcard_constraints: + showcase=r"[^/]+", + init_time=r"\d{12}", leadtime=r"\d+", - baseline=r"[^/]+", resources: slurm_partition="postproc", cpus_per_task=4, + mem_mb=32_000, runtime="30m", params: - fcst_grib=lambda wc: ( - Path(OUT_ROOT) / f"data/runs/{wc.run_id}/{wc.init_time}/grib" - ).resolve(), - fcst_label=lambda wc: RUN_CONFIGS[wc.run_id]["label"], - baseline_root=lambda wc: BASELINE_CONFIGS[resolve_baseline_id(wc.baseline)][ - "root" - ], + sources=_comparison_sources, + truth=( + [ + "--truth", + config["truth"]["root"], + "--truth_label", + config["truth"]["label"], + ] + if COMPARISONS_CONFIG["truth"] + else [] + ), regions_json=json.dumps(COMPARISONS_CONFIG["regions"]), outdir=lambda wc: str( ( - Path(OUT_ROOT) - / f"results/{wc.showcase}/{wc.run_id}/{wc.init_time}/compare/{wc.baseline}" + Path(OUT_ROOT) / f"results/{wc.showcase}/{wc.init_time}/compare" ).resolve() ), shell: """ export ECCODES_DEFINITION_PATH=$(realpath .venv/share/eccodes-cosmo-resources/definitions) + SOURCE_ARGS=() + for s in {params.sources:q}; do + SOURCE_ARGS+=(--source "$s") + done uv run python {input.script} \ - --forecast {params.fcst_grib:q} --forecast_label {params.fcst_label:q} \ - --baseline {params.baseline_root:q} --baseline_label {wildcards.baseline:q} \ + "${{SOURCE_ARGS[@]}}" {params.truth:q} \ --date {wildcards.init_time:q} --leadtime {wildcards.leadtime:q} \ --param {wildcards.param:q} --regions_json {params.regions_json:q} \ --outdir {params.outdir:q} >{log} 2>&1 diff --git a/workflow/scripts/compare_forecast_frames.py b/workflow/scripts/compare_forecast_frames.py index 3aa88539..b65fd6b0 100644 --- a/workflow/scripts/compare_forecast_frames.py +++ b/workflow/scripts/compare_forecast_frames.py @@ -1,7 +1,8 @@ -"""Plot side-by-side maps of an ML run and a baseline at a common valid time.""" +"""Plot side-by-side maps of the truth, ML runs and baselines at a common valid time.""" import json import logging +import math from argparse import ArgumentParser from datetime import datetime from datetime import timedelta @@ -11,6 +12,7 @@ import numpy as np from data_input import load_forecast_data +from data_input import load_truth_data from data_input import parse_aggregated_param from plotting import resolve_region from plotting import StatePlotter @@ -24,9 +26,7 @@ ) -def load_source(root, init_time, lead_time, param): - """Return (lon, lat, values) as flat arrays.""" - ds = load_forecast_data(root, init_time, [lead_time], [param]).squeeze() +def _flat(ds, param): base_param, _ = parse_aggregated_param(param) values, _ = prepare_values(base_param, np.asarray(ds[param].values).ravel()) lon = np.asarray(ds["longitude"].values).ravel() @@ -34,12 +34,38 @@ def load_source(root, init_time, lead_time, param): return lon, lat, values +def load_source(root, init_time, lead_time, param): + """Return (lon, lat, values) as flat arrays.""" + ds = load_forecast_data(root, init_time, [lead_time], [param]).squeeze() + return _flat(ds, param) + + +def load_truth(root, init_time, lead_time, param): + """Return (lon, lat, values) of the truth at the valid time as flat arrays.""" + ds = load_truth_data(Path(root), init_time, [lead_time], [param]) + valid_time = np.datetime64(init_time + timedelta(hours=lead_time)) + return _flat(ds.sel(time=valid_time).squeeze(), param) + + +def parse_source(text): + """'label=root' -> (label, root); the label may not contain '='.""" + label, sep, root = text.partition("=") + if not sep or not label or not root: + raise ValueError(f"Expected 'label=root', got {text!r}") + return label, root + + def main(): parser = ArgumentParser() - parser.add_argument("--forecast", type=str, help="Directory to ML run grib data") - parser.add_argument("--forecast_label", type=str, default="Forecast") - parser.add_argument("--baseline", type=str, help="Baseline root") - parser.add_argument("--baseline_label", type=str, default="Baseline") + parser.add_argument( + "--source", + action="append", + default=[], + type=parse_source, + help="Panel as 'label=root' (ML run grib dir or baseline root), repeatable.", + ) + parser.add_argument("--truth", type=str, default=None, help="Truth root (zarr)") + parser.add_argument("--truth_label", type=str, default="Truth") parser.add_argument("--date", type=str, help="reference datetime") parser.add_argument("--leadtime", type=int, help="leadtime") parser.add_argument("--param", type=str, help="parameter") @@ -59,15 +85,14 @@ def main(): outdir.mkdir(parents=True, exist_ok=True) valid_time = init_time + timedelta(hours=lead_time) - sources = [ - (args.forecast_label, args.forecast), - (args.baseline_label, args.baseline), - ] + panels = [(label, root, load_source) for label, root in args.source] + if args.truth: + panels.insert(0, (args.truth_label, args.truth, load_truth)) data = {} - for label, root in sources: + for label, root, loader in panels: LOG.info("Loading '%s' from %s", label, root) try: - lon, lat, values = load_source(root, init_time, lead_time, param) + lon, lat, values = loader(root, init_time, lead_time, param) except Exception: # Loaders fail in different ways for lead times a source does not # provide (outside its steps, no accumulation window, ...). @@ -88,15 +113,17 @@ def main(): for region_name, region_cfg in regions.items(): LOG.info("Plotting region %s", region_name) projection, extent = resolve_region(region_name, region_cfg, lead_time) + columns = min(len(panels), 4) + rows = math.ceil(len(panels) / columns) fig = ekp.Figure( crs=projection, domain=extent, - rows=1, - columns=2, - size=(10, 5), + rows=rows, + columns=columns, + size=(5 * columns, 5 * rows), ) - for col, (label, _) in enumerate(sources): - subplot = fig.add_map(row=0, column=col) + for i, (label, _, _) in enumerate(panels): + subplot = fig.add_map(row=i // columns, column=i % columns) if label not in data: subplot.standard_layers() subplot.title(f"{label} (not available)") @@ -105,12 +132,12 @@ def main(): plotter.plot_field( subplot, values, - title=label, + title=f"{label} (max {np.nanmax(values):.1f} {units})", colorbar=False, **style, ) if data: - # One colorbar shared by both panels, which use identical levels. + # One colorbar shared by all panels, which use identical levels. fig.legend(location="bottom") fig.title( f"{param} valid: {valid_time:%Y-%m-%d %H:%M} UTC " diff --git a/workflow/tools/config.schema.json b/workflow/tools/config.schema.json index 850ec4b3..b464d6ad 100644 --- a/workflow/tools/config.schema.json +++ b/workflow/tools/config.schema.json @@ -90,16 +90,22 @@ }, "ComparisonsConfig": { "additionalProperties": false, - "description": "Configuration for side-by-side map comparisons of runs and baselines.", + "description": "Configuration for side-by-side map comparisons of truth, runs and baselines.", "properties": { "enabled": { "default": false, - "description": "Whether to generate side-by-side maps of each run and each selected baseline.", + "description": "Whether to generate side-by-side maps: one figure per init time, lead time, parameter and domain, with a panel for the truth (if gridded), each run and each selected baseline.", "title": "Enabled", "type": "boolean" }, + "truth": { + "default": true, + "description": "Whether to include a truth panel. Ignored for station truth (jretrievedwh), which cannot be drawn as a map.", + "title": "Truth", + "type": "boolean" + }, "baselines": { - "description": "Labels of the baselines to compare each run with, one plot per baseline. Required (non-empty) when comparisons are enabled.", + "description": "Labels of the baselines to add as panels.", "items": { "type": "string" }, From f6d8aa5633673cef440c8162ec5e610939b98097 Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Wed, 21 Oct 2026 16:42:13 +0200 Subject: [PATCH 07/15] Remove duplicated fix --- workflow/rules/inference.smk | 6 ------ 1 file changed, 6 deletions(-) diff --git a/workflow/rules/inference.smk b/workflow/rules/inference.smk index fcb88d96..da517e2c 100644 --- a/workflow/rules/inference.smk +++ b/workflow/rules/inference.smk @@ -305,12 +305,6 @@ rule inference_execute: local VENV=/user-environment source "$VENV/bin/activate" - # Importing earthkit.data in the Snakemake process (e.g. via data_input) - # exports ECCODES_ECKIT_GEO=1, which srun would pass on. It makes ecCodes - # derive unstructured-grid coordinates that differ slightly from the - # checkpoint/dataset ones, and anemoi-inference then refuses to combine states. - unset ECCODES_ECKIT_GEO - if [ "{params.disable_local_definitions}" = "False" ]; then export ECCODES_DEFINITION_PATH="$VENV/share/eccodes-cosmo-resources/definitions" fi From 20585f5727ed64d0241b5fa677db30d0d0133b2a Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Tue, 10 Nov 2026 13:55:31 +0100 Subject: [PATCH 08/15] fix: resolve comparison baselines by baseline label only --- src/evalml/config.py | 13 +++++++++++++ tests/unit/test_config.py | 24 ++++++++++++++++++++++++ workflow/rules/plot.smk | 4 +++- 3 files changed, 40 insertions(+), 1 deletion(-) diff --git a/src/evalml/config.py b/src/evalml/config.py index bdeb5831..919d62cd 100644 --- a/src/evalml/config.py +++ b/src/evalml/config.py @@ -1081,6 +1081,19 @@ def validate_scorecard_references(self) -> "ConfigModel": ) return self + @model_validator(mode="after") + def validate_comparison_baselines(self) -> "ConfigModel": + labels = [ + item.baseline.label for item in self.runs if isinstance(item, BaselineItem) + ] + unknown = [b for b in self.showcase.comparisons.baselines if b not in labels] + if unknown: + raise ValueError( + f"showcase.comparisons.baselines {unknown} do not match the label of " + f"any baseline in `runs`: {labels}" + ) + return self + model_config = { "extra": "forbid", # fail on misspelled keys "populate_by_name": True, diff --git a/tests/unit/test_config.py b/tests/unit/test_config.py index 0363d6de..11dd16f5 100644 --- a/tests/unit/test_config.py +++ b/tests/unit/test_config.py @@ -106,6 +106,30 @@ def test_comparisons_need_a_reference(example_config): ConfigModel.model_validate(example_config) +def test_comparison_baselines_must_be_baselines(example_config): + baseline_label = next( + item["baseline"]["label"] + for item in example_config["runs"] + if "baseline" in item + ) + forecaster_label = next( + item[k]["label"] + for item in example_config["runs"] + for k in item + if k in ("forecaster", "temporal_downscaler") + ) + example_config["showcase"]["comparisons"] = { + "enabled": True, + "baselines": [baseline_label], + } + ConfigModel.model_validate(example_config) + + for label in ("no-such-label", forecaster_label): + example_config["showcase"]["comparisons"]["baselines"] = [label] + with pytest.raises(ValueError, match="do not match the label of any baseline"): + ConfigModel.model_validate(example_config) + + def test_score_on_common_points(example_config): """score_on_common_points accepts a subset of experiment.params only.""" diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 81848a43..57c78aa1 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -232,8 +232,10 @@ def _comparison_sources(wc): f"{(Path(OUT_ROOT) / f'data/runs/{run_id}/{wc.init_time}/grib').resolve()}" for run_id in CANDIDATES ] + # baselines only: config validation checks the labels against them + baselines_by_label = {cfg["label"]: cfg for cfg in BASELINE_CONFIGS.values()} sources += [ - f"{label}={BASELINE_CONFIGS[resolve_baseline_id(label)]['root']}" + f"{label}={baselines_by_label[label]['root']}" for label in COMPARISONS_CONFIG["baselines"] ] return sources From 0b47406f65ef0e5c50a490dbe6c08188d0569df3 Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Tue, 10 Nov 2026 14:22:14 +0100 Subject: [PATCH 09/15] fix: clip cloud cover in percent in prepare_values --- src/plotting/styles.py | 45 +++++++++++++++++------------------------- 1 file changed, 18 insertions(+), 27 deletions(-) diff --git a/src/plotting/styles.py b/src/plotting/styles.py index c7bd5108..afa0c171 100644 --- a/src/plotting/styles.py +++ b/src/plotting/styles.py @@ -60,37 +60,27 @@ def get_style(param, units_override=None, accu=1): def preprocess_field(param: str, state: dict): """ - - Temperatures: K -> °C - Wind speed: sqrt(u^2 + v^2) - - Precipitation: m -> mm + - Everything else: see :func:`prepare_values` Returns: (field_array, units_override or None) """ fields = state["fields"] - if param in ("T_2M", "TD_2M", "TW_2M", "SKT", "T", "TD"): - return kelvin_to_celsius(fields[param]), "°C" if param == "SP_10M": - return ekm_wind.speed(fields["U_10M"], fields["V_10M"]), "m/s" - if param == "SP": - return ekm_wind.speed(fields["U"], fields["V"]), "m/s" - if param == "TOT_PREC": - return np.maximum(fields[param], 0), "mm" - if param in ("CLCT", "CLCL", "CLCM", "CLCH"): - # Avoid exact 0/100 plateaus breaking tricontourf on orthographic - # projections (tmp/reproduce_clct_bug.py). Pair with extend="neither". - # Any new bounded field with silent-blank or GeometryCollection-crash - # globe frames likely needs the same clip-away-from-boundary fix. - return np.clip(fields[param], 1e-4, 100 - 1e-4), None - if param == "SSRD": - # Same issue, bottom boundary only (night-side plateau). - return np.maximum(fields[param], 1e-6), None - return fields[param], None + values = ekm_wind.speed(fields["U_10M"], fields["V_10M"]) + elif param == "SP": + values = ekm_wind.speed(fields["U"], fields["V"]) + else: + values = fields[param] + return prepare_values(param, values) def prepare_values(param: str, values): """Unit conversion and clipping of an already derived ``param`` field. - Same as :func:`preprocess_field`, for callers that load derived params - (e.g. ``SP_10M``) directly. Returns: (field_array, units_override or None) + - Temperatures: K -> °C + - Precipitation: m -> mm + Used directly by callers that load derived params (e.g. ``SP_10M``). + Returns: (field_array, units_override or None) """ if param in ("T_2M", "TD_2M", "TW_2M", "SKT", "T", "TD"): return kelvin_to_celsius(values), "°C" @@ -98,12 +88,13 @@ def prepare_values(param: str, values): return values, "m/s" if param == "TOT_PREC": return np.maximum(values, 0), "mm" - if param in ("CLCT", "CLCL"): - # Avoid exact 0/1 plateaus breaking tricontourf on orthographic - # projections (tmp/reproduce_clct_bug.py). Pair with extend="neither". - # Any new bounded field with silent-blank or GeometryCollection-crash - # globe frames likely needs the same clip-away-from-boundary fix. - return np.clip(values, 1e-6, 1 - 1e-6), None + if param in ("CLCT", "CLCL", "CLCM", "CLCH"): + # Cloud cover is in percent. Avoid exact 0/100 plateaus breaking + # tricontourf on orthographic projections (tmp/reproduce_clct_bug.py). + # Pair with extend="neither". Any new bounded field with silent-blank + # or GeometryCollection-crash globe frames likely needs the same + # clip-away-from-boundary fix. + return np.clip(values, 1e-4, 100 - 1e-4), None if param == "SSRD": # Same issue, bottom boundary only (night-side plateau). return np.maximum(values, 1e-6), None From eaa416a7371231aab09546f1089a2202b36dfe10 Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Tue, 10 Nov 2026 14:46:12 +0100 Subject: [PATCH 10/15] fix: fail compare_forecast_frames on real loading errors --- src/data_input/__init__.py | 2 +- tests/unit/test_compare_forecast_frames.py | 51 +++++++++++++++++++++ workflow/rules/plot.smk | 16 +++++-- workflow/scripts/compare_forecast_frames.py | 37 +++++++++++---- 4 files changed, 91 insertions(+), 15 deletions(-) diff --git a/src/data_input/__init__.py b/src/data_input/__init__.py index 39df43bd..3c8c29ce 100644 --- a/src/data_input/__init__.py +++ b/src/data_input/__init__.py @@ -323,7 +323,7 @@ def _select_valid_times(ds, times: np.datetime64, strict: bool = False): ) return ds.sel(time=times_np[times_included]) else: - raise ValueError( + raise LookupError( "Valid times are not included in the dataset. " "Please check the valid times and the dataset." ) diff --git a/tests/unit/test_compare_forecast_frames.py b/tests/unit/test_compare_forecast_frames.py index cb3e9a79..ce836523 100644 --- a/tests/unit/test_compare_forecast_frames.py +++ b/tests/unit/test_compare_forecast_frames.py @@ -1,9 +1,11 @@ import sys from pathlib import Path +import numpy as np import pytest sys.path.insert(0, str(Path(__file__).parents[2] / "workflow" / "scripts")) +import compare_forecast_frames # noqa: E402 from compare_forecast_frames import parse_source # noqa: E402 @@ -15,3 +17,52 @@ def test_parse_source(): def test_parse_source_rejects(text): with pytest.raises(ValueError): parse_source(text) + + +@pytest.fixture +def run_main(monkeypatch, tmp_path): + """Run main() without plotting; loaders raise for roots listed in `errors`.""" + loaded = [] + + def run(args, errors): + def loader(root, *_): + if root in errors: + raise errors[root] + loaded.append(root) + return np.zeros(1), np.zeros(1), np.zeros(1) + + monkeypatch.setattr(compare_forecast_frames, "load_source", loader) + monkeypatch.setattr(compare_forecast_frames, "load_truth", loader) + monkeypatch.setattr(compare_forecast_frames, "StatePlotter", lambda *_: None) + monkeypatch.setattr( + sys, + "argv", + ["compare_forecast_frames.py", *args] + + ["--date", "202401010000", "--leadtime", "6", "--param", "T_2M"] + + ["--regions_json", "{}", "--outdir", str(tmp_path)], + ) + compare_forecast_frames.main() + return loaded + + return run + + +def test_missing_baseline_and_truth_are_skipped(run_main): + args = ["--source", "ML=run", "--baseline", "ICON=icon", "--truth", "truth"] + errors = {"icon": FileNotFoundError("icon"), "truth": LookupError("truth")} + assert run_main(args, errors) == ["run"] + + +def test_candidate_failure_raises(run_main): + with pytest.raises(FileNotFoundError): + run_main(["--source", "ML=run"], {"run": FileNotFoundError("run")}) + + +def test_unexpected_baseline_failure_raises(run_main): + with pytest.raises(ValueError): + run_main(["--baseline", "ICON=icon"], {"icon": ValueError("corrupt")}) + + +def test_nothing_loaded_raises(run_main): + with pytest.raises(RuntimeError): + run_main(["--baseline", "ICON=icon"], {"icon": FileNotFoundError("icon")}) diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 57c78aa1..7b41fbab 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -225,20 +225,22 @@ use rule plot_scoremaps as plot_scoremaps_baseline with: def _comparison_sources(wc): - """'label=root' for each panel after the truth: all candidate runs, then the - selected baselines.""" - sources = [ + """'label=root' for each candidate run panel, after the truth.""" + return [ f"{RUN_CONFIGS[run_id]['label']}=" f"{(Path(OUT_ROOT) / f'data/runs/{run_id}/{wc.init_time}/grib').resolve()}" for run_id in CANDIDATES ] + + +def _comparison_baselines(wc): + """'label=root' for each selected baseline panel, after the runs.""" # baselines only: config validation checks the labels against them baselines_by_label = {cfg["label"]: cfg for cfg in BASELINE_CONFIGS.values()} - sources += [ + return [ f"{label}={baselines_by_label[label]['root']}" for label in COMPARISONS_CONFIG["baselines"] ] - return sources rule compare_forecast_frames: @@ -270,6 +272,7 @@ rule compare_forecast_frames: runtime="30m", params: sources=_comparison_sources, + baselines=_comparison_baselines, truth=( [ "--truth", @@ -293,6 +296,9 @@ rule compare_forecast_frames: for s in {params.sources:q}; do SOURCE_ARGS+=(--source "$s") done + for s in {params.baselines:q}; do + SOURCE_ARGS+=(--baseline "$s") + done uv run python {input.script} \ "${{SOURCE_ARGS[@]}}" {params.truth:q} \ --date {wildcards.init_time:q} --leadtime {wildcards.leadtime:q} \ diff --git a/workflow/scripts/compare_forecast_frames.py b/workflow/scripts/compare_forecast_frames.py index b65fd6b0..fd383400 100644 --- a/workflow/scripts/compare_forecast_frames.py +++ b/workflow/scripts/compare_forecast_frames.py @@ -26,6 +26,12 @@ ) +# What the loaders raise when a baseline or the truth does not cover the lead +# time: missing ICON archive file, INCA step outside the file, valid time +# outside the truth zarr. Anything else is a real failure. +NOT_AVAILABLE = (FileNotFoundError, LookupError) + + def _flat(ds, param): base_param, _ = parse_aggregated_param(param) values, _ = prepare_values(base_param, np.asarray(ds[param].values).ravel()) @@ -62,7 +68,15 @@ def main(): action="append", default=[], type=parse_source, - help="Panel as 'label=root' (ML run grib dir or baseline root), repeatable.", + help="ML run panel as 'label=grib dir', repeatable. Must load.", + ) + parser.add_argument( + "--baseline", + action="append", + default=[], + type=parse_source, + help="Baseline panel as 'label=root', repeatable. " + "Left empty if it has no data at the lead time.", ) parser.add_argument("--truth", type=str, default=None, help="Truth root (zarr)") parser.add_argument("--truth_label", type=str, default="Truth") @@ -85,19 +99,22 @@ def main(): outdir.mkdir(parents=True, exist_ok=True) valid_time = init_time + timedelta(hours=lead_time) - panels = [(label, root, load_source) for label, root in args.source] + panels = [(label, root, load_source, False) for label, root in args.source] + panels += [(label, root, load_source, True) for label, root in args.baseline] if args.truth: - panels.insert(0, (args.truth_label, args.truth, load_truth)) + panels.insert(0, (args.truth_label, args.truth, load_truth, True)) + if not panels: + parser.error("no panels to plot") data = {} - for label, root, loader in panels: + for label, root, loader, optional in panels: LOG.info("Loading '%s' from %s", label, root) try: lon, lat, values = loader(root, init_time, lead_time, param) - except Exception: - # Loaders fail in different ways for lead times a source does not - # provide (outside its steps, no accumulation window, ...). + except NOT_AVAILABLE: + if not optional: + raise LOG.warning( - "Could not load '%s' at +%sh, plotting an empty panel", + "'%s' has no data at +%sh, plotting an empty panel", label, lead_time, exc_info=True, @@ -105,6 +122,8 @@ def main(): continue # Triangulate each grid once, shared by all regions. data[label] = (StatePlotter(lon, lat, outdir), values) + if not data: + raise RuntimeError(f"No source has data at +{lead_time}h") base_param, accu = parse_aggregated_param(param) _, units = prepare_values(base_param, np.empty(0)) @@ -122,7 +141,7 @@ def main(): columns=columns, size=(5 * columns, 5 * rows), ) - for i, (label, _, _) in enumerate(panels): + for i, (label, *_) in enumerate(panels): subplot = fig.add_map(row=i // columns, column=i % columns) if label not in data: subplot.standard_layers() From 2a93c7342d7d3b2de05afcf62792094be1a76908 Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Tue, 10 Nov 2026 14:50:04 +0100 Subject: [PATCH 11/15] fix: compute panel max over the plotted extent --- workflow/scripts/compare_forecast_frames.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/workflow/scripts/compare_forecast_frames.py b/workflow/scripts/compare_forecast_frames.py index fd383400..26c54dbc 100644 --- a/workflow/scripts/compare_forecast_frames.py +++ b/workflow/scripts/compare_forecast_frames.py @@ -8,6 +8,7 @@ from datetime import timedelta from pathlib import Path +import cartopy.crs as ccrs import earthkit.plots as ekp import numpy as np @@ -53,6 +54,16 @@ def load_truth(root, init_time, lead_time, param): return _flat(ds.sel(time=valid_time).squeeze(), param) +def max_in_view(subplot, lon, lat, values): + """Max of the values at points inside the subplot's visible map extent.""" + xyz = subplot._crs.transform_points(ccrs.PlateCarree(), lon, lat) + x, y = xyz[:, 0], xyz[:, 1] + x0, x1, y0, y1 = subplot.ax.get_extent() + # Points on the far side of an orthographic globe project to non-finite x/y. + inside = (x >= x0) & (x <= x1) & (y >= y0) & (y <= y1) & np.isfinite(values) + return np.max(values[inside]) if inside.any() else np.nan + + def parse_source(text): """'label=root' -> (label, root); the label may not contain '='.""" label, sep, root = text.partition("=") @@ -148,10 +159,11 @@ def main(): subplot.title(f"{label} (not available)") continue plotter, values = data[label] + vmax = max_in_view(subplot, plotter.lon, plotter.lat, values) plotter.plot_field( subplot, values, - title=f"{label} (max {np.nanmax(values):.1f} {units})", + title=f"{label} (max {vmax:.1f} {units})", colorbar=False, **style, ) From 2f09c0bc54f536cf996cf3acf5fd477dd158cbf9 Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Tue, 10 Nov 2026 15:01:32 +0100 Subject: [PATCH 12/15] fix: keep compare figure title clear of panel titles --- workflow/scripts/compare_forecast_frames.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/workflow/scripts/compare_forecast_frames.py b/workflow/scripts/compare_forecast_frames.py index 26c54dbc..a1cd6396 100644 --- a/workflow/scripts/compare_forecast_frames.py +++ b/workflow/scripts/compare_forecast_frames.py @@ -150,7 +150,7 @@ def main(): domain=extent, rows=rows, columns=columns, - size=(5 * columns, 5 * rows), + figsize=(5 * columns, 5 * rows), ) for i, (label, *_) in enumerate(panels): subplot = fig.add_map(row=i // columns, column=i % columns) @@ -170,7 +170,9 @@ def main(): if data: # One colorbar shared by all panels, which use identical levels. fig.legend(location="bottom") - fig.title( + # Figure.title() passes an explicit y, which constrained layout ignores + # when reserving space, so the suptitle can overlap the panel titles. + fig.fig.suptitle( f"{param} valid: {valid_time:%Y-%m-%d %H:%M} UTC " f"(init: {init_time:%Y-%m-%d %H:%M} UTC, +{lead_time}h)" ) From daee0007d0e7f29abe5c873bd785644fc141ceb3 Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Tue, 10 Nov 2026 15:20:45 +0100 Subject: [PATCH 13/15] Revert config changes --- config/varda-single-1.0.yaml | 7 ------- 1 file changed, 7 deletions(-) diff --git a/config/varda-single-1.0.yaml b/config/varda-single-1.0.yaml index 484e5112..d1736595 100644 --- a/config/varda-single-1.0.yaml +++ b/config/varda-single-1.0.yaml @@ -133,13 +133,6 @@ showcase: # - alps - icon-ch - switzerland - comparisons: - enabled: true - baselines: - - ICON-CH1-CTRL - domains: - - alps - lead_times: [24] locations: output_root: output/ From 6d3705c0f787b241b723fe97863c22d255cbd1a3 Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Tue, 10 Nov 2026 15:29:13 +0100 Subject: [PATCH 14/15] More fixes from review --- src/plotting/styles.py | 17 +++- ...es.py => test_plot_forecast_comparison.py} | 14 +-- workflow/Snakefile | 2 +- workflow/rules/common.smk | 96 +++---------------- workflow/rules/plot.smk | 10 +- ..._frames.py => plot_forecast_comparison.py} | 0 6 files changed, 45 insertions(+), 94 deletions(-) rename tests/unit/{test_compare_forecast_frames.py => test_plot_forecast_comparison.py} (80%) rename workflow/scripts/{compare_forecast_frames.py => plot_forecast_comparison.py} (100%) diff --git a/src/plotting/styles.py b/src/plotting/styles.py index afa0c171..68965bcb 100644 --- a/src/plotting/styles.py +++ b/src/plotting/styles.py @@ -9,9 +9,24 @@ from plotting.colormap_defaults import CMAP_DEFAULTS +def _precip_lookup(accu): + """Key of the TOT_PREC palette for an accumulation window of ``accu`` hours. + + Palettes exist only for some windows (1h, 6h); others use the longest + defined window not exceeding ``accu`` (e.g. 3h -> 1h, 24h -> 6h). + """ + windows = sorted( + int(key[len("TOT_PREC_") : -1]) + for key in CMAP_DEFAULTS + if key.startswith("TOT_PREC_") and key.endswith("H") + ) + window = max((w for w in windows if w <= accu), default=windows[0]) + return f"TOT_PREC_{window}H" + + def get_style(param, units_override=None, accu=1): """Get style and colormap settings for the plot.""" - lookup = f"{param}_{accu}H" if param == "TOT_PREC" else param + lookup = _precip_lookup(accu) if param == "TOT_PREC" else param cfg = CMAP_DEFAULTS[lookup] units = units_override if units_override is not None else cfg.get("units", "") diff --git a/tests/unit/test_compare_forecast_frames.py b/tests/unit/test_plot_forecast_comparison.py similarity index 80% rename from tests/unit/test_compare_forecast_frames.py rename to tests/unit/test_plot_forecast_comparison.py index ce836523..0ee06968 100644 --- a/tests/unit/test_compare_forecast_frames.py +++ b/tests/unit/test_plot_forecast_comparison.py @@ -5,8 +5,8 @@ import pytest sys.path.insert(0, str(Path(__file__).parents[2] / "workflow" / "scripts")) -import compare_forecast_frames # noqa: E402 -from compare_forecast_frames import parse_source # noqa: E402 +import plot_forecast_comparison # noqa: E402 +from plot_forecast_comparison import parse_source # noqa: E402 def test_parse_source(): @@ -31,17 +31,17 @@ def loader(root, *_): loaded.append(root) return np.zeros(1), np.zeros(1), np.zeros(1) - monkeypatch.setattr(compare_forecast_frames, "load_source", loader) - monkeypatch.setattr(compare_forecast_frames, "load_truth", loader) - monkeypatch.setattr(compare_forecast_frames, "StatePlotter", lambda *_: None) + monkeypatch.setattr(plot_forecast_comparison, "load_source", loader) + monkeypatch.setattr(plot_forecast_comparison, "load_truth", loader) + monkeypatch.setattr(plot_forecast_comparison, "StatePlotter", lambda *_: None) monkeypatch.setattr( sys, "argv", - ["compare_forecast_frames.py", *args] + ["plot_forecast_comparison.py", *args] + ["--date", "202401010000", "--leadtime", "6", "--param", "T_2M"] + ["--regions_json", "{}", "--outdir", str(tmp_path)], ) - compare_forecast_frames.main() + plot_forecast_comparison.main() return loaded return run diff --git a/workflow/Snakefile b/workflow/Snakefile index e3e9e4cf..7a3809bc 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -285,7 +285,7 @@ rule showcase_all: ), ( expand( - rules.compare_forecast_frames.output, + rules.plot_forecast_comparison.output, init_time=[t.strftime("%Y%m%d%H%M") for t in REFTIMES], param=SHOWCASE_CONFIG["params"], leadtime=COMPARISONS_CONFIG["lead_times"], diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index e1931a25..2ad240ec 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -8,7 +8,7 @@ from urllib.parse import urlparse from evalml.config import MODEL_REGISTRY_ROOT from data_input import inca, truth_type -from data_input import parse_aggregated_param +from evalml.leadtimes import resolve_leadtimes CONFIG_ROOT = Path("config").resolve() OUT_ROOT = Path(config["locations"]["output_root"]) @@ -137,11 +137,7 @@ def parse_showcase_regions(domains=None): (full-globe domains). """ if domains is None: - domains = ( - config.get("showcase", {}) - .get("animations", {}) - .get("domains", ["globe", "europe", "switzerland"]) - ) + domains = config["showcase"]["animations"]["domains"] result = {} for r in domains: if isinstance(r, str): @@ -153,10 +149,10 @@ def parse_showcase_regions(domains=None): } else: result[r["name"]] = { - "extent": r.get("extent"), - "projection": r.get("projection", "orthographic"), - "rotate": r.get("rotate", False), - "hours_per_revolution": r.get("hours_per_revolution", 96.0), + "extent": r["extent"], + "projection": r["projection"], + "rotate": r["rotate"], + "hours_per_revolution": r["hours_per_revolution"], } return result @@ -526,21 +522,20 @@ HOLDOUT_FILE = ( ) HOLDOUT_WORKDIR_NAME = "holdout_stations.csv" REGIONS = parse_regions() -_showcase = config.get("showcase", {}) +_showcase = config["showcase"] SHOWCASE_CONFIG = { "regions": parse_showcase_regions(), - "params": _showcase.get("params", ["T_2M", "SP_10M"]), - "fps": _showcase.get("animations", {}).get("frames_per_second", 2.0), + "params": _showcase["params"], + "fps": _showcase["animations"]["frames_per_second"], } -_comparisons = _showcase.get("comparisons", {}) +_comparisons = _showcase["comparisons"] COMPARISONS_CONFIG = { - "enabled": _comparisons.get("enabled", False), - "regions": parse_showcase_regions(_comparisons.get("domains", ["alps"])), - "baselines": _comparisons.get("baselines", []), - "lead_times": _comparisons.get("lead_times", [24]), + "enabled": _comparisons["enabled"], + "regions": parse_showcase_regions(_comparisons["domains"]), + "baselines": _comparisons["baselines"], + "lead_times": _comparisons["lead_times"], # station truth (jretrievedwh) cannot be drawn as a map - "truth": _comparisons.get("truth", True) - and "jretrieve" not in str(config["truth"]["root"]), + "truth": _comparisons["truth"] and TRUTH_TYPE != "jretrieve", } EXPERIMENT_PARAMS = EXPERIMENT_CONFIG.get( "params", ["T_2M", "TD_2M", "SP_10M", "PS", "PMSL", "TOT_PREC6"] @@ -571,31 +566,6 @@ SCORECARD_CONFIGS = ( ) -# Params with no value at lead time 0. Two distinct reasons land here: -# - period-accumulated params (TOT_PREC/tp) verify a [lead - period, lead] -# window, so they have no value at lead times shorter than one step spacing -# (e.g. no 0h precip map). -# - diagnostic params (e.g. CLCT/tcc in AIFS) aren't part of the model's input state, -# so they're simply absent from the initial-state GRIB file written at step 0. -# Cloud cover is prognostic in some models (Varda-single RC); listing it here -# only drops their 0h animation frame. -# Short and canonical names both appear across the workflow (showcases vs maps). -PARAMS_WITHOUT_STEP_ZERO_VALUE = { - "TOT_PREC", - "tp", - "CLCT", - "tcc", - "CLCL", - "lcc", - "CLCM", - "mcc", - "CLCH", - "hcc", - "SSRD", - "ssrd", -} - - def balfrin_cpus(mem_mb, actual_cpus=1, node_mem_mb=512_000, node_cpus=256): """Minimum cores to request on Balfrin based on memory usage. @@ -612,42 +582,6 @@ def balfrin_cpus(mem_mb, actual_cpus=1, node_mem_mb=512_000, node_cpus=256): return max(actual_cpus, mem_cpus) -def resolve_leadtimes(steps_spec, requested="all", param=None): - """Lead times to compute for a single participant. - - A run or baseline produces only the lead times in its own ``steps`` spec - (``start/stop/step``, hours). This returns those of the ``requested`` - selection that the participant actually produces — the literal ``"all"`` - (every produced lead time) or an explicit list of ints — so a 36h lead is - never requested of an ICON-CH1 baseline (steps ``0/33/6``), nor a >120h - lead of ICON-CH2. Explicitly requested lead times the participant cannot - produce are skipped with a warning. For accumulated ``param``s, lead times - shorter than one step spacing are dropped (no accumulation window). - """ - start, end, step = map(int, steps_spec.split("/")) - supported = set(range(start, end + 1, step)) - wanted = supported if requested == "all" else set(requested) - - unsupported = sorted(wanted - supported) - if unsupported: - logging.getLogger("snakemake").warning( - "Skipping lead time(s) %sh: not produced by forecast steps '%s'.", - unsupported, - steps_spec, - ) - - valid = wanted & supported - base_param, agg_hours = parse_aggregated_param(param) if param else (None, None) - if agg_hours: - # Aggregated params (e.g. TOT_PREC6) need the window start [lead - agg, lead] - valid = { - lt for lt in valid if lt - agg_hours == 0 or lt - agg_hours in supported - } - elif base_param in PARAMS_WITHOUT_STEP_ZERO_VALUE: - valid = {lt for lt in valid if lt >= step} - return sorted(valid) - - def scoremap_runtime(wc): """SLURM runtime (minutes) for a scoremap job, which loops over all reference times. Accumulated parameters (TOT_PREC*) re-accumulate hourly analysis diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 7b41fbab..7b30124e 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -168,6 +168,8 @@ rule make_forecast_animation: for f in input if Path(f).stat().st_size > 0 ] + if not frames: + raise ValueError(f"No non-empty frames to animate for {output[0]}") frames[0].save( output[0], save_all=True, @@ -243,9 +245,9 @@ def _comparison_baselines(wc): ] -rule compare_forecast_frames: +rule plot_forecast_comparison: input: - script="workflow/scripts/compare_forecast_frames.py", + script="workflow/scripts/plot_forecast_comparison.py", inference_okfiles=lambda wc: expand( rules.inference_execute.output.okfile, run_id=CANDIDATES, @@ -260,14 +262,14 @@ rule compare_forecast_frames: ), log: OUT_ROOT - / "logs/{showcase}/{init_time}/compare_forecast_frames_{leadtime}_{param}.log", + / "logs/plot_forecast_comparison/{showcase}-{init_time}-{leadtime}-{param}.log", wildcard_constraints: showcase=r"[^/]+", init_time=r"\d{12}", leadtime=r"\d+", resources: slurm_partition="postproc", - cpus_per_task=4, + cpus_per_task=balfrin_cpus(32_000, actual_cpus=4), mem_mb=32_000, runtime="30m", params: diff --git a/workflow/scripts/compare_forecast_frames.py b/workflow/scripts/plot_forecast_comparison.py similarity index 100% rename from workflow/scripts/compare_forecast_frames.py rename to workflow/scripts/plot_forecast_comparison.py From 0c79b94e1f87076e2afbd0c1d8f7e1ee255e16df Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Tue, 10 Nov 2026 15:30:39 +0100 Subject: [PATCH 15/15] Revert config change --- config/varda-single-1.0.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config/varda-single-1.0.yaml b/config/varda-single-1.0.yaml index d1736595..78c8c15f 100644 --- a/config/varda-single-1.0.yaml +++ b/config/varda-single-1.0.yaml @@ -122,7 +122,7 @@ showcase: - SP_10M - TOT_PREC1 meteograms: - enabled: true + enabled: false stations: [JUN] #, COV, GOR, WFJ, SAE, SAM, DAV, ZER, ANT, VSBAS, BRT, LTB, GOS, CEV, BIA] animations: enabled: true