[dnerini]

Add side-by-side map comparisons of runs and baselines to the showcase

New showcase.comparisons config section and plot_compare_maps rule that
plot each ML run next to selected baselines at the same valid time, with a
shared colour scale. Defaults: disabled, alps domain, SP_10M, and the lead
times produced by both the run and every selected baseline.

Move get_style() and preprocess_field() from plot_forecast_frame.py to
src/plotting/styles.py so both scripts share them.

Example output:

[Louis-Frey]

Hi Daniele,
Nice functionality, thanks! I'm testing the PR currently (feedback maybe next week).
At some point, Carlos was working on something similar, have you picked up from him?
Have a nice weekend!

[lclanzi]

Thanks @dnerini, this is indeed a very cool feature.

I have tested the PR and here are some initial comments I have

1. If we use two forecasters (say A and B) and an ICON baseline, the rule makes two figures, [A | ICON] and [B | ICON]. It would also be nice to have the option for a single [A | B | ICON] figure. This additional flexibility would require introducing a configurable option to control how the figures are generated.

2. The current output is a PNG at fixed lead times. An option to make a GIF over several lead times would be useful.

[Louis-Frey]

Hi Daniele,
I tested this and gave it to Claude for a review, this is what came out. Feel free to take anything from it that seems useful.

Test setup

• Config: copy of config/varda-single-1.0.yaml with comparisons of Varda-single-1.0 vs [ICON-CH1-CTRL, ICON-CH2-CTRL], domain icon-ch, param TOT_PREC1, lead_times: [0, 1, ..., 120], two inits (2025-07-02 06 UTC and 2025-12-02 06 UTC). Meteograms and animations disabled.
• Requesting 0-120h with ICON-CH1 (0-33h) in the comparison was on purpose, to test the lead-time intersection.

What works

• Lead-time intersection: only 0-33h are scheduled, and 34-120h are skipped with a clear warning.
• 66 of 68 plots were produced (+1h to +33h for both inits), about 2 min per job on a normal postproc node.

Bugs

1. TOT_PREC1 at +0h crashes the job (confirmed in the run).
Both +0h jobs fail with KeyError: "No variable named 'TOT_PREC1'. Variables on the dataset include []" in load_source (plot_compare_maps.py:31), which fails the whole showcase_all target. Cause: get_compare_leadtimes (plot.smk:231 and :235) passes the aggregated name TOT_PREC1 to resolve_leadtimes, which only drops step 0 if param in PARAMS_WITHOUT_STEP_ZERO_VALUE. That set holds base names (TOT_PREC), so lead 0 is kept. Passing the base name (via parse_aggregated_param) should fix it. The existing animation helper (plot.smk:141) has the same pattern, so this is not new in this PR, but the PR inherits it.

2. Case mismatch in the committed config (config/varda-single-1.0.yaml:136).
baselines: [ICON-CH1-ctrl], but the baseline label is ICON-CH1-CTRL. The filter is an exact match, so no baseline is selected and the "comparison" shows Varda alone, without any error or warning. Suggest fixing the label, and maybe raising an error (or at least a warning) for names in comparisons.baselines that match no baseline.

Probable issues (from reading the code, not tested)

3. Accumulations shorter than the window. Even with the base name fixed, resolve_leadtimes (common.smk:526) only drops lead times below the output step, not below the accumulation window. So for example TOT_PREC6 at +1h to +5h on hourly output, or TOT_PREC1 against a 6-hourly baseline, would likely hit the same empty-data crash as in 1.

4. Default baselines: None can silently remove all lead times. With baselines unset, get_compare_leadtimes (plot.smk:233-235) intersects over all baselines, including INCA (0-6h). With the default lead_times: [24] the result is empty, so nothing is plotted and nothing is reported. Related: the loop runs over BASELINE_CONFIGS, while _get_compare_baselines only keeps baselines with available data, so a baseline that is never plotted can still restrict the lead times.

Smaller points and suggestions

• load_source has its own preprocessing instead of reusing plotting.styles.preprocess_field, which this PR moves into styles.py. Today the two differ, for example in clipping negative precipitation and in the CLCT/CLCL/SSRD clipping that avoids blank or crashing globe maps. Reusing it would keep animations and comparisons consistent.
• resolve_region duplicates the region logic of plot_forecast_frame.py (without rotate). The schema says comparison domains use the same format as animations.domains, so a shared helper would avoid the two drifting apart.
• StatePlotter (which triangulates the full grid) is created inside the region x source loop (plot_compare_maps.py:101). Building it once per source, before the region loop, avoids repeating the triangulation for every region.
• One job per lead time means each job reopens all GRIB archives to draw one map. On a busy postproc node in my first run, jobs took about 11 min instead of 2 and two hit the 30 min limit. One job per (run, init, param) that loads all lead times once would be cheaper and more robust.
• Feature request: comparison animations. There are no GIFs for the comparison PNGs yet. Reusing the existing GIF step on the compare/ frames would make this consistent with the single-panel animations. Fine as a follow-up PR.

[dnerini]

If we use two forecasters (say A and B) and an ICON baseline, the rule makes two figures, [A | ICON] and [B | ICON]. It would also be nice to have the option for a single [A | B | ICON] figure. This additional flexibility would require introducing a configurable option to control how the figures are generated.

The current output is a PNG at fixed lead times. An option to make a GIF over several lead times would be useful.

thanks @lclanzi for your inputs! While I agree that both ideas could be useful, I'd suggest to implement them separately once we encounter the concrete need. This PR implements the plot I actually needed for a presentation and nothing more and I believe the code could be easily extended when needed.

[dnerini]

thanks @Louis-Frey for running this PR through Claude! I did my best to include the feedback, although some inputs were out of scope I think. See f0788cd

[lclanzi]

Thanks @dnerini!

[Louis-Frey]

Hi @dnerini, thanks for implementing the changes. I discussed the PR in its current state with Claude and we still have the following feedback:

1. Silent empty panels. compare_forecast_frames.py catches every error while loading data and just draws an empty panel. That's fine for "lead time not produced", but it also hides real problems: a wrong baseline path, a corrupt GRIB, or the Varda output missing. The run then reports success with blank maps. Maybe it's worth to at least fail when the run itself (not a baseline) can't be loaded, or to catch only the expected error types. This is a judgment call on your part I guess.

2. A side effect in f159e18 (minor, more of a question). common.smk now imports data_input, which loads earthkit, and that sets ECCODES_ECKIT_GEO=1 for the whole Snakemake process. You unset it only for inference; other rules would still inherit it. I haven't checked whether that matters. Main already imports data_input when truth is jretrieve, so this isn't completely new.

[dnerini]

updated code to allow for arbitrary runs, baselines, and truth, eg

[radiradev]

Would be nice to have this! There is a small issue with the max value in the plots - it is computed from the entire field, so it could be outside of the the region plotted.

[frazane]

@copilot resolve the merge conflicts in this pull request

[Copilot]

@copilot resolve the merge conflicts in this pull request

Resolved the merge conflicts in b2a56c6. Python syntax checks passed; targeted tests could not run because pytest and project dependencies are not installed.

[dnerini]

I tested this minimal config

# yaml-language-server: $schema=../workflow/tools/config.schema.json
description: |
  Minimal showcase to test side-by-side map comparisons (truth, runs, baselines)
  and animations of aggregated params, based on varda-single-1.0.yaml.

config_label: minimal-showcase-compare

dates:
  - 2025-03-01T00:00

runs:
  - temporal_downscaler:
      checkpoint: https://service.meteoswiss.ch/mlstore#/models/sruc-m-2-interpolator/versions/3
      label: Varda-single-1.0
      steps: 0/24/1
      config: resources/inference/configs/sgm-temporal-downscaler-global_trimedge_multi.yaml
      extra_requirements:
        - anemoi-datasets==0.5.35
        - git+https://github.com/ecmwf/anemoi-inference.git@88c1ec13632c570f246cc2580ce324e45f36a03b
        - git+https://github.com/MeteoSwiss/anemoi-plugins-meteoswiss.git@e5f95f3eb7b8907407faa7ac2e2771c1d7fe5aa4
      forecaster:
        checkpoint: https://service.meteoswiss.ch/mlstore#/models/varda-forecaster-det-sgm/versions/1
        config: resources/inference/configs/sgm-multidataset-forecaster-global-ich1-oper.yaml
        steps: 0/24/6
        extra_requirements:
          - anemoi-datasets==0.5.35
          - git+https://github.com/ecmwf/anemoi-inference.git@88c1ec13632c570f246cc2580ce324e45f36a03b
          - git+https://github.com/MeteoSwiss/anemoi-plugins-meteoswiss.git@e5f95f3eb7b8907407faa7ac2e2771c1d7fe5aa4
  - baseline:
      label: ICON-CH1-CTRL
      root: /store_new/mch/msopr/osm/ICON-CH1-EPS
      steps: 0/24/1

# Gridded truth, so that the comparison figures get a truth panel
# (station truth such as SwissMetNet is skipped).
truth:
  label: KENDA-CH1
  root: /store_new/mch/msopr/ml/datasets/mch-ich1-1km-2024-2025-1h-pl13-v1.0.zarr

# Required by the schema; only the showcase workflow is used for this config.
experiment:
  params:
    - SP_10M
  stratification:
    regions:
      - icon
    root: /store_new/mch/msopr/ml/regions/Prognoseregionen_LV95_20220517
  dashboard:
    stratification:
      - season

showcase:
  params:
    - SP_10M
    - TOT_PREC1
  meteograms:
    enabled: false
  animations:
    enabled: true
    frames_per_second: 5
    domains:
      - switzerland
  comparisons:
    enabled: true
    baselines: [ICON-CH1-CTRL]
    domains:
      - alps
    lead_times: [6, 24]

locations:
  output_root: output/

profile:
  executor: slurm
  global_resources:
    gpus: 16
  default_resources:
    slurm_partition: "postproc"
    cpus_per_task: 1
    mem_mb_per_cpu: 1800
    runtime: "1h"
    gpus: 0
    slurm_account: s83
  jobs: 50

which produced this comparison plot

[frazane]

Good to go from my side. Thanks for addressing the comments!