From 2b49526d073faa7d70484d79206c2df85a66c93e Mon Sep 17 00:00:00 2001 From: Jonas Bhend Date: Fri, 19 Dec 2025 14:44:06 +0100 Subject: [PATCH 01/58] Simplified way of implementing fields --- src/verification/__init__.py | 6 +++++- workflow/scripts/verif_single_init.py | 1 + 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/verification/__init__.py b/src/verification/__init__.py index 5f7c7a00..f0df171f 100644 --- a/src/verification/__init__.py +++ b/src/verification/__init__.py @@ -203,7 +203,11 @@ def verify( scores = _merge_metrics(scores) statistics = _merge_metrics(statistics) - out = xr.merge([scores, statistics], join="outer", compat="no_conflicts") + out = xr.merge( + [scores, statistics, fcst_aligned - obs_aligned], + join="outer", + compat="no_conflicts", + ) LOG.info("Computed metrics in %.2f seconds", time.time() - start) LOG.info("Metrics dataset: \n%s", out) return out diff --git a/workflow/scripts/verif_single_init.py b/workflow/scripts/verif_single_init.py index 57b1e2c4..52d9e08f 100644 --- a/workflow/scripts/verif_single_init.py +++ b/workflow/scripts/verif_single_init.py @@ -114,6 +114,7 @@ def main(args: ScriptConfig): # compute metrics and statistics results = verify(fcst, analysis, args.label, args.analysis_label, args.regions) + LOG.info("Verification results:\n%s", results) # save results to NetCDF args.output.parent.mkdir(parents=True, exist_ok=True) From a456b1d5e3806f58178f83c74f9263a05589ca9e Mon Sep 17 00:00:00 2001 From: Jonas Bhend Date: Fri, 19 Dec 2025 15:56:16 +0100 Subject: [PATCH 02/58] Exclude spatial data from being plotted and included in dashboard --- workflow/scripts/report_experiment_dashboard.py | 2 +- workflow/scripts/verif_plot_metrics.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/workflow/scripts/report_experiment_dashboard.py b/workflow/scripts/report_experiment_dashboard.py index 327547b0..a2abe766 100644 --- a/workflow/scripts/report_experiment_dashboard.py +++ b/workflow/scripts/report_experiment_dashboard.py @@ -39,7 +39,7 @@ def main(args): LOG.info("Loaded verification netcdf: \n%s", ds) # extract only non-spatial variables to pd.DataFrame - nonspatial_vars = [d for d in ds.data_vars if "spatial" not in d] + nonspatial_vars = [d for d in ds.data_vars if "x" not in ds[d].dims] df = ds[nonspatial_vars].to_array("stack").to_dataframe(name="value").reset_index() df[["param", "metric"]] = df["stack"].str.split(".", n=1, expand=True) df.drop(columns=["stack"], inplace=True) diff --git a/workflow/scripts/verif_plot_metrics.py b/workflow/scripts/verif_plot_metrics.py index a5a50bba..ae6d9d66 100644 --- a/workflow/scripts/verif_plot_metrics.py +++ b/workflow/scripts/verif_plot_metrics.py @@ -84,7 +84,7 @@ def main(args: Namespace) -> None: ds = xr.concat(dfs, dim="source", join="outer") # extract only non-spatial variables to pd.DataFrame - nonspatial_vars = [d for d in ds.data_vars if "spatial" not in d] + nonspatial_vars = [d for d in ds.data_vars if "x" not in ds[d].dims] all_df = ( ds[nonspatial_vars].to_array("stack").to_dataframe(name="value").reset_index() ) From 0ec286f4b24fe81f226bbd5ac9c56273ec6bceb8 Mon Sep 17 00:00:00 2001 From: Jonas Bhend Date: Wed, 7 Jan 2026 11:28:55 +0100 Subject: [PATCH 03/58] delete intermeidate verification files --- workflow/rules/verif.smk | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/workflow/rules/verif.smk b/workflow/rules/verif.smk index cba4a307..872f43eb 100644 --- a/workflow/rules/verif.smk +++ b/workflow/rules/verif.smk @@ -27,7 +27,7 @@ rule verif_metrics_baseline: analysis_label=config["analysis"].get("label"), regions=REGION_TXT, output: - OUT_ROOT / "data/baselines/{baseline_id}/{init_time}/verif.nc", + temp(OUT_ROOT / "data/baselines/{baseline_id}/{init_time}/verif.nc"), log: OUT_ROOT / "logs/verif_metrics_baseline/{baseline_id}-{init_time}.log", resources: @@ -76,7 +76,7 @@ rule verif_metrics: Path(OUT_ROOT) / f"data/runs/{wc.run_id}/{wc.init_time}/grib" ).resolve(), log: - OUT_ROOT / "logs/verif_metrics/{run_id}-{init_time}.log", + temp(OUT_ROOT / "logs/verif_metrics/{run_id}-{init_time}.log"), resources: cpus_per_task=24, mem_mb=50_000, From 28692d6c6adf890ab88554f6832d9ba9ded457c0 Mon Sep 17 00:00:00 2001 From: Jonas Bhend Date: Wed, 7 Jan 2026 12:02:58 +0100 Subject: [PATCH 04/58] Fix typo --- workflow/rules/verif.smk | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/workflow/rules/verif.smk b/workflow/rules/verif.smk index 872f43eb..0566ccf2 100644 --- a/workflow/rules/verif.smk +++ b/workflow/rules/verif.smk @@ -63,7 +63,7 @@ rule verif_metrics: inference_okfile=rules.execute_inference.output.okfile, analysis_zarr=config["analysis"].get("analysis_zarr"), output: - OUT_ROOT / "data/runs/{run_id}/{init_time}/verif.nc", + temp(OUT_ROOT / "data/runs/{run_id}/{init_time}/verif.nc"), # wildcard_constraints: # run_id="^" # to avoid ambiguitiy with run_baseline_verif # TODO: implement logic to use experiment name instead of run_id as wildcard @@ -76,7 +76,7 @@ rule verif_metrics: Path(OUT_ROOT) / f"data/runs/{wc.run_id}/{wc.init_time}/grib" ).resolve(), log: - temp(OUT_ROOT / "logs/verif_metrics/{run_id}-{init_time}.log"), + OUT_ROOT / "logs/verif_metrics/{run_id}-{init_time}.log", resources: cpus_per_task=24, mem_mb=50_000, From f3dcf0db6b5b53f6a2190590182a380deca0ab14 Mon Sep 17 00:00:00 2001 From: Jonas Bhend Date: Wed, 7 Jan 2026 12:04:04 +0100 Subject: [PATCH 05/58] include score components for maps --- src/verification/__init__.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/verification/__init__.py b/src/verification/__init__.py index f0df171f..51e483a4 100644 --- a/src/verification/__init__.py +++ b/src/verification/__init__.py @@ -198,13 +198,24 @@ def verify( score = xr.concat(score, dim="region") fcst_statistics = xr.concat(fcst_statistics, dim="region") obs_statistics = xr.concat(obs_statistics, dim="region") + score_spatial = _compute_scores( + fcst_aligned[param], + obs_aligned[param], + prefix=param + ".", + suffix=".spatial", + dim=[], + ) statistics.append(xr.concat([fcst_statistics, obs_statistics], dim="source")) - scores.append(score) + scores.append( + xr.merge([score, score_spatial], join="outer", compat="no_conflicts") + ) scores = _merge_metrics(scores) statistics = _merge_metrics(statistics) + LOG.info("Computed scores dataset: \n%s", scores) + LOG.info("Computed statistics dataset: \n%s", statistics) out = xr.merge( - [scores, statistics, fcst_aligned - obs_aligned], + [scores, statistics], join="outer", compat="no_conflicts", ) From 99dac523cdb71355c670088d4d15f72a1f8074bb Mon Sep 17 00:00:00 2001 From: Jonas Bhend Date: Wed, 7 Jan 2026 17:09:53 +0100 Subject: [PATCH 06/58] Revert "Exclude spatial data from being plotted and included in dashboard" This reverts commit cdefa16a7f0cc8d9721b7598e49a88800ab61d23. --- workflow/scripts/report_experiment_dashboard.py | 2 +- workflow/scripts/verif_plot_metrics.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/workflow/scripts/report_experiment_dashboard.py b/workflow/scripts/report_experiment_dashboard.py index a2abe766..327547b0 100644 --- a/workflow/scripts/report_experiment_dashboard.py +++ b/workflow/scripts/report_experiment_dashboard.py @@ -39,7 +39,7 @@ def main(args): LOG.info("Loaded verification netcdf: \n%s", ds) # extract only non-spatial variables to pd.DataFrame - nonspatial_vars = [d for d in ds.data_vars if "x" not in ds[d].dims] + nonspatial_vars = [d for d in ds.data_vars if "spatial" not in d] df = ds[nonspatial_vars].to_array("stack").to_dataframe(name="value").reset_index() df[["param", "metric"]] = df["stack"].str.split(".", n=1, expand=True) df.drop(columns=["stack"], inplace=True) diff --git a/workflow/scripts/verif_plot_metrics.py b/workflow/scripts/verif_plot_metrics.py index ae6d9d66..a5a50bba 100644 --- a/workflow/scripts/verif_plot_metrics.py +++ b/workflow/scripts/verif_plot_metrics.py @@ -84,7 +84,7 @@ def main(args: Namespace) -> None: ds = xr.concat(dfs, dim="source", join="outer") # extract only non-spatial variables to pd.DataFrame - nonspatial_vars = [d for d in ds.data_vars if "x" not in ds[d].dims] + nonspatial_vars = [d for d in ds.data_vars if "spatial" not in d] all_df = ( ds[nonspatial_vars].to_array("stack").to_dataframe(name="value").reset_index() ) From edcca5b6caff4cd470700c170ee3fe661e77fcaf Mon Sep 17 00:00:00 2001 From: Jonas Bhend Date: Wed, 7 Jan 2026 20:17:29 +0100 Subject: [PATCH 07/58] remove source dimension from scores --- src/verification/__init__.py | 33 ++++++++++++++++++++------------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/src/verification/__init__.py b/src/verification/__init__.py index 51e483a4..e9951b8a 100644 --- a/src/verification/__init__.py +++ b/src/verification/__init__.py @@ -84,17 +84,26 @@ def _compute_scores( """ dim = ["x", "y"] if "x" in fcst.dims and "y" in fcst.dims else ["values"] error = fcst - obs - scores = xr.Dataset( - { - f"{prefix}BIAS{suffix}": error.mean(dim=dim, skipna=True), - f"{prefix}MSE{suffix}": (error**2).mean(dim=dim, skipna=True), - f"{prefix}MAE{suffix}": abs(error).mean(dim=dim, skipna=True), - f"{prefix}VAR{suffix}": error.var(dim=dim, skipna=True), - f"{prefix}CORR{suffix}": xr.corr(fcst, obs, dim=dim), - f"{prefix}R2{suffix}": xr.corr(fcst, obs, dim=dim) ** 2, - } - ) - scores = scores.expand_dims({"source": [source]}) + if dim == []: + scores = xr.Dataset( + { + f"{prefix}BIAS{suffix}": error, + f"{prefix}MSE{suffix}": (error**2), + f"{prefix}MAE{suffix}": abs(error), + } + ) + else: + scores = xr.Dataset( + { + f"{prefix}BIAS{suffix}": error.mean(dim=dim, skipna=True), + f"{prefix}MSE{suffix}": (error**2).mean(dim=dim, skipna=True), + f"{prefix}MAE{suffix}": abs(error).mean(dim=dim, skipna=True), + f"{prefix}VAR{suffix}": error.var(dim=dim, skipna=True), + f"{prefix}CORR{suffix}": xr.corr(fcst, obs, dim=dim), + f"{prefix}R2{suffix}": xr.corr(fcst, obs, dim=dim) ** 2, + } + ) + # scores = scores.expand_dims({"source": [source]}) return scores @@ -212,8 +221,6 @@ def verify( scores = _merge_metrics(scores) statistics = _merge_metrics(statistics) - LOG.info("Computed scores dataset: \n%s", scores) - LOG.info("Computed statistics dataset: \n%s", statistics) out = xr.merge( [scores, statistics], join="outer", From 0ec8a0f926bd4b641c1436aa4793e1c8592da515 Mon Sep 17 00:00:00 2001 From: Jonas Bhend Date: Thu, 8 Jan 2026 08:40:39 +0100 Subject: [PATCH 08/58] clean up --- src/verification/__init__.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/verification/__init__.py b/src/verification/__init__.py index e9951b8a..8f05349c 100644 --- a/src/verification/__init__.py +++ b/src/verification/__init__.py @@ -221,11 +221,7 @@ def verify( scores = _merge_metrics(scores) statistics = _merge_metrics(statistics) - out = xr.merge( - [scores, statistics], - join="outer", - compat="no_conflicts", - ) + out = xr.merge([scores, statistics], join="outer", compat="no_conflicts") LOG.info("Computed metrics in %.2f seconds", time.time() - start) LOG.info("Metrics dataset: \n%s", out) return out From 51f22a6e44cf5f2660a7813cba1bbf73e7ab4614 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Mon, 12 Jan 2026 16:28:01 +0100 Subject: [PATCH 09/58] New rule and plotting file for plotting maps of summary statistics. (No changes to code yet.) --- workflow/rules/plot.smk | 32 +++ workflow/scripts/plot_summary_stat_maps.mo.py | 239 ++++++++++++++++++ 2 files changed, 271 insertions(+) create mode 100644 workflow/scripts/plot_summary_stat_maps.mo.py diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index d4def1ba..7fbd3320 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -67,3 +67,35 @@ rule make_forecast_animation: """ convert -delay {params.delay} -loop 0 {input} {output} """ + + +rule plot_summary_stat_maps: + input: + script="workflow/scripts/plot_forecast_frame.mo.py", + inference_okfile=rules.execute_inference.output.okfile, + output: + temp( + OUT_ROOT + / "showcases/{run_id}/{init_time}/frames/{init_time}_{leadtime}_{param}_{region}.png" + ), + wildcard_constraints: + leadtime=r"\d+", # only digits + resources: + slurm_partition="postproc", + cpus_per_task=1, + runtime="10m", + params: + grib_out_dir=lambda wc: ( + Path(OUT_ROOT) / f"data/runs/{wc.run_id}/{wc.init_time}/grib" + ).resolve(), + shell: + """ + export ECCODES_DEFINITION_PATH=$(realpath .venv/share/eccodes-cosmo-resources/definitions) + python {input.script} \ + --input {params.grib_out_dir} --date {wildcards.init_time} --outfn {output[0]} \ + --param {wildcards.param} --leadtime {wildcards.leadtime} --region {wildcards.region} \ + # interactive editing (needs to set localrule: True and use only one core) + # marimo edit {input.script} -- \ + # --input {params.grib_out_dir} --date {wildcards.init_time} --outfn {output[0]}\ + # --param {wildcards.param} --leadtime {wildcards.leadtime} --region {wildcards.region}\ + """ \ No newline at end of file diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py new file mode 100644 index 00000000..b8592141 --- /dev/null +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -0,0 +1,239 @@ +import marimo + +__generated_with = "0.16.5" +app = marimo.App(width="medium") + + +@app.cell +def _(): + import logging + from argparse import ArgumentParser + from pathlib import Path + + import cartopy.crs as ccrs + import earthkit.plots as ekp + import numpy as np + + from plotting import DOMAINS + from plotting import StatePlotter + from plotting.colormap_defaults import CMAP_DEFAULTS + from plotting.compat import load_state_from_grib + + return ( + ArgumentParser, + CMAP_DEFAULTS, + Path, + StatePlotter, + ekp, + load_state_from_grib, + logging, + np, + DOMAINS, + ccrs, + ) + + +@app.cell +def _(logging): + LOG = logging.getLogger(__name__) + LOG_FMT = "%(asctime)s - %(name)s - %(levelname)s - %(message)s" + logging.basicConfig(level=logging.INFO, format=LOG_FMT) + return (LOG,) + + +@app.cell +def _(ArgumentParser, Path): + parser = ArgumentParser() + + parser.add_argument( + "--input", type=str, default=None, help="Directory to grib data" + ) + parser.add_argument("--date", type=str, default=None, help="reference datetime") + parser.add_argument("--outfn", type=str, help="output filename") + parser.add_argument("--leadtime", type=str, help="leadtime") + parser.add_argument("--param", type=str, help="parameter") + parser.add_argument("--region", type=str, help="name of region") + + args = parser.parse_args() + grib_dir = Path(args.input) + init_time = args.date + outfn = Path(args.outfn) + lead_time = args.leadtime + param = args.param + region = args.region + return ( + args, + grib_dir, + init_time, + lead_time, + outfn, + param, + region, + ) + + +@app.cell +def _(grib_dir, init_time, lead_time, load_state_from_grib, param): + # load grib file + grib_file = grib_dir / f"{init_time}_{lead_time}.grib" + if param == "SP_10M": + paramlist = ["U_10M", "V_10M"] + elif param == "SP": + paramlist = ["U", "V"] + else: + paramlist = [param] + state = load_state_from_grib(grib_file, paramlist=paramlist) + return (state,) + + +@app.cell +def _(CMAP_DEFAULTS, ekp): + def get_style(param, units_override=None): + """Get style and colormap settings for the plot. + Needed because cmap/norm does not work in Style(colors=cmap), + still needs to be passed as arguments to tripcolor()/tricontourf(). + """ + cfg = CMAP_DEFAULTS[param] + units = units_override if units_override is not None else cfg.get("units", "") + return { + "style": ekp.styles.Style( + levels=cfg.get("bounds", cfg.get("levels", None)), + extend="both", + units=units, + colors=cfg.get("colors", None), + ), + "norm": cfg.get("norm", None), + "cmap": cfg.get("cmap", None), + "levels": cfg.get("levels", None), + "vmin": cfg.get("vmin", None), + "vmax": cfg.get("vmax", None), + "colors": cfg.get("colors", None), + } + + return (get_style,) + + +@app.cell +def _(LOG, np): + """Preprocess fields with pint-based unit conversion and derived quantities.""" + try: + import pint # type: ignore + + _ureg = pint.UnitRegistry() + + def _k_to_c(arr): + # robust conversion with pint, fallback if dtype unsupported + try: + return (_ureg.Quantity(arr, _ureg.kelvin).to(_ureg.degC)).magnitude + except Exception: + return arr - 273.15 + + def _ms_to_knots(arr): + # robust conversion with pint, fallback if dtype unsupported + try: + return ( + _ureg.Quantity(arr, _ureg.meter / _ureg.second).to(_ureg.knot) + ).magnitude + except Exception: + return arr * 1.943844 + + def _m_to_mm(arr): + # robust conversion with pint, fallback if dtype unsupported + try: + return (_ureg.Quantity(arr, _ureg.meter).to(_ureg.millimeter)).magnitude + except Exception: + return arr * 1000 + + except Exception: + LOG.warning("pint not available; falling back hardcoded conversions") + + def _k_to_c(arr): + return arr - 273.15 + + def _ms_to_knots(arr): + return arr * 1.943844 + + def _m_to_mm(arr): + return arr * 1000 + + 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"] + # temperature variables + if param in ("T_2M", "TD_2M", "T", "TD"): + return _k_to_c(fields[param]), "°C" + # 10m wind speed (allow legacy 'uv' alias) + if param == "SP_10M": + u = fields["U_10M"] + v = fields["V_10M"] + return np.sqrt(u**2 + v**2), "m/s" + # wind speed from standard-level components + if param == "SP": + u = fields["U"] + v = fields["V"] + return np.sqrt(u**2 + v**2), "m/s" + if param == "TOT_PREC": + return _m_to_mm(fields[param]), "mm" + # default: passthrough + return fields[param], None + + return (preprocess_field,) + + +@app.cell +def _( + LOG, + StatePlotter, + args, + get_style, + outfn, + param, + preprocess_field, + region, + state, + DOMAINS, + ccrs, +): + # plot individual fields + plotter = StatePlotter( + state["longitudes"], + state["latitudes"], + outfn.parent, + ) + fig = plotter.init_geoaxes( + nrows=1, + ncols=1, + projection=DOMAINS[region]["projection"], + bbox=DOMAINS[region]["extent"], + name=region, + size=(6, 6), + ) + subplot = fig.add_map(row=0, column=0) + + # preprocess field (unit conversion, derived quantities) + field, units_override = preprocess_field(param, state) + + plotter.plot_field(subplot, field, **get_style(args.param, units_override)) + subplot.ax.add_geometries( + state["lam_envelope"], + edgecolor="black", + facecolor="none", + crs=ccrs.PlateCarree(), + ) + validtime = state["valid_time"].strftime("%Y%m%d%H%M") + # leadtime = int(state["lead_time"].total_seconds() // 3600) + + fig.title(f"{param}, time: {validtime}") + + fig.save(outfn, bbox_inches="tight", dpi=200) + LOG.info(f"saved: {outfn}") + return + + +if __name__ == "__main__": + app.run() From 366249d5563caa9dc1d795d25ee84f65b35fe952 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Mon, 12 Jan 2026 17:01:42 +0100 Subject: [PATCH 10/58] Obvious changes to the new plotting rule for maps of summary statistics. --- workflow/rules/plot.smk | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 7fbd3320..fee1dfe3 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -71,13 +71,10 @@ rule make_forecast_animation: rule plot_summary_stat_maps: input: - script="workflow/scripts/plot_forecast_frame.mo.py", + script="workflow/scripts/plot_summary_stat_maps.mo.py", inference_okfile=rules.execute_inference.output.okfile, output: - temp( - OUT_ROOT - / "showcases/{run_id}/{init_time}/frames/{init_time}_{leadtime}_{param}_{region}.png" - ), + OUT_ROOT / "results/summary_stats/maps/{run_id}/{leadtime}/{metric}_{param}_{region}.png", wildcard_constraints: leadtime=r"\d+", # only digits resources: @@ -85,9 +82,7 @@ rule plot_summary_stat_maps: cpus_per_task=1, runtime="10m", params: - grib_out_dir=lambda wc: ( - Path(OUT_ROOT) / f"data/runs/{wc.run_id}/{wc.init_time}/grib" - ).resolve(), + # What do I need here? shell: """ export ECCODES_DEFINITION_PATH=$(realpath .venv/share/eccodes-cosmo-resources/definitions) From f53139519f83168321f6db96f19e228772958aec Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Mon, 12 Jan 2026 17:30:14 +0100 Subject: [PATCH 11/58] Some more changes (preliminary, to be continued). --- workflow/rules/plot.smk | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index fee1dfe3..8090a213 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -82,12 +82,14 @@ rule plot_summary_stat_maps: cpus_per_task=1, runtime="10m", params: - # What do I need here? + nc_out_dir=lambda wc: ( + Path(OUT_ROOT) / f"data/runs/{wc.run_id}/{wc.init_time}/grib" # to be adjusted. + ).resolve(), shell: """ export ECCODES_DEFINITION_PATH=$(realpath .venv/share/eccodes-cosmo-resources/definitions) python {input.script} \ - --input {params.grib_out_dir} --date {wildcards.init_time} --outfn {output[0]} \ + --input {params.nc_out_dir} --date {wildcards.init_time} --outfn {output[0]} \ --param {wildcards.param} --leadtime {wildcards.leadtime} --region {wildcards.region} \ # interactive editing (needs to set localrule: True and use only one core) # marimo edit {input.script} -- \ From 32c7ecf6ff676c9c978090f44208bfd1e53ab558 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Tue, 13 Jan 2026 18:02:48 +0100 Subject: [PATCH 12/58] Further changes to plotting scripts. --- workflow/rules/plot.smk | 4 +++- workflow/scripts/plot_summary_stat_maps.mo.py | 14 +++++++++++++- 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 8090a213..2014e000 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -83,7 +83,9 @@ rule plot_summary_stat_maps: runtime="10m", params: nc_out_dir=lambda wc: ( - Path(OUT_ROOT) / f"data/runs/{wc.run_id}/{wc.init_time}/grib" # to be adjusted. + Path(OUT_ROOT) / f"data/runs/{wc.run_id}/{wc.init_time}/grib" + # not sure how to do this, because the baselines are in, e.g., output/data/baselines/COSMO-E/verif_aggregated.nc + # and the runs are in output/data/runs/runID/verif_aggregated.nc ).resolve(), shell: """ diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index b8592141..4b8eb38f 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -6,18 +6,30 @@ @app.cell def _(): + + # this sure stays the same. import logging from argparse import ArgumentParser from pathlib import Path + # this sure stays the same. import cartopy.crs as ccrs import earthkit.plots as ekp import numpy as np + # this stays the same as well. from plotting import DOMAINS + + # no changes to StatePlotter required according to ChatGPT. from plotting import StatePlotter + + # Probably need some new colour maps. + # at least one for biases (diverging), maybe different diverging ones for the different variables. + # from plotting.colormap_defaults import CMAP_DEFAULTS - from plotting.compat import load_state_from_grib + + # need to load nc files... + from plotting.compat import load_state_from_grib # TODO: load state from nc? return ( ArgumentParser, From c2ab6451d822271cdd4a9598c3dec91ae118134c Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Wed, 14 Jan 2026 15:12:43 +0100 Subject: [PATCH 13/58] First version of colour maps finished. For Bias, RMSE and MAE map plots. --- src/plotting/colormap_defaults.py | 45 ++++++++++++++++++++++++++++++- 1 file changed, 44 insertions(+), 1 deletion(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index 52ff0502..4a85f67e 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -4,12 +4,23 @@ from matplotlib import pyplot as plt import warnings from .colormap_loader import load_ncl_colormap - +from matplotlib.colors import BoundaryNorm +import numpy as np def _fallback(): warnings.warn("No colormap found for this parameter, using fallback.", UserWarning) return {"cmap": plt.get_cmap("viridis"), "norm": None, "units": ""} +def symmetric_boundary_norm(nlevels): + """ + Returns a callable that creates a symmetric BoundaryNorm + around zero with `nlevels` discrete colors. Used for creating colormaps for bias. + """ + def _norm(data): + vmax = np.nanmax(np.abs(data)) + boundaries = np.linspace(-vmax, vmax, nlevels + 1) + return BoundaryNorm(boundaries=boundaries, ncolors=nlevels) + return _norm _CMAP_DEFAULTS = { "SP": {"cmap": plt.get_cmap("coolwarm", 11), "vmin": 800 * 100, "vmax": 1100 * 100}, @@ -44,6 +55,38 @@ def _fallback(): "units": "mm", "levels": [0, 0.05, 0.1, 0.25, 0.5, 1, 1.5, 2, 3, 4, 5, 6, 7, 100], }, + + # hard-code this for the moment, can still make smarter later on: + # RMSE and MAE first (is all the same). Use Reds colormap to indicate 'error'. + # Use a limited number of levels so that absolute values of error can be read from the map. + # always start at 0 so that the saturation of the colour corresponds to the error magnitude. + + # RMSE: + "U_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "m/s"}, + "V_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "m/s"}, + "TD_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "°C"}, + "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "°C"}, + "PMSL.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "Pa"}, + "PS.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "Pa"}, + "TOT_PREC.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "mm"}, + + # MAE: + "U_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "m/s"}, + "V_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "m/s"}, + "TD_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "°C"}, + "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "°C"}, + "PMSL.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "Pa"}, + "PS.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "Pa"}, + "TOT_PREC.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "mm"}, + + # Bias: + "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "m/s"}, + "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "m/s"}, + "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "°C"}, + "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "°C"}, + "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "Pa"}, + "PS.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "Pa"}, + "TOT_PREC.BIAS.spatial": {"cmap": plt.get_cmap("BrBG", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "mm"} } CMAP_DEFAULTS = defaultdict(_fallback, _CMAP_DEFAULTS) From 27b91e2d123af3a1b142f73e9cf65266db8ec91e Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Wed, 14 Jan 2026 17:55:36 +0100 Subject: [PATCH 14/58] Better comments in the colour map code. --- src/plotting/colormap_defaults.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index 4a85f67e..5e22f72f 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -57,7 +57,8 @@ def _norm(data): }, # hard-code this for the moment, can still make smarter later on: - # RMSE and MAE first (is all the same). Use Reds colormap to indicate 'error'. + # RMSE and MAE first (is all the same). Sequential colour map to reflect the nature of the data (error, all positive). + # Red is suggestive of 'bad' (high error). # Use a limited number of levels so that absolute values of error can be read from the map. # always start at 0 so that the saturation of the colour corresponds to the error magnitude. @@ -80,6 +81,8 @@ def _norm(data): "TOT_PREC.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "mm"}, # Bias: + # diverging colour scheme for the Bias to reflect the nature of the data (can be positive or negative, symmetric). + # Red-Blue colour scheme for all variables except precipitation, where a Brown-Green scheme is more suggestive. "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "m/s"}, "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "m/s"}, "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "°C"}, From 1b2e6702e4bdf4f976da32d5b96981e372073f65 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Wed, 14 Jan 2026 17:56:39 +0100 Subject: [PATCH 15/58] Better Comments, some further changes to code. --- workflow/scripts/plot_summary_stat_maps.mo.py | 25 ++++++++----------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index 4b8eb38f..14eb182a 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -23,13 +23,12 @@ def _(): # no changes to StatePlotter required according to ChatGPT. from plotting import StatePlotter - # Probably need some new colour maps. - # at least one for biases (diverging), maybe different diverging ones for the different variables. - # + # Added some new colour maps for the Bias / MAE / RMSE map plots. from plotting.colormap_defaults import CMAP_DEFAULTS - # need to load nc files... - from plotting.compat import load_state_from_grib # TODO: load state from nc? + # need to load nc files. But this statement is not needed any more because + # the .nc files can just be read with xr.open_dataset + # from plotting.compat import load_state_from_grib return ( ArgumentParser, @@ -37,7 +36,7 @@ def _(): Path, StatePlotter, ekp, - load_state_from_grib, + # load_state_from_grib, logging, np, DOMAINS, @@ -58,25 +57,23 @@ def _(ArgumentParser, Path): parser = ArgumentParser() parser.add_argument( - "--input", type=str, default=None, help="Directory to grib data" + "--input", type=str, default=None, help="Directory to .nc data containing the error fields" ) - parser.add_argument("--date", type=str, default=None, help="reference datetime") + parser.add_argument("--outfn", type=str, help="output filename") parser.add_argument("--leadtime", type=str, help="leadtime") parser.add_argument("--param", type=str, help="parameter") parser.add_argument("--region", type=str, help="name of region") args = parser.parse_args() - grib_dir = Path(args.input) - init_time = args.date + nc_dir = Path(args.input) outfn = Path(args.outfn) lead_time = args.leadtime param = args.param region = args.region return ( args, - grib_dir, - init_time, + nc_dir, lead_time, outfn, param, @@ -85,9 +82,9 @@ def _(ArgumentParser, Path): @app.cell -def _(grib_dir, init_time, lead_time, load_state_from_grib, param): +def _(nc_dir, init_time, lead_time, load_state_from_grib, param): # load grib file - grib_file = grib_dir / f"{init_time}_{lead_time}.grib" + grib_file = nc_dir / f"{init_time}_{lead_time}.grib" if param == "SP_10M": paramlist = ["U_10M", "V_10M"] elif param == "SP": From 52dd1e71e7a84f00e8197e2b68a78176270691c3 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 15 Jan 2026 15:46:47 +0100 Subject: [PATCH 16/58] Added back instances of lead time. --- workflow/scripts/plot_summary_stat_maps.mo.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index 14eb182a..2a60de34 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -59,7 +59,7 @@ def _(ArgumentParser, Path): parser.add_argument( "--input", type=str, default=None, help="Directory to .nc data containing the error fields" ) - + # parser.add_argument("--date", type=str, default=None, help="reference datetime") # to be deleted? parser.add_argument("--outfn", type=str, help="output filename") parser.add_argument("--leadtime", type=str, help="leadtime") parser.add_argument("--param", type=str, help="parameter") @@ -67,6 +67,7 @@ def _(ArgumentParser, Path): args = parser.parse_args() nc_dir = Path(args.input) + # init_time = args.date # to be deleted? outfn = Path(args.outfn) lead_time = args.leadtime param = args.param @@ -74,6 +75,7 @@ def _(ArgumentParser, Path): return ( args, nc_dir, + # init_time, # to be deleted? lead_time, outfn, param, @@ -83,8 +85,8 @@ def _(ArgumentParser, Path): @app.cell def _(nc_dir, init_time, lead_time, load_state_from_grib, param): - # load grib file - grib_file = nc_dir / f"{init_time}_{lead_time}.grib" + # load .nc verification file + nc_file = nc_dir / "verif_aggregated.nc" if param == "SP_10M": paramlist = ["U_10M", "V_10M"] elif param == "SP": From ddc5883ce3e953af50a93ef10eb7d53f3c349e3e Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 15 Jan 2026 17:55:59 +0100 Subject: [PATCH 17/58] New function for loading netCDF files added to compat.py --- src/plotting/compat.py | 76 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) diff --git a/src/plotting/compat.py b/src/plotting/compat.py index 665287e0..fbe26f3e 100644 --- a/src/plotting/compat.py +++ b/src/plotting/compat.py @@ -90,3 +90,79 @@ def load_state_from_raw( if key.startswith("field_"): state["fields"][key.removeprefix("field_")] = value return state + + +def load_state_from_netcdf( + file: Path, + paramlist: list[str], + *, + season: str = "all", + init_hour: int = -999, + lead_time: int | None = None, # hours +) -> dict: + """ + NetCDF analogue of load_state_from_grib(), restricted to spatial variables. + """ + + ds = xr.open_dataset(file) + + # --- normalize lead_time to hours (float) --- + if ds["lead_time"].dtype.kind == "m": + ds = ds.assign_coords( + lead_time=ds["lead_time"].dt.total_seconds() / 3600 + ) + + # --- select season / init_hour / lead_time --- + ds = ds.sel(season=season, init_hour=init_hour) + if lead_time is not None: + ds = ds.sel(lead_time=lead_time) + + # --- infer reference + valid time --- + # Assumption: forecast_reference_time is not explicitly stored + # We reconstruct something consistent with GRIB usage + forecast_reference_time = None + valid_time = None + if lead_time is not None: + valid_time = pd.to_datetime(lead_time, unit="h", origin="unix") + + # --- get lat / lon (assumed present as coordinates) --- + lat = ds["lat"].values if "lat" in ds.coords else ds["latitude"].values + lon = ds["lon"].values if "lon" in ds.coords else ds["longitude"].values + + lon2d, lat2d = np.meshgrid(lon, lat) + lats = lat2d.flatten() + lons = lon2d.flatten() + + state = { + "forecast_reference_time": forecast_reference_time, + "valid_time": valid_time, + "latitudes": lats, + "longitudes": lons, + "fields": {}, + } + + # --- LAM envelope (convex hull) --- + lam_hull = MultiPoint(list(zip(lons.tolist(), lats.tolist()))).convex_hull + state["lam_envelope"] = gpd.GeoSeries([lam_hull], crs="EPSG:4326") + + # --- extract spatial fields --- + for param in paramlist: + # e.g. U_10M.MAE.spatial + matching_vars = [ + v for v in ds.data_vars + if v.startswith(f"{param}.") and v.endswith(".spatial") + ] + + if not matching_vars: + state["fields"][param] = np.full(lats.size, np.nan, dtype=float) + continue + + # If multiple metrics exist, concatenate them + arrays = [] + for var in matching_vars: + arr = ds[var].values # (lead_time, y, x) or (y, x) + arrays.append(arr.reshape(-1)) + + state["fields"][param] = np.concatenate(arrays) + + return state From 35396a70e4b0603bd5b1bb8cdf73fa3610816cab Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 15 Jan 2026 18:06:05 +0100 Subject: [PATCH 18/58] Marimo app cell for loading data from .nc --- src/plotting/compat.py | 1 + workflow/scripts/plot_summary_stat_maps.mo.py | 13 +++++++++---- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/plotting/compat.py b/src/plotting/compat.py index fbe26f3e..06609d1d 100644 --- a/src/plotting/compat.py +++ b/src/plotting/compat.py @@ -5,6 +5,7 @@ import geopandas as gpd import numpy as np import pandas as pd +import xarray as xr from meteodatalab import data_source from meteodatalab import grib_decoder from shapely.geometry import MultiPoint diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index 2a60de34..27cc957c 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -84,16 +84,21 @@ def _(ArgumentParser, Path): @app.cell -def _(nc_dir, init_time, lead_time, load_state_from_grib, param): - # load .nc verification file - nc_file = nc_dir / "verif_aggregated.nc" +def _(nc_file, param, lead_time, load_state_from_netcdf): + # load .nc verification file: if param == "SP_10M": paramlist = ["U_10M", "V_10M"] elif param == "SP": paramlist = ["U", "V"] else: paramlist = [param] - state = load_state_from_grib(grib_file, paramlist=paramlist) + + state = load_state_from_netcdf( + nc_file, + paramlist=paramlist, + lead_time=lead_time, + ) + return (state,) From e8a92aaa738e8c5f9e44268057ecf21a0ea4e8d6 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Mon, 19 Jan 2026 15:50:02 +0100 Subject: [PATCH 19/58] Remove .nc-loading function again, do it with earthkit instead. --- src/plotting/compat.py | 75 ------------------------------------------ 1 file changed, 75 deletions(-) diff --git a/src/plotting/compat.py b/src/plotting/compat.py index 06609d1d..93634aa3 100644 --- a/src/plotting/compat.py +++ b/src/plotting/compat.py @@ -92,78 +92,3 @@ def load_state_from_raw( state["fields"][key.removeprefix("field_")] = value return state - -def load_state_from_netcdf( - file: Path, - paramlist: list[str], - *, - season: str = "all", - init_hour: int = -999, - lead_time: int | None = None, # hours -) -> dict: - """ - NetCDF analogue of load_state_from_grib(), restricted to spatial variables. - """ - - ds = xr.open_dataset(file) - - # --- normalize lead_time to hours (float) --- - if ds["lead_time"].dtype.kind == "m": - ds = ds.assign_coords( - lead_time=ds["lead_time"].dt.total_seconds() / 3600 - ) - - # --- select season / init_hour / lead_time --- - ds = ds.sel(season=season, init_hour=init_hour) - if lead_time is not None: - ds = ds.sel(lead_time=lead_time) - - # --- infer reference + valid time --- - # Assumption: forecast_reference_time is not explicitly stored - # We reconstruct something consistent with GRIB usage - forecast_reference_time = None - valid_time = None - if lead_time is not None: - valid_time = pd.to_datetime(lead_time, unit="h", origin="unix") - - # --- get lat / lon (assumed present as coordinates) --- - lat = ds["lat"].values if "lat" in ds.coords else ds["latitude"].values - lon = ds["lon"].values if "lon" in ds.coords else ds["longitude"].values - - lon2d, lat2d = np.meshgrid(lon, lat) - lats = lat2d.flatten() - lons = lon2d.flatten() - - state = { - "forecast_reference_time": forecast_reference_time, - "valid_time": valid_time, - "latitudes": lats, - "longitudes": lons, - "fields": {}, - } - - # --- LAM envelope (convex hull) --- - lam_hull = MultiPoint(list(zip(lons.tolist(), lats.tolist()))).convex_hull - state["lam_envelope"] = gpd.GeoSeries([lam_hull], crs="EPSG:4326") - - # --- extract spatial fields --- - for param in paramlist: - # e.g. U_10M.MAE.spatial - matching_vars = [ - v for v in ds.data_vars - if v.startswith(f"{param}.") and v.endswith(".spatial") - ] - - if not matching_vars: - state["fields"][param] = np.full(lats.size, np.nan, dtype=float) - continue - - # If multiple metrics exist, concatenate them - arrays = [] - for var in matching_vars: - arr = ds[var].values # (lead_time, y, x) or (y, x) - arrays.append(arr.reshape(-1)) - - state["fields"][param] = np.concatenate(arrays) - - return state From 06be6fa853dc3536e6271b0bb3d453bcc6149f69 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Wed, 21 Jan 2026 18:04:30 +0100 Subject: [PATCH 20/58] All kinds of changes. Co-Development session with Francesco. Got a long way towards the png plots. Co-authored-by: Francesco Zanetta --- workflow/Snakefile | 7 + workflow/rules/plot.smk | 19 +- workflow/scripts/plot_summary_stat_maps.mo.py | 205 +++++------------- 3 files changed, 77 insertions(+), 154 deletions(-) diff --git a/workflow/Snakefile b/workflow/Snakefile index 7f430cf9..7a6e08b3 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -52,6 +52,13 @@ rule experiment_all: OUT_ROOT / "data/runs/{run_id}/summary.md", run_id=collect_all_candidates(), ), + expand( + OUT_ROOT / "results/experiment/metrics_spatial/{run_id}/{param}_{metric}_{leadtime}.png", + run_id=collect_all_candidates(), + leadtime=[6, 12], + metric=["BIAS", "RMSE", "MAE"], + param=["TOT_PREC", "T_2M"], + ) rule showcase_all: diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 2014e000..6b827dfc 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -70,11 +70,12 @@ rule make_forecast_animation: rule plot_summary_stat_maps: + localrule: True input: script="workflow/scripts/plot_summary_stat_maps.mo.py", - inference_okfile=rules.execute_inference.output.okfile, + verif_file=OUT_ROOT / "data/runs/{run_id}/verif_aggregated.nc" output: - OUT_ROOT / "results/summary_stats/maps/{run_id}/{leadtime}/{metric}_{param}_{region}.png", + OUT_ROOT / "results/experiment/metrics_spatial/{run_id}/{param}_{metric}_{leadtime}.png", wildcard_constraints: leadtime=r"\d+", # only digits resources: @@ -83,18 +84,18 @@ rule plot_summary_stat_maps: runtime="10m", params: nc_out_dir=lambda wc: ( - Path(OUT_ROOT) / f"data/runs/{wc.run_id}/{wc.init_time}/grib" + Path(OUT_ROOT) / f"data/runs/{wc.run_id}/verif_aggregated.nc" # not sure how to do this, because the baselines are in, e.g., output/data/baselines/COSMO-E/verif_aggregated.nc # and the runs are in output/data/runs/runID/verif_aggregated.nc ).resolve(), shell: """ export ECCODES_DEFINITION_PATH=$(realpath .venv/share/eccodes-cosmo-resources/definitions) - python {input.script} \ - --input {params.nc_out_dir} --date {wildcards.init_time} --outfn {output[0]} \ - --param {wildcards.param} --leadtime {wildcards.leadtime} --region {wildcards.region} \ + # python {input.script} \ + # --input {input.verif_file} --outfn {output[0]} \ + # --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} # interactive editing (needs to set localrule: True and use only one core) - # marimo edit {input.script} -- \ - # --input {params.grib_out_dir} --date {wildcards.init_time} --outfn {output[0]}\ - # --param {wildcards.param} --leadtime {wildcards.leadtime} --region {wildcards.region}\ + marimo edit {input.script} -- \ + --input {input.verif_file} --outfn {output[0]}\ + --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} """ \ No newline at end of file diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index 27cc957c..54bf7863 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -1,12 +1,12 @@ import marimo -__generated_with = "0.16.5" +__generated_with = "0.19.4" app = marimo.App(width="medium") @app.cell def _(): - + # this sure stays the same. import logging from argparse import ArgumentParser @@ -16,10 +16,11 @@ def _(): import cartopy.crs as ccrs import earthkit.plots as ekp import numpy as np + import xarray as xr # this stays the same as well. from plotting import DOMAINS - + # no changes to StatePlotter required according to ChatGPT. from plotting import StatePlotter @@ -28,19 +29,16 @@ def _(): # need to load nc files. But this statement is not needed any more because # the .nc files can just be read with xr.open_dataset - # from plotting.compat import load_state_from_grib - return ( ArgumentParser, CMAP_DEFAULTS, + DOMAINS, Path, StatePlotter, ekp, - # load_state_from_grib, logging, np, - DOMAINS, - ccrs, + xr, ) @@ -49,57 +47,51 @@ def _(logging): LOG = logging.getLogger(__name__) LOG_FMT = "%(asctime)s - %(name)s - %(levelname)s - %(message)s" logging.basicConfig(level=logging.INFO, format=LOG_FMT) - return (LOG,) + return @app.cell -def _(ArgumentParser, Path): +def _(ArgumentParser, Path, np): parser = ArgumentParser() parser.add_argument( "--input", type=str, default=None, help="Directory to .nc data containing the error fields" ) - # parser.add_argument("--date", type=str, default=None, help="reference datetime") # to be deleted? parser.add_argument("--outfn", type=str, help="output filename") parser.add_argument("--leadtime", type=str, help="leadtime") parser.add_argument("--param", type=str, help="parameter") - parser.add_argument("--region", type=str, help="name of region") + # parser.add_argument("--region", type=str, help="name of region") + parser.add_argument("--metric", type=str, help = "Evaluation Metric. So far Bias, RMSE or MAE are implemented.") + parser.add_argument("--season", type=str, default="all", help="season filter") + parser.add_argument("--init_hour", type=str, default="all", help="initialization hour filter") args = parser.parse_args() - nc_dir = Path(args.input) - # init_time = args.date # to be deleted? + verif_file = Path(args.input) outfn = Path(args.outfn) lead_time = args.leadtime param = args.param - region = args.region - return ( - args, - nc_dir, - # init_time, # to be deleted? - lead_time, - outfn, - param, - region, - ) + # region = args.region + season = args.season + init_hour = args.init_hour + metric = args.metric + if isinstance(init_hour, str): + if init_hour == "all": + init_hour = -999 + else: + raise ValueError("init_hour must be 'all' or an integer hour") -@app.cell -def _(nc_file, param, lead_time, load_state_from_netcdf): - # load .nc verification file: - if param == "SP_10M": - paramlist = ["U_10M", "V_10M"] - elif param == "SP": - paramlist = ["U", "V"] - else: - paramlist = [param] - - state = load_state_from_netcdf( - nc_file, - paramlist=paramlist, - lead_time=lead_time, - ) + lead_time = np.timedelta64(lead_time, 'h') + return init_hour, lead_time, metric, outfn, param, season, verif_file - return (state,) + +@app.cell +def _(init_hour, lead_time, metric, param, season, verif_file, xr): + ds = xr.open_dataset(verif_file) + var = f"{param}.{metric}.spatial" + ds = ds[var].sel(init_hour=init_hour, lead_time=lead_time, season=season) + ds + return ds, var @app.cell @@ -125,129 +117,52 @@ def get_style(param, units_override=None): "vmax": cfg.get("vmax", None), "colors": cfg.get("colors", None), } - return (get_style,) @app.cell -def _(LOG, np): - """Preprocess fields with pint-based unit conversion and derived quantities.""" - try: - import pint # type: ignore - - _ureg = pint.UnitRegistry() - - def _k_to_c(arr): - # robust conversion with pint, fallback if dtype unsupported - try: - return (_ureg.Quantity(arr, _ureg.kelvin).to(_ureg.degC)).magnitude - except Exception: - return arr - 273.15 - - def _ms_to_knots(arr): - # robust conversion with pint, fallback if dtype unsupported - try: - return ( - _ureg.Quantity(arr, _ureg.meter / _ureg.second).to(_ureg.knot) - ).magnitude - except Exception: - return arr * 1.943844 - - def _m_to_mm(arr): - # robust conversion with pint, fallback if dtype unsupported - try: - return (_ureg.Quantity(arr, _ureg.meter).to(_ureg.millimeter)).magnitude - except Exception: - return arr * 1000 - - except Exception: - LOG.warning("pint not available; falling back hardcoded conversions") - - def _k_to_c(arr): - return arr - 273.15 - - def _ms_to_knots(arr): - return arr * 1.943844 - - def _m_to_mm(arr): - return arr * 1000 - - 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"] - # temperature variables - if param in ("T_2M", "TD_2M", "T", "TD"): - return _k_to_c(fields[param]), "°C" - # 10m wind speed (allow legacy 'uv' alias) - if param == "SP_10M": - u = fields["U_10M"] - v = fields["V_10M"] - return np.sqrt(u**2 + v**2), "m/s" - # wind speed from standard-level components - if param == "SP": - u = fields["U"] - v = fields["V"] - return np.sqrt(u**2 + v**2), "m/s" - if param == "TOT_PREC": - return _m_to_mm(fields[param]), "mm" - # default: passthrough - return fields[param], None - - return (preprocess_field,) - - -@app.cell -def _( - LOG, - StatePlotter, - args, - get_style, - outfn, - param, - preprocess_field, - region, - state, - DOMAINS, - ccrs, -): +def _(DOMAINS, StatePlotter, ds, get_style, outfn, var): # plot individual fields + import matplotlib.pyplot as plt + plotter = StatePlotter( - state["longitudes"], - state["latitudes"], + ds["longitude"].values.ravel(), + ds["latitude"].values.ravel(), outfn.parent, ) fig = plotter.init_geoaxes( nrows=1, ncols=1, - projection=DOMAINS[region]["projection"], - bbox=DOMAINS[region]["extent"], - name=region, + projection=DOMAINS["switzerland"]["projection"], + bbox=DOMAINS["switzerland"]["extent"], + name="switzerland", size=(6, 6), ) subplot = fig.add_map(row=0, column=0) - # preprocess field (unit conversion, derived quantities) - field, units_override = preprocess_field(param, state) + # # preprocess field (unit conversion, derived quantities) + # field, units_override = preprocess_field(param, state) - plotter.plot_field(subplot, field, **get_style(args.param, units_override)) - subplot.ax.add_geometries( - state["lam_envelope"], - edgecolor="black", - facecolor="none", - crs=ccrs.PlateCarree(), - ) - validtime = state["valid_time"].strftime("%Y%m%d%H%M") - # leadtime = int(state["lead_time"].total_seconds() // 3600) + plotter.plot_field(subplot, ds.values.ravel(), **get_style(var)) + # subplot.ax.add_geometries( + # state["lam_envelope"], + # edgecolor="black", + # facecolor="none", + # crs=ccrs.PlateCarree(), + # ) + plt.show() + # validtime = state["valid_time"].strftime("%Y%m%d%H%M") + # # leadtime = int(state["lead_time"].total_seconds() // 3600) + + # fig.title(f"{param}, time: {validtime}") + + # fig.save(outfn, bbox_inches="tight", dpi=200) + # LOG.info(f"saved: {outfn}") + return - fig.title(f"{param}, time: {validtime}") - fig.save(outfn, bbox_inches="tight", dpi=200) - LOG.info(f"saved: {outfn}") +@app.cell +def _(): return From 61c728ec66ea55bfbcbe7164920d44115d590fa0 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 22 Jan 2026 17:14:21 +0100 Subject: [PATCH 21/58] Generalized to the other non-trivial (non-wind) variables. --- workflow/Snakefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/workflow/Snakefile b/workflow/Snakefile index 7a6e08b3..41b54415 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -57,7 +57,7 @@ rule experiment_all: run_id=collect_all_candidates(), leadtime=[6, 12], metric=["BIAS", "RMSE", "MAE"], - param=["TOT_PREC", "T_2M"], + param=["TD_2M", "T_2M", "PMSL", "PS", "TOT_PREC"], ) From 761f8d94cff59a721172260f54f5781d4ea64e37 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 22 Jan 2026 17:24:11 +0100 Subject: [PATCH 22/58] Some changes to plotting script. --- workflow/scripts/plot_summary_stat_maps.mo.py | 22 +++++++++++++------ 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index 54bf7863..52ea116d 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -26,9 +26,6 @@ def _(): # Added some new colour maps for the Bias / MAE / RMSE map plots. from plotting.colormap_defaults import CMAP_DEFAULTS - - # need to load nc files. But this statement is not needed any more because - # the .nc files can just be read with xr.open_dataset return ( ArgumentParser, CMAP_DEFAULTS, @@ -47,7 +44,7 @@ def _(logging): LOG = logging.getLogger(__name__) LOG_FMT = "%(asctime)s - %(name)s - %(levelname)s - %(message)s" logging.basicConfig(level=logging.INFO, format=LOG_FMT) - return + return (LOG,) @app.cell @@ -121,7 +118,18 @@ def get_style(param, units_override=None): @app.cell -def _(DOMAINS, StatePlotter, ds, get_style, outfn, var): +def _( + DOMAINS, + LOG, + StatePlotter, + ds, + get_style, + lead_time, + metric, + outfn, + param, + var, +): # plot individual fields import matplotlib.pyplot as plt @@ -150,11 +158,11 @@ def _(DOMAINS, StatePlotter, ds, get_style, outfn, var): # facecolor="none", # crs=ccrs.PlateCarree(), # ) - plt.show() + # validtime = state["valid_time"].strftime("%Y%m%d%H%M") # # leadtime = int(state["lead_time"].total_seconds() // 3600) - # fig.title(f"{param}, time: {validtime}") + fig.title(f"{metric} of {param}, Lead Time: {lead_time}") # fig.save(outfn, bbox_inches="tight", dpi=200) # LOG.info(f"saved: {outfn}") From 1cd4dbf22e87464b69cd52d7126109d5722abd8d Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Fri, 23 Jan 2026 13:25:07 +0100 Subject: [PATCH 23/58] Plotting region now dynamical (but not yet properly working). Output written to .png now working. --- workflow/Snakefile | 3 ++- workflow/rules/plot.smk | 19 +++++++++++-------- workflow/scripts/plot_summary_stat_maps.mo.py | 13 +++++++------ 3 files changed, 20 insertions(+), 15 deletions(-) diff --git a/workflow/Snakefile b/workflow/Snakefile index 41b54415..0c96d3c8 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -53,11 +53,12 @@ rule experiment_all: run_id=collect_all_candidates(), ), expand( - OUT_ROOT / "results/experiment/metrics_spatial/{run_id}/{param}_{metric}_{leadtime}.png", + OUT_ROOT / "results/experiment/metrics_spatial/{run_id}/{param}_{metric}_{leadtime}_{region}.png", run_id=collect_all_candidates(), leadtime=[6, 12], metric=["BIAS", "RMSE", "MAE"], param=["TD_2M", "T_2M", "PMSL", "PS", "TOT_PREC"], + region=["europe", "switzerland"] ) diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 6b827dfc..cad66bd5 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -73,9 +73,12 @@ rule plot_summary_stat_maps: localrule: True input: script="workflow/scripts/plot_summary_stat_maps.mo.py", - verif_file=OUT_ROOT / "data/runs/{run_id}/verif_aggregated.nc" + verif_file=OUT_ROOT / "data/runs/{run_id}/verif_aggregated.nc", + # verif_file=EXPERIMENT_PARTICIPANTS.values(), + # copied from rule report_experiment_dashboard - should be correct here, + # but needs adjustments in the output as well - tbd. output: - OUT_ROOT / "results/experiment/metrics_spatial/{run_id}/{param}_{metric}_{leadtime}.png", + OUT_ROOT / "results/experiment/metrics_spatial/{run_id}/{param}_{metric}_{leadtime}_{region}.png", wildcard_constraints: leadtime=r"\d+", # only digits resources: @@ -91,11 +94,11 @@ rule plot_summary_stat_maps: shell: """ export ECCODES_DEFINITION_PATH=$(realpath .venv/share/eccodes-cosmo-resources/definitions) - # python {input.script} \ - # --input {input.verif_file} --outfn {output[0]} \ - # --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} + python {input.script} \ + --input {input.verif_file} --outfn {output[0]} --region {wildcards.region} \ + --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} \ # interactive editing (needs to set localrule: True and use only one core) - marimo edit {input.script} -- \ - --input {input.verif_file} --outfn {output[0]}\ - --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} + # marimo edit {input.script} -- \ + # --input {input.verif_file} --outfn {output[0]} --region {wildcards.region} \ + # --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} \ """ \ No newline at end of file diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index 52ea116d..0f68c203 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -57,7 +57,7 @@ def _(ArgumentParser, Path, np): parser.add_argument("--outfn", type=str, help="output filename") parser.add_argument("--leadtime", type=str, help="leadtime") parser.add_argument("--param", type=str, help="parameter") - # parser.add_argument("--region", type=str, help="name of region") + parser.add_argument("--region", type=str, help="name of region") parser.add_argument("--metric", type=str, help = "Evaluation Metric. So far Bias, RMSE or MAE are implemented.") parser.add_argument("--season", type=str, default="all", help="season filter") parser.add_argument("--init_hour", type=str, default="all", help="initialization hour filter") @@ -128,6 +128,7 @@ def _( metric, outfn, param, + region, var, ): # plot individual fields @@ -141,9 +142,9 @@ def _( fig = plotter.init_geoaxes( nrows=1, ncols=1, - projection=DOMAINS["switzerland"]["projection"], - bbox=DOMAINS["switzerland"]["extent"], - name="switzerland", + projection=DOMAINS[region]["projection"], + bbox=DOMAINS[region]["extent"], + name=region, size=(6, 6), ) subplot = fig.add_map(row=0, column=0) @@ -164,8 +165,8 @@ def _( fig.title(f"{metric} of {param}, Lead Time: {lead_time}") - # fig.save(outfn, bbox_inches="tight", dpi=200) - # LOG.info(f"saved: {outfn}") + fig.save(outfn, bbox_inches="tight", dpi=200) + LOG.info(f"saved: {outfn}") return From cdb2ccd8f91890b55386c9dbaabfbbd1b4eb967b Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Fri, 23 Jan 2026 13:31:00 +0100 Subject: [PATCH 24/58] Dynamic Regions now working. --- workflow/scripts/plot_summary_stat_maps.mo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index 0f68c203..a2be8ca1 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -67,7 +67,7 @@ def _(ArgumentParser, Path, np): outfn = Path(args.outfn) lead_time = args.leadtime param = args.param - # region = args.region + region = args.region season = args.season init_hour = args.init_hour metric = args.metric From 558689c980b886c53d4a99022c477c76e7f14149 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Fri, 23 Jan 2026 14:36:07 +0100 Subject: [PATCH 25/58] Store results under experiment hash. --- workflow/Snakefile | 5 +++-- workflow/rules/plot.smk | 5 +---- workflow/scripts/plot_summary_stat_maps.mo.py | 11 ++++++++++- 3 files changed, 14 insertions(+), 7 deletions(-) diff --git a/workflow/Snakefile b/workflow/Snakefile index 0c96d3c8..f8c9214f 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -53,12 +53,13 @@ rule experiment_all: run_id=collect_all_candidates(), ), expand( - OUT_ROOT / "results/experiment/metrics_spatial/{run_id}/{param}_{metric}_{leadtime}_{region}.png", + OUT_ROOT / "results/{experiment}/metrics/spatial/{run_id}/{param}_{metric}_{leadtime}_{region}.png", run_id=collect_all_candidates(), leadtime=[6, 12], metric=["BIAS", "RMSE", "MAE"], param=["TD_2M", "T_2M", "PMSL", "PS", "TOT_PREC"], - region=["europe", "switzerland"] + region=["europe", "switzerland"], + experiment=EXPERIMENT_HASH ) diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index cad66bd5..298784a2 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -74,11 +74,8 @@ rule plot_summary_stat_maps: input: script="workflow/scripts/plot_summary_stat_maps.mo.py", verif_file=OUT_ROOT / "data/runs/{run_id}/verif_aggregated.nc", - # verif_file=EXPERIMENT_PARTICIPANTS.values(), - # copied from rule report_experiment_dashboard - should be correct here, - # but needs adjustments in the output as well - tbd. output: - OUT_ROOT / "results/experiment/metrics_spatial/{run_id}/{param}_{metric}_{leadtime}_{region}.png", + OUT_ROOT / "results/{experiment}/metrics/spatial/{run_id}/{param}_{metric}_{leadtime}_{region}.png", wildcard_constraints: leadtime=r"\d+", # only digits resources: diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index a2be8ca1..49ce4893 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -79,7 +79,16 @@ def _(ArgumentParser, Path, np): raise ValueError("init_hour must be 'all' or an integer hour") lead_time = np.timedelta64(lead_time, 'h') - return init_hour, lead_time, metric, outfn, param, season, verif_file + return ( + init_hour, + lead_time, + metric, + outfn, + param, + region, + season, + verif_file, + ) @app.cell From 986a7ee0ce75ba9bfa1f91407c83de7d579dea80 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Mon, 26 Jan 2026 15:05:39 +0100 Subject: [PATCH 26/58] Introduced new domain "switzerland_small" for more detailed inspection of results at smaller spatial scale. --- src/plotting/__init__.py | 4 ++++ workflow/Snakefile | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/plotting/__init__.py b/src/plotting/__init__.py index ce5e4e63..520e35ae 100644 --- a/src/plotting/__init__.py +++ b/src/plotting/__init__.py @@ -37,6 +37,10 @@ "extent": [0, 17.5, 40.5, 53.0], "projection": _PROJECTIONS["orthographic"], }, + "switzerland_small": { + "extent": [5.5, 11.0, 45.5, 48.0], + "projection": _PROJECTIONS["orthographic"], + }, } diff --git a/workflow/Snakefile b/workflow/Snakefile index f8c9214f..64a30694 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -58,7 +58,7 @@ rule experiment_all: leadtime=[6, 12], metric=["BIAS", "RMSE", "MAE"], param=["TD_2M", "T_2M", "PMSL", "PS", "TOT_PREC"], - region=["europe", "switzerland"], + region=["switzerland", "switzerland_small"], experiment=EXPERIMENT_HASH ) From 7f70fb453c49340e377e8224ccbeb1fce08b3ce6 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Tue, 27 Jan 2026 16:04:46 +0100 Subject: [PATCH 27/58] Reverse Red-Blue colour maps for bias. --- src/plotting/colormap_defaults.py | 34 +++++++++++++++---------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index 5e22f72f..5539b363 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -11,16 +11,16 @@ def _fallback(): warnings.warn("No colormap found for this parameter, using fallback.", UserWarning) return {"cmap": plt.get_cmap("viridis"), "norm": None, "units": ""} -def symmetric_boundary_norm(nlevels): - """ - Returns a callable that creates a symmetric BoundaryNorm - around zero with `nlevels` discrete colors. Used for creating colormaps for bias. - """ - def _norm(data): - vmax = np.nanmax(np.abs(data)) - boundaries = np.linspace(-vmax, vmax, nlevels + 1) - return BoundaryNorm(boundaries=boundaries, ncolors=nlevels) - return _norm +# def symmetric_boundary_norm(nlevels): +# """ +# Returns a callable that creates a symmetric BoundaryNorm +# around zero with `nlevels` discrete colors. Used for creating colormaps for bias. +# """ +# def _norm(data): +# vmax = np.nanmax(np.abs(data)) +# boundaries = np.linspace(-vmax, vmax, nlevels + 1) +# return BoundaryNorm(boundaries=boundaries, ncolors=nlevels) +# return _norm _CMAP_DEFAULTS = { "SP": {"cmap": plt.get_cmap("coolwarm", 11), "vmin": 800 * 100, "vmax": 1100 * 100}, @@ -83,13 +83,13 @@ def _norm(data): # Bias: # diverging colour scheme for the Bias to reflect the nature of the data (can be positive or negative, symmetric). # Red-Blue colour scheme for all variables except precipitation, where a Brown-Green scheme is more suggestive. - "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "m/s"}, - "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "m/s"}, - "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "°C"}, - "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "°C"}, - "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "Pa"}, - "PS.BIAS.spatial": {"cmap": plt.get_cmap("RdBu", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "Pa"}, - "TOT_PREC.BIAS.spatial": {"cmap": plt.get_cmap("BrBG", 11), "norm": symmetric_boundary_norm(nlevels=11)} | {"units": "mm"} + "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, + "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, + "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "°C"}, + "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "°C"}, + "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, + "PS.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, + "TOT_PREC.BIAS.spatial": {"cmap": plt.get_cmap("BrBG", 11)} | {"units": "mm"} } CMAP_DEFAULTS = defaultdict(_fallback, _CMAP_DEFAULTS) From 50e899fdce05ffe395b98a0d107ff986cee79fe4 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Tue, 27 Jan 2026 16:08:51 +0100 Subject: [PATCH 28/58] Preliminary changes to plotting script for getting symmetric colour map for bias. --- workflow/scripts/plot_summary_stat_maps.mo.py | 25 ++++++++++++++----- 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index 49ce4893..2d81ac8d 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -102,25 +102,38 @@ def _(init_hour, lead_time, metric, param, season, verif_file, xr): @app.cell def _(CMAP_DEFAULTS, ekp): - def get_style(param, units_override=None): + def get_style(param, metric, units_override=None): """Get style and colormap settings for the plot. Needed because cmap/norm does not work in Style(colors=cmap), still needs to be passed as arguments to tripcolor()/tricontourf(). """ cfg = CMAP_DEFAULTS[param] units = units_override if units_override is not None else cfg.get("units", "") + + levels = cfg.get("levels", None) + + if (metric == "BIAS"): + # For bias, we want to use a diverging colormap with symmetric bounds around zero. + max_abs_val = max(abs(levels)) + levels = numpy.arange(-max_abs_val, max_abs_val, dtype=None) + vmin = -max_abs_val + vmax = max_abs_val + else: + vmin = cfg.get("vmin", None) + vmax = cfg.get("vmax", None) + return { "style": ekp.styles.Style( - levels=cfg.get("bounds", cfg.get("levels", None)), + levels=cfg.get("bounds", levels), extend="both", units=units, colors=cfg.get("colors", None), ), "norm": cfg.get("norm", None), "cmap": cfg.get("cmap", None), - "levels": cfg.get("levels", None), - "vmin": cfg.get("vmin", None), - "vmax": cfg.get("vmax", None), + "levels": levels, + "vmin": vmin, + "vmax": vmax, "colors": cfg.get("colors", None), } return (get_style,) @@ -161,7 +174,7 @@ def _( # # preprocess field (unit conversion, derived quantities) # field, units_override = preprocess_field(param, state) - plotter.plot_field(subplot, ds.values.ravel(), **get_style(var)) + plotter.plot_field(subplot, ds.values.ravel(), **get_style(var, metric)) # subplot.ax.add_geometries( # state["lam_envelope"], # edgecolor="black", From a0aefc0ad065e272b4cbc28b41fd21de30f1eb1b Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Tue, 27 Jan 2026 16:23:57 +0100 Subject: [PATCH 29/58] Temporarily changed plotting script back to original to see if all of it still works. --- workflow/scripts/plot_summary_stat_maps.mo.py | 62 +++++++++++++------ 1 file changed, 43 insertions(+), 19 deletions(-) diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index 2d81ac8d..9b51c920 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -99,45 +99,69 @@ def _(init_hour, lead_time, metric, param, season, verif_file, xr): ds return ds, var - @app.cell def _(CMAP_DEFAULTS, ekp): - def get_style(param, metric, units_override=None): + def get_style(param, units_override=None): """Get style and colormap settings for the plot. Needed because cmap/norm does not work in Style(colors=cmap), still needs to be passed as arguments to tripcolor()/tricontourf(). """ cfg = CMAP_DEFAULTS[param] units = units_override if units_override is not None else cfg.get("units", "") - - levels = cfg.get("levels", None) - - if (metric == "BIAS"): - # For bias, we want to use a diverging colormap with symmetric bounds around zero. - max_abs_val = max(abs(levels)) - levels = numpy.arange(-max_abs_val, max_abs_val, dtype=None) - vmin = -max_abs_val - vmax = max_abs_val - else: - vmin = cfg.get("vmin", None) - vmax = cfg.get("vmax", None) - return { "style": ekp.styles.Style( - levels=cfg.get("bounds", levels), + levels=cfg.get("bounds", cfg.get("levels", None)), extend="both", units=units, colors=cfg.get("colors", None), ), "norm": cfg.get("norm", None), "cmap": cfg.get("cmap", None), - "levels": levels, - "vmin": vmin, - "vmax": vmax, + "levels": cfg.get("levels", None), + "vmin": cfg.get("vmin", None), + "vmax": cfg.get("vmax", None), "colors": cfg.get("colors", None), } + return (get_style,) +# @app.cell +# def _(CMAP_DEFAULTS, ekp): +# def get_style(param, metric, units_override=None): +# """Get style and colormap settings for the plot. +# Needed because cmap/norm does not work in Style(colors=cmap), +# still needs to be passed as arguments to tripcolor()/tricontourf(). +# """ +# cfg = CMAP_DEFAULTS[param] +# units = units_override if units_override is not None else cfg.get("units", "") + +# levels = cfg.get("levels", None) + +# if (metric == "BIAS"): +# # For bias, we want to use a diverging colormap with symmetric bounds around zero. +# max_abs_val = max(abs(levels)) +# levels = numpy.arange(-max_abs_val, max_abs_val, dtype=None) +# vmin = -max_abs_val +# vmax = max_abs_val +# else: +# vmin = cfg.get("vmin", None) +# vmax = cfg.get("vmax", None) + +# return { +# "style": ekp.styles.Style( +# levels=cfg.get("bounds", levels), +# extend="both", +# units=units, +# colors=cfg.get("colors", None), +# ), +# "norm": cfg.get("norm", None), +# "cmap": cfg.get("cmap", None), +# "levels": levels, +# "vmin": vmin, +# "vmax": vmax, +# "colors": cfg.get("colors", None), +# } +# return (get_style,) @app.cell def _( From be6fa354094d9a1c8dc88da10d94b64603b26f7a Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Tue, 27 Jan 2026 18:04:52 +0100 Subject: [PATCH 30/58] Fix to the functioning of _compute_scores and _compute_statistics with argument dim. --- src/verification/__init__.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/verification/__init__.py b/src/verification/__init__.py index 8f05349c..747d3b4e 100644 --- a/src/verification/__init__.py +++ b/src/verification/__init__.py @@ -77,12 +77,13 @@ def _compute_scores( prefix="", suffix="", source="", + dim=["x", "y"], ) -> xr.Dataset: """ Compute basic verification metrics between two xarray DataArrays (fcst and obs). Returns a xarray Dataset with the computed metrics. """ - dim = ["x", "y"] if "x" in fcst.dims and "y" in fcst.dims else ["values"] + error = fcst - obs if dim == []: scores = xr.Dataset( @@ -112,12 +113,13 @@ def _compute_statistics( prefix="", suffix="", source="", + dim=["x", "y"], ) -> xr.Dataset: """ Compute basic statistics of a xarray DataArray (data). Returns a xarray Dataset with the computed statistics. """ - dim = ["x", "y"] if "x" in data.dims and "y" in data.dims else ["values"] + stats = xr.Dataset( { f"{prefix}mean{suffix}": data.mean(dim=dim, skipna=True), @@ -161,6 +163,8 @@ def verify( """ start = time.time() + dim = ["x", "y"] if "x" in fcst.dims and "y" in fcst.dims else ["values"] + # rewrite the verification to use dask and xarray # chunk the data to avoid memory issues # compute the metrics in parallel @@ -188,19 +192,19 @@ def verify( # scores vs time (reduce spatially) score.append( _compute_scores( - fcst_param, obs_param, prefix=param + ".", source=fcst_label + fcst_param, obs_param, prefix=param + ".", source=fcst_label, dim=dim ).expand_dims(region=[region]) ) # statistics vs time (reduce spatially) fcst_statistics.append( _compute_statistics( - fcst_param, prefix=param + ".", source=fcst_label + fcst_param, prefix=param + ".", source=fcst_label, dim=dim ).expand_dims(region=[region]) ) obs_statistics.append( _compute_statistics( - obs_param, prefix=param + ".", source=obs_label + obs_param, prefix=param + ".", source=obs_label, dim=dim ).expand_dims(region=[region]) ) From 2fcd1f92a0ecd0ae8de2644ea1b9d3b04a90713a Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 29 Jan 2026 15:32:25 +0100 Subject: [PATCH 31/58] Working version for colour scale for bias that is symmetric about zero. vmin and vmax values for variables other than T_2M yet to be defined. --- src/plotting/colormap_defaults.py | 13 +--- workflow/scripts/plot_summary_stat_maps.mo.py | 63 ++++++------------- 2 files changed, 19 insertions(+), 57 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index 5539b363..10c461a4 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -11,17 +11,6 @@ def _fallback(): warnings.warn("No colormap found for this parameter, using fallback.", UserWarning) return {"cmap": plt.get_cmap("viridis"), "norm": None, "units": ""} -# def symmetric_boundary_norm(nlevels): -# """ -# Returns a callable that creates a symmetric BoundaryNorm -# around zero with `nlevels` discrete colors. Used for creating colormaps for bias. -# """ -# def _norm(data): -# vmax = np.nanmax(np.abs(data)) -# boundaries = np.linspace(-vmax, vmax, nlevels + 1) -# return BoundaryNorm(boundaries=boundaries, ncolors=nlevels) -# return _norm - _CMAP_DEFAULTS = { "SP": {"cmap": plt.get_cmap("coolwarm", 11), "vmin": 800 * 100, "vmax": 1100 * 100}, "TD_2M": load_ncl_colormap("t2m_29lev.ct"), @@ -86,7 +75,7 @@ def _fallback(): "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "°C"}, - "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "°C"}, + "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "vmin": -5.5, "vmax": 5.5} | {"units": "°C"}, "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "PS.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "TOT_PREC.BIAS.spatial": {"cmap": plt.get_cmap("BrBG", 11)} | {"units": "mm"} diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index 9b51c920..349ea300 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -100,68 +100,41 @@ def _(init_hour, lead_time, metric, param, season, verif_file, xr): return ds, var @app.cell -def _(CMAP_DEFAULTS, ekp): - def get_style(param, units_override=None): +def _(CMAP_DEFAULTS, ekp, np): + def get_style(param, metric, units_override=None): """Get style and colormap settings for the plot. Needed because cmap/norm does not work in Style(colors=cmap), still needs to be passed as arguments to tripcolor()/tricontourf(). """ - cfg = CMAP_DEFAULTS[param] + + metric_key = f"{param}.{metric}.spatial" + cfg = CMAP_DEFAULTS[metric_key] if metric_key in CMAP_DEFAULTS else CMAP_DEFAULTS[param] units = units_override if units_override is not None else cfg.get("units", "") + + vmin = cfg.get("vmin", None) + vmax = cfg.get("vmax", None) + levels = cfg.get("levels", None) + + # For the case of Bias, construct symmetric levels: + if metric == "BIAS": + levels = np.linspace(start = vmin, stop = vmax, num=12) + return { "style": ekp.styles.Style( - levels=cfg.get("bounds", cfg.get("levels", None)), + levels=cfg.get("bounds", levels), extend="both", units=units, colors=cfg.get("colors", None), ), "norm": cfg.get("norm", None), "cmap": cfg.get("cmap", None), - "levels": cfg.get("levels", None), - "vmin": cfg.get("vmin", None), - "vmax": cfg.get("vmax", None), + "levels": levels, + "vmin": vmin, + "vmax": vmax, "colors": cfg.get("colors", None), } - return (get_style,) -# @app.cell -# def _(CMAP_DEFAULTS, ekp): -# def get_style(param, metric, units_override=None): -# """Get style and colormap settings for the plot. -# Needed because cmap/norm does not work in Style(colors=cmap), -# still needs to be passed as arguments to tripcolor()/tricontourf(). -# """ -# cfg = CMAP_DEFAULTS[param] -# units = units_override if units_override is not None else cfg.get("units", "") - -# levels = cfg.get("levels", None) - -# if (metric == "BIAS"): -# # For bias, we want to use a diverging colormap with symmetric bounds around zero. -# max_abs_val = max(abs(levels)) -# levels = numpy.arange(-max_abs_val, max_abs_val, dtype=None) -# vmin = -max_abs_val -# vmax = max_abs_val -# else: -# vmin = cfg.get("vmin", None) -# vmax = cfg.get("vmax", None) - -# return { -# "style": ekp.styles.Style( -# levels=cfg.get("bounds", levels), -# extend="both", -# units=units, -# colors=cfg.get("colors", None), -# ), -# "norm": cfg.get("norm", None), -# "cmap": cfg.get("cmap", None), -# "levels": levels, -# "vmin": vmin, -# "vmax": vmax, -# "colors": cfg.get("colors", None), -# } -# return (get_style,) @app.cell def _( From 7864b5d3d68a1541be11bcb951d1856b20004616 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 29 Jan 2026 17:21:22 +0100 Subject: [PATCH 32/58] New colour map defaults for T_2M --- src/plotting/colormap_defaults.py | 28 +++++++++---------- workflow/scripts/plot_summary_stat_maps.mo.py | 4 +++ 2 files changed, 18 insertions(+), 14 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index 10c461a4..e3ed2e22 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -52,22 +52,22 @@ def _fallback(): # always start at 0 so that the saturation of the colour corresponds to the error magnitude. # RMSE: - "U_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "m/s"}, - "V_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "m/s"}, - "TD_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "°C"}, - "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "°C"}, - "PMSL.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "Pa"}, - "PS.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "Pa"}, - "TOT_PREC.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "mm"}, + "U_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, + "V_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, + "TD_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "°C"}, + "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0, "vmax": 6} | {"units": "°C"}, + "PMSL.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, + "PS.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, + "TOT_PREC.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "mm"}, # MAE: - "U_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "m/s"}, - "V_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "m/s"}, - "TD_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "°C"}, - "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "°C"}, - "PMSL.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "Pa"}, - "PS.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "Pa"}, - "TOT_PREC.MAE.spatial": {"cmap": plt.get_cmap("Reds", 11), "vmin": 0} | {"units": "mm"}, + "U_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, + "V_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, + "TD_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "°C"}, + "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0, "vmax": 6} | {"units": "°C"}, + "PMSL.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, + "PS.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, + "TOT_PREC.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "mm"}, # Bias: # diverging colour scheme for the Bias to reflect the nature of the data (can be positive or negative, symmetric). diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index 349ea300..af4be700 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -118,6 +118,10 @@ def get_style(param, metric, units_override=None): # For the case of Bias, construct symmetric levels: if metric == "BIAS": levels = np.linspace(start = vmin, stop = vmax, num=12) + elif metric == "RMSE": + levels = np.linspace(start = vmin, stop = vmax, num=7) + elif metric == "MAE": + levels = np.linspace(start = vmin, stop = vmax, num=7) return { "style": ekp.styles.Style( From de0f2794ed845e44d243b54fec44a3fdc5247576 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Tue, 3 Feb 2026 16:26:46 +0100 Subject: [PATCH 33/58] Hard-Code levels for colour breaks. This way, everything is in the same place and can be changed there when needed. Accompanying changes in the plotting script (marimo cell that gets the colour map defaults). --- src/plotting/colormap_defaults.py | 6 ++--- workflow/scripts/plot_summary_stat_maps.mo.py | 24 ++++--------------- 2 files changed, 8 insertions(+), 22 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index e3ed2e22..c9968240 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -55,7 +55,7 @@ def _fallback(): "U_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "V_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "TD_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "°C"}, - "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0, "vmax": 6} | {"units": "°C"}, + "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 1, 2, 3, 4, 5, 6]} | {"units": "°C"}, "PMSL.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "PS.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "TOT_PREC.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "mm"}, @@ -64,7 +64,7 @@ def _fallback(): "U_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "V_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "TD_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "°C"}, - "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0, "vmax": 6} | {"units": "°C"}, + "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 1, 2, 3, 4, 5, 6]} | {"units": "°C"}, "PMSL.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "PS.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "TOT_PREC.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "mm"}, @@ -75,7 +75,7 @@ def _fallback(): "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "°C"}, - "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "vmin": -5.5, "vmax": 5.5} | {"units": "°C"}, + "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : [-5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5]} | {"units": "°C"}, "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "PS.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "TOT_PREC.BIAS.spatial": {"cmap": plt.get_cmap("BrBG", 11)} | {"units": "mm"} diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index af4be700..4ee6e9da 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -100,41 +100,27 @@ def _(init_hour, lead_time, metric, param, season, verif_file, xr): return ds, var @app.cell -def _(CMAP_DEFAULTS, ekp, np): +def _(CMAP_DEFAULTS, ekp): def get_style(param, metric, units_override=None): """Get style and colormap settings for the plot. Needed because cmap/norm does not work in Style(colors=cmap), still needs to be passed as arguments to tripcolor()/tricontourf(). """ - metric_key = f"{param}.{metric}.spatial" cfg = CMAP_DEFAULTS[metric_key] if metric_key in CMAP_DEFAULTS else CMAP_DEFAULTS[param] units = units_override if units_override is not None else cfg.get("units", "") - - vmin = cfg.get("vmin", None) - vmax = cfg.get("vmax", None) - levels = cfg.get("levels", None) - - # For the case of Bias, construct symmetric levels: - if metric == "BIAS": - levels = np.linspace(start = vmin, stop = vmax, num=12) - elif metric == "RMSE": - levels = np.linspace(start = vmin, stop = vmax, num=7) - elif metric == "MAE": - levels = np.linspace(start = vmin, stop = vmax, num=7) - return { "style": ekp.styles.Style( - levels=cfg.get("bounds", levels), + levels=cfg.get("bounds", cfg.get("levels", None)), extend="both", units=units, colors=cfg.get("colors", None), ), "norm": cfg.get("norm", None), "cmap": cfg.get("cmap", None), - "levels": levels, - "vmin": vmin, - "vmax": vmax, + "levels": cfg.get("levels", None), + "vmin": cfg.get("vmin", None), + "vmax": cfg.get("vmax", None), "colors": cfg.get("colors", None), } return (get_style,) From 69e7b36dce3537ff12c2986e6b0f3d853f5a5967 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Tue, 3 Feb 2026 17:18:29 +0100 Subject: [PATCH 34/58] New Rule for Map plots of Baselines. --- workflow/rules/plot.smk | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 298784a2..ee113266 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -75,7 +75,7 @@ rule plot_summary_stat_maps: script="workflow/scripts/plot_summary_stat_maps.mo.py", verif_file=OUT_ROOT / "data/runs/{run_id}/verif_aggregated.nc", output: - OUT_ROOT / "results/{experiment}/metrics/spatial/{run_id}/{param}_{metric}_{leadtime}_{region}.png", + OUT_ROOT / "results/{experiment}/metrics/spatial/{run_id}/{param}_{metric}_{region}_{leadtime}.png", wildcard_constraints: leadtime=r"\d+", # only digits resources: @@ -98,4 +98,16 @@ rule plot_summary_stat_maps: # marimo edit {input.script} -- \ # --input {input.verif_file} --outfn {output[0]} --region {wildcards.region} \ # --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} \ - """ \ No newline at end of file + """ + +use rule plot_summary_stat_maps as plot_summary_stat_maps_baseline with: + input: + script="workflow/scripts/plot_summary_stat_maps.mo.py", + verif_file=OUT_ROOT / "data/baselines/{baseline_id}/verif_aggregated.nc", + output: + OUT_ROOT / "results/{experiment}/metrics/spatial/{baseline_id}/{param}_{metric}_{region}_{leadtime}.png", + params: + nc_out_dir=lambda wc: ( + Path(OUT_ROOT) / f"data/baselines/{wc.baseline_id}/verif_aggregated.nc" + # not sure if this is actually needed. Verification file is already specified above as input. Leave it for the time being. + ).resolve() \ No newline at end of file From 1353788e151fdcedafca79f786bdc29a30cca03a Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Wed, 4 Feb 2026 15:14:09 +0100 Subject: [PATCH 35/58] Different File naming (Region first -> Order) Required accompanying change in Snakefile. --- workflow/Snakefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/workflow/Snakefile b/workflow/Snakefile index 64a30694..d21469bb 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -53,7 +53,7 @@ rule experiment_all: run_id=collect_all_candidates(), ), expand( - OUT_ROOT / "results/{experiment}/metrics/spatial/{run_id}/{param}_{metric}_{leadtime}_{region}.png", + OUT_ROOT / "results/{experiment}/metrics/spatial/{run_id}/{param}_{metric}_{region}_{leadtime}.png", run_id=collect_all_candidates(), leadtime=[6, 12], metric=["BIAS", "RMSE", "MAE"], From 569759a1f46a6730cf8ac7843025d2bb182d5b48 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 5 Feb 2026 16:00:40 +0100 Subject: [PATCH 36/58] Temporarily remove rule for maps for baselines. --- workflow/rules/plot.smk | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index ee113266..54977ac8 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -100,14 +100,14 @@ rule plot_summary_stat_maps: # --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} \ """ -use rule plot_summary_stat_maps as plot_summary_stat_maps_baseline with: - input: - script="workflow/scripts/plot_summary_stat_maps.mo.py", - verif_file=OUT_ROOT / "data/baselines/{baseline_id}/verif_aggregated.nc", - output: - OUT_ROOT / "results/{experiment}/metrics/spatial/{baseline_id}/{param}_{metric}_{region}_{leadtime}.png", - params: - nc_out_dir=lambda wc: ( - Path(OUT_ROOT) / f"data/baselines/{wc.baseline_id}/verif_aggregated.nc" - # not sure if this is actually needed. Verification file is already specified above as input. Leave it for the time being. - ).resolve() \ No newline at end of file +# use rule plot_summary_stat_maps as plot_summary_stat_maps_baseline with: +# input: +# script="workflow/scripts/plot_summary_stat_maps.mo.py", +# verif_file=OUT_ROOT / "data/baselines/{baseline_id}/verif_aggregated.nc", +# output: +# OUT_ROOT / "results/{experiment}/metrics/spatial/{baseline_id}/{param}_{metric}_{region}_{leadtime}.png", +# params: +# nc_out_dir=lambda wc: ( +# Path(OUT_ROOT) / f"data/baselines/{wc.baseline_id}/verif_aggregated.nc" +# # not sure if this is actually needed. Verification file is already specified above as input. Leave it for the time being. +# ).resolve() \ No newline at end of file From bda0e0d51a2f1d8d0adb104a64d37493eed976a3 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 5 Feb 2026 16:03:04 +0100 Subject: [PATCH 37/58] Renamed domain "switzerland_small". This has caused problems before. --- src/plotting/__init__.py | 2 +- workflow/Snakefile | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/plotting/__init__.py b/src/plotting/__init__.py index 520e35ae..cb808df8 100644 --- a/src/plotting/__init__.py +++ b/src/plotting/__init__.py @@ -37,7 +37,7 @@ "extent": [0, 17.5, 40.5, 53.0], "projection": _PROJECTIONS["orthographic"], }, - "switzerland_small": { + "switzerlandsmall": { "extent": [5.5, 11.0, 45.5, 48.0], "projection": _PROJECTIONS["orthographic"], }, diff --git a/workflow/Snakefile b/workflow/Snakefile index d21469bb..3bcee554 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -56,9 +56,9 @@ rule experiment_all: OUT_ROOT / "results/{experiment}/metrics/spatial/{run_id}/{param}_{metric}_{region}_{leadtime}.png", run_id=collect_all_candidates(), leadtime=[6, 12], - metric=["BIAS", "RMSE", "MAE"], - param=["TD_2M", "T_2M", "PMSL", "PS", "TOT_PREC"], - region=["switzerland", "switzerland_small"], + metric=["BIAS", "RMSE"], + param=["T_2M", "TOT_PREC"], + region=["switzerland", "switzerlandsmall"], experiment=EXPERIMENT_HASH ) From ccf6279e50ac762a55055c35946fb56875fd6cb1 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 5 Feb 2026 16:04:24 +0100 Subject: [PATCH 38/58] Generate Colour breaks more elegantly. --- src/plotting/colormap_defaults.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index c9968240..dd6530a3 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -75,7 +75,7 @@ def _fallback(): "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "°C"}, - "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : [-5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5]} | {"units": "°C"}, + "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -5.5, stop = 6, step = 1)} | {"units": "°C"}, "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "PS.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "TOT_PREC.BIAS.spatial": {"cmap": plt.get_cmap("BrBG", 11)} | {"units": "mm"} From 258b2dea37c0376944cf92dcb86f48c1ccc359ab Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 5 Feb 2026 16:24:41 +0100 Subject: [PATCH 39/58] Colour breaks from function for all T2m metrics. --- src/plotting/colormap_defaults.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index dd6530a3..d7b35745 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -55,7 +55,7 @@ def _fallback(): "U_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "V_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "TD_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "°C"}, - "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 1, 2, 3, 4, 5, 6]} | {"units": "°C"}, + "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": np.arange(start = 0, stop = 6.1, step = 1)} | {"units": "°C"}, "PMSL.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "PS.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "TOT_PREC.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "mm"}, @@ -64,7 +64,7 @@ def _fallback(): "U_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "V_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "TD_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "°C"}, - "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 1, 2, 3, 4, 5, 6]} | {"units": "°C"}, + "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": np.arange(start = 0, stop = 6.1, step = 1)} | {"units": "°C"}, "PMSL.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "PS.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "TOT_PREC.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "mm"}, @@ -75,7 +75,7 @@ def _fallback(): "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "°C"}, - "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -5.5, stop = 6, step = 1)} | {"units": "°C"}, + "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -5.5, stop = 5.6, step = 1)} | {"units": "°C"}, "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "PS.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "TOT_PREC.BIAS.spatial": {"cmap": plt.get_cmap("BrBG", 11)} | {"units": "mm"} From 5c0151b8e78d9573ea1eabffc8054de285dfdbef Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Mon, 9 Feb 2026 19:09:45 +0100 Subject: [PATCH 40/58] Map plots now also for baselines. --- workflow/Snakefile | 11 ++++++++++- workflow/rules/plot.smk | 33 ++++++++++++++++++++++++++++++++- 2 files changed, 42 insertions(+), 2 deletions(-) diff --git a/workflow/Snakefile b/workflow/Snakefile index 3bcee554..991c288d 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -53,13 +53,22 @@ rule experiment_all: run_id=collect_all_candidates(), ), expand( - OUT_ROOT / "results/{experiment}/metrics/spatial/{run_id}/{param}_{metric}_{region}_{leadtime}.png", + OUT_ROOT / "results/{experiment}/metrics/spatial/runs/{run_id}/{param}_{metric}_{region}_{leadtime}.png", run_id=collect_all_candidates(), leadtime=[6, 12], metric=["BIAS", "RMSE"], param=["T_2M", "TOT_PREC"], region=["switzerland", "switzerlandsmall"], experiment=EXPERIMENT_HASH + ), + expand( + OUT_ROOT / "results/{experiment}/metrics/spatial/baselines/{baseline_id}/{param}_{metric}_{region}_{leadtime}.png", + baseline_id=collect_all_baselines(), + leadtime=[6, 12], + metric=["BIAS", "RMSE"], + param=["T_2M", "TOT_PREC"], + region=["switzerland", "switzerlandsmall"], + experiment=EXPERIMENT_HASH ) diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 54977ac8..0190761c 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -75,7 +75,7 @@ rule plot_summary_stat_maps: script="workflow/scripts/plot_summary_stat_maps.mo.py", verif_file=OUT_ROOT / "data/runs/{run_id}/verif_aggregated.nc", output: - OUT_ROOT / "results/{experiment}/metrics/spatial/{run_id}/{param}_{metric}_{region}_{leadtime}.png", + OUT_ROOT / "results/{experiment}/metrics/spatial/runs/{run_id}/{param}_{metric}_{region}_{leadtime}.png", wildcard_constraints: leadtime=r"\d+", # only digits resources: @@ -100,6 +100,37 @@ rule plot_summary_stat_maps: # --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} \ """ +rule plot_summary_stat_maps_baseline: + localrule: True + input: + script="workflow/scripts/plot_summary_stat_maps.mo.py", + verif_file=OUT_ROOT / "data/baselines/{baseline_id}/verif_aggregated.nc", + output: + OUT_ROOT / "results/{experiment}/metrics/spatial/baselines/{baseline_id}/{param}_{metric}_{region}_{leadtime}.png", + wildcard_constraints: + leadtime=r"\d+", # only digits + resources: + slurm_partition="postproc", + cpus_per_task=1, + runtime="10m", + params: + nc_out_dir=lambda wc: ( + Path(OUT_ROOT) / f"data/baselines/{wc.baseline_id}/verif_aggregated.nc" + # not sure how to do this, because the baselines are in, e.g., output/data/baselines/COSMO-E/verif_aggregated.nc + # and the runs are in output/data/runs/runID/verif_aggregated.nc + ).resolve(), + shell: + """ + export ECCODES_DEFINITION_PATH=$(realpath .venv/share/eccodes-cosmo-resources/definitions) + python {input.script} \ + --input {input.verif_file} --outfn {output[0]} --region {wildcards.region} \ + --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} \ + # interactive editing (needs to set localrule: True and use only one core) + # marimo edit {input.script} -- \ + # --input {input.verif_file} --outfn {output[0]} --region {wildcards.region} \ + # --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} \ + """ + # use rule plot_summary_stat_maps as plot_summary_stat_maps_baseline with: # input: # script="workflow/scripts/plot_summary_stat_maps.mo.py", From 7b50809e8cc9fdcd4e684cfed1dd90196ea0235a Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Mon, 9 Feb 2026 19:10:28 +0100 Subject: [PATCH 41/58] Fix memory leak occurring during creation of html. --- .../scripts/report_experiment_dashboard.py | 68 ++++++++++++++++++- 1 file changed, 65 insertions(+), 3 deletions(-) diff --git a/workflow/scripts/report_experiment_dashboard.py b/workflow/scripts/report_experiment_dashboard.py index 327547b0..fc6ed6a9 100644 --- a/workflow/scripts/report_experiment_dashboard.py +++ b/workflow/scripts/report_experiment_dashboard.py @@ -28,19 +28,81 @@ def program_summary_log(args): LOG.info("=" * 80) +def dataset_metrics_to_dataframe( + ds: xr.Dataset, + forbidden_dims=("values",), + metric_dims=("source", "season", "init_hour", "region", "lead_time", "eps"), +): + """ + Convert a verification xarray.Dataset to a pandas.DataFrame compatible with + the old `to_array("stack").to_dataframe()` workflow. + + Returns columns: + - metric_dims... + - stack (e.g. "T_2M.MAE") + - value + """ + + # keep only non-spatial metric variables + metric_vars = [] + for v in ds.data_vars: + da = ds[v] + + # skip spatial metrics + if "spatial" in v: + continue + + # skip forbidden dimensions + if any(d in da.dims for d in forbidden_dims): + continue + + # only keep vars whose dims are a subset of expected metric dims + if not set(da.dims).issubset(metric_dims): + continue + + metric_vars.append(v) + + if not metric_vars: + raise ValueError("No compatible metric variables found in dataset") + + # stack variables exactly like the original code + df = ( + ds[metric_vars] + .to_array("stack") + .to_dataframe(name="value") + .reset_index() + ) + + return df + def main(args): program_summary_log(args) # Load, de-duplicate lead_time, and keep best provider per source (same logic as verif_plot_metrics) - dfs = [xr.open_dataset(f) for f in args.verif_files] + drop_variables = ["TOT_PREC.MAE.spatial", "TOT_PREC.RMSE.spatial", "TOT_PREC.BIAS.spatial", + "PMSL.MAE.spatial", "PMSL.RMSE.spatial", "PMSL.BIAS.spatial", + "PS.MAE.spatial", "PS.RMSE.spatial", "PS.BIAS.spatial", + "V_10M.MAE.spatial", "V_10M.RMSE.spatial", "V_10M.BIAS.spatial", + "U_10M.MAE.spatial", "U_10M.RMSE.spatial", "U_10M.BIAS.spatial", + "TD_2M.MAE.spatial", "TD_2M.RMSE.spatial", "TD_2M.BIAS.spatial", + "T_2M.MAE.spatial", "T_2M.RMSE.spatial", "T_2M.BIAS.spatial"] + + dfs = [xr.open_dataset(f, drop_variables = drop_variables) for f in args.verif_files] dfs = [_ensure_unique_lead_time(d) for d in dfs] dfs = _select_best_sources(dfs) ds = xr.concat(dfs, dim="source", join="outer") LOG.info("Loaded verification netcdf: \n%s", ds) # extract only non-spatial variables to pd.DataFrame - nonspatial_vars = [d for d in ds.data_vars if "spatial" not in d] - df = ds[nonspatial_vars].to_array("stack").to_dataframe(name="value").reset_index() + # nonspatial_vars = [d for d in ds.data_vars if "spatial" not in d] + # df = ds[nonspatial_vars].to_array("stack").to_dataframe(name="value").reset_index() + + df = dataset_metrics_to_dataframe( + ds, + forbidden_dims=("values",), # critical! + metric_dims=("source", "season", "init_hour", "region", "lead_time", "eps"), + ) + df[["param", "metric"]] = df["stack"].str.split(".", n=1, expand=True) df.drop(columns=["stack"], inplace=True) df["lead_time"] = df["lead_time"].dt.total_seconds() / 3600 From c3d37a20fd3db2ff5815003e462ca623054dcac3 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Tue, 10 Feb 2026 17:01:08 +0100 Subject: [PATCH 42/58] "Alias" Rule for Plotting maps of Baselines This does the same thing with less code. --- workflow/rules/plot.smk | 37 +++---------------------------------- 1 file changed, 3 insertions(+), 34 deletions(-) diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 0190761c..2924201e 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -100,45 +100,14 @@ rule plot_summary_stat_maps: # --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} \ """ -rule plot_summary_stat_maps_baseline: - localrule: True +use rule plot_summary_stat_maps as plot_summary_stat_maps_baseline with: input: script="workflow/scripts/plot_summary_stat_maps.mo.py", verif_file=OUT_ROOT / "data/baselines/{baseline_id}/verif_aggregated.nc", output: OUT_ROOT / "results/{experiment}/metrics/spatial/baselines/{baseline_id}/{param}_{metric}_{region}_{leadtime}.png", - wildcard_constraints: - leadtime=r"\d+", # only digits - resources: - slurm_partition="postproc", - cpus_per_task=1, - runtime="10m", params: nc_out_dir=lambda wc: ( Path(OUT_ROOT) / f"data/baselines/{wc.baseline_id}/verif_aggregated.nc" - # not sure how to do this, because the baselines are in, e.g., output/data/baselines/COSMO-E/verif_aggregated.nc - # and the runs are in output/data/runs/runID/verif_aggregated.nc - ).resolve(), - shell: - """ - export ECCODES_DEFINITION_PATH=$(realpath .venv/share/eccodes-cosmo-resources/definitions) - python {input.script} \ - --input {input.verif_file} --outfn {output[0]} --region {wildcards.region} \ - --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} \ - # interactive editing (needs to set localrule: True and use only one core) - # marimo edit {input.script} -- \ - # --input {input.verif_file} --outfn {output[0]} --region {wildcards.region} \ - # --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} \ - """ - -# use rule plot_summary_stat_maps as plot_summary_stat_maps_baseline with: -# input: -# script="workflow/scripts/plot_summary_stat_maps.mo.py", -# verif_file=OUT_ROOT / "data/baselines/{baseline_id}/verif_aggregated.nc", -# output: -# OUT_ROOT / "results/{experiment}/metrics/spatial/{baseline_id}/{param}_{metric}_{region}_{leadtime}.png", -# params: -# nc_out_dir=lambda wc: ( -# Path(OUT_ROOT) / f"data/baselines/{wc.baseline_id}/verif_aggregated.nc" -# # not sure if this is actually needed. Verification file is already specified above as input. Leave it for the time being. -# ).resolve() \ No newline at end of file + # not sure if this is actually needed. Verification file is already specified above as input. Leave it for the time being. + ).resolve() \ No newline at end of file From 8323fbda7017e29b80d46cdba5831d1923332a15 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Wed, 11 Feb 2026 10:19:44 +0100 Subject: [PATCH 43/58] Reorganisation of Domains. --- src/plotting/__init__.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/plotting/__init__.py b/src/plotting/__init__.py index cb808df8..c905fcc3 100644 --- a/src/plotting/__init__.py +++ b/src/plotting/__init__.py @@ -29,15 +29,15 @@ "extent": [-16.0, 25.0, 30.0, 65.0], "projection": _PROJECTIONS["orthographic"], }, + # The domains which are originally called "centraleurope" and "switzerland" + # are mostly the same. I suggest making domain "switzerland" much smaller, + # so that more spatial detail can be seen, especially in the complex + # topography of the alps. "centraleurope": { "extent": [-2.6, 19.5, 40.2, 52.3], "projection": _PROJECTIONS["orthographic"], }, "switzerland": { - "extent": [0, 17.5, 40.5, 53.0], - "projection": _PROJECTIONS["orthographic"], - }, - "switzerlandsmall": { "extent": [5.5, 11.0, 45.5, 48.0], "projection": _PROJECTIONS["orthographic"], }, From d03781cf86b0f558097b688f0effbc5cb6fef363 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Wed, 11 Feb 2026 11:57:17 +0100 Subject: [PATCH 44/58] Introduce Seasonal Stratification. --- workflow/Snakefile | 22 ++++++++++--------- workflow/rules/plot.smk | 6 +++-- workflow/scripts/plot_summary_stat_maps.mo.py | 3 ++- 3 files changed, 18 insertions(+), 13 deletions(-) diff --git a/workflow/Snakefile b/workflow/Snakefile index 991c288d..51d13b41 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -53,21 +53,23 @@ rule experiment_all: run_id=collect_all_candidates(), ), expand( - OUT_ROOT / "results/{experiment}/metrics/spatial/runs/{run_id}/{param}_{metric}_{region}_{leadtime}.png", + OUT_ROOT / "results/{experiment}/metrics/spatial/runs/{run_id}/{param}_{metric}_{region}_{season}_{leadtime}.png", run_id=collect_all_candidates(), - leadtime=[6, 12], - metric=["BIAS", "RMSE"], - param=["T_2M", "TOT_PREC"], - region=["switzerland", "switzerlandsmall"], + leadtime=[6], + metric=["BIAS"], + param=["T_2M"], + region=["centraleurope", "switzerland"], + season=["all", "DJF", "MAM", "JJA", "SON"], experiment=EXPERIMENT_HASH ), expand( - OUT_ROOT / "results/{experiment}/metrics/spatial/baselines/{baseline_id}/{param}_{metric}_{region}_{leadtime}.png", + OUT_ROOT / "results/{experiment}/metrics/spatial/baselines/{baseline_id}/{param}_{metric}_{region}_{season}_{leadtime}.png", baseline_id=collect_all_baselines(), - leadtime=[6, 12], - metric=["BIAS", "RMSE"], - param=["T_2M", "TOT_PREC"], - region=["switzerland", "switzerlandsmall"], + leadtime=[6], + metric=["BIAS"], + param=["T_2M"], + region=["centraleurope", "switzerland"], + season=["all", "DJF", "MAM", "JJA", "SON"], experiment=EXPERIMENT_HASH ) diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 2924201e..10f85c47 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -75,7 +75,7 @@ rule plot_summary_stat_maps: script="workflow/scripts/plot_summary_stat_maps.mo.py", verif_file=OUT_ROOT / "data/runs/{run_id}/verif_aggregated.nc", output: - OUT_ROOT / "results/{experiment}/metrics/spatial/runs/{run_id}/{param}_{metric}_{region}_{leadtime}.png", + OUT_ROOT / "results/{experiment}/metrics/spatial/runs/{run_id}/{param}_{metric}_{region}_{season}_{leadtime}.png", wildcard_constraints: leadtime=r"\d+", # only digits resources: @@ -94,10 +94,12 @@ rule plot_summary_stat_maps: python {input.script} \ --input {input.verif_file} --outfn {output[0]} --region {wildcards.region} \ --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} \ + --season {wildcards.season} \ # interactive editing (needs to set localrule: True and use only one core) # marimo edit {input.script} -- \ # --input {input.verif_file} --outfn {output[0]} --region {wildcards.region} \ # --param {wildcards.param} --leadtime {wildcards.leadtime} --metric {wildcards.metric} \ + # --season {wildcards.season} \ """ use rule plot_summary_stat_maps as plot_summary_stat_maps_baseline with: @@ -105,7 +107,7 @@ use rule plot_summary_stat_maps as plot_summary_stat_maps_baseline with: script="workflow/scripts/plot_summary_stat_maps.mo.py", verif_file=OUT_ROOT / "data/baselines/{baseline_id}/verif_aggregated.nc", output: - OUT_ROOT / "results/{experiment}/metrics/spatial/baselines/{baseline_id}/{param}_{metric}_{region}_{leadtime}.png", + OUT_ROOT / "results/{experiment}/metrics/spatial/baselines/{baseline_id}/{param}_{metric}_{region}_{season}_{leadtime}.png", params: nc_out_dir=lambda wc: ( Path(OUT_ROOT) / f"data/baselines/{wc.baseline_id}/verif_aggregated.nc" diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index 4ee6e9da..4ff57749 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -138,6 +138,7 @@ def _( outfn, param, region, + season, var, ): # plot individual fields @@ -172,7 +173,7 @@ def _( # validtime = state["valid_time"].strftime("%Y%m%d%H%M") # # leadtime = int(state["lead_time"].total_seconds() // 3600) - fig.title(f"{metric} of {param}, Lead Time: {lead_time}") + fig.title(f"{metric} of {param}, Season: {season}, Lead Time: {lead_time}") fig.save(outfn, bbox_inches="tight", dpi=200) LOG.info(f"saved: {outfn}") From 48e67892cb1d690aa16b31951d50c08da24db5a9 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 12 Feb 2026 12:33:04 +0100 Subject: [PATCH 45/58] Definitive Colour Levels for 2m-Temperature --- src/plotting/colormap_defaults.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index d7b35745..12dd9420 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -55,7 +55,7 @@ def _fallback(): "U_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "V_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "TD_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "°C"}, - "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": np.arange(start = 0, stop = 6.1, step = 1)} | {"units": "°C"}, + "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": np.arange(start = 0, stop = 3.1, step = 0.5)} | {"units": "°C"}, "PMSL.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "PS.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "TOT_PREC.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "mm"}, @@ -64,7 +64,7 @@ def _fallback(): "U_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "V_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "TD_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "°C"}, - "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": np.arange(start = 0, stop = 6.1, step = 1)} | {"units": "°C"}, + "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": np.arange(start = 0, stop = 3.1, step = 0.5)} | {"units": "°C"}, "PMSL.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "PS.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "TOT_PREC.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "mm"}, @@ -75,7 +75,7 @@ def _fallback(): "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "°C"}, - "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -5.5, stop = 5.6, step = 1)} | {"units": "°C"}, + "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "PS.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "TOT_PREC.BIAS.spatial": {"cmap": plt.get_cmap("BrBG", 11)} | {"units": "mm"} From 1c3e62984566d039b7af8ae0ba47b81aa6c53184 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Fri, 13 Feb 2026 16:17:09 +0100 Subject: [PATCH 46/58] Colour levels for Precipitation complete. --- src/plotting/colormap_defaults.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index 12dd9420..bda0d86c 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -58,8 +58,10 @@ def _fallback(): "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": np.arange(start = 0, stop = 3.1, step = 0.5)} | {"units": "°C"}, "PMSL.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "PS.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, - "TOT_PREC.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "mm"}, - + "TOT_PREC.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.001, 0.0015, 0.002, 0.003, 0.004]} | {"units": "mm"}, + # would ideally want a 6th colour on the high end of the colour scale, but somehow + # matplotlib does not do that -> ? + # MAE: "U_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "V_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, @@ -67,7 +69,8 @@ def _fallback(): "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": np.arange(start = 0, stop = 3.1, step = 0.5)} | {"units": "°C"}, "PMSL.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "PS.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, - "TOT_PREC.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "mm"}, + "TOT_PREC.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.001, 0.0015, 0.002, 0.003, 0.004]} | {"units": "mm"}, + # the levels for precip are a bit on the bright side, but still worth keeping consistent with RMSE. # Bias: # diverging colour scheme for the Bias to reflect the nature of the data (can be positive or negative, symmetric). @@ -78,7 +81,7 @@ def _fallback(): "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "PS.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, - "TOT_PREC.BIAS.spatial": {"cmap": plt.get_cmap("BrBG", 11)} | {"units": "mm"} + "TOT_PREC.BIAS.spatial": {"cmap": plt.get_cmap("BrBG", 9), "levels": [-0.001, -0.0005, -0.00025, -0.0001, 0.0001, 0.00025, 0.0005, 0.001]} | {"units": "mm"} } CMAP_DEFAULTS = defaultdict(_fallback, _CMAP_DEFAULTS) From a4a5607fb4454b382dc3ffe58adbd9cb67a38b07 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Tue, 17 Feb 2026 12:14:42 +0100 Subject: [PATCH 47/58] Colour Levels for Dew-Point Temperature and similar style for MAE and RMSE of T2m. --- src/plotting/colormap_defaults.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index bda0d86c..612265b9 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -54,8 +54,8 @@ def _fallback(): # RMSE: "U_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "V_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, - "TD_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "°C"}, - "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": np.arange(start = 0, stop = 3.1, step = 0.5)} | {"units": "°C"}, + "TD_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, + "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "PMSL.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "PS.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "TOT_PREC.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.001, 0.0015, 0.002, 0.003, 0.004]} | {"units": "mm"}, @@ -65,8 +65,8 @@ def _fallback(): # MAE: "U_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "V_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, - "TD_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "°C"}, - "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": np.arange(start = 0, stop = 3.1, step = 0.5)} | {"units": "°C"}, + "TD_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, + "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "PMSL.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "PS.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "TOT_PREC.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.001, 0.0015, 0.002, 0.003, 0.004]} | {"units": "mm"}, @@ -77,7 +77,7 @@ def _fallback(): # Red-Blue colour scheme for all variables except precipitation, where a Brown-Green scheme is more suggestive. "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, - "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "°C"}, + "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "PS.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, From f245ea38785ad2901746dfd4baf68051d26c9484 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 19 Feb 2026 15:46:58 +0100 Subject: [PATCH 48/58] Calculate and Evaluate Wind Speed too --- src/plotting/colormap_defaults.py | 3 +++ workflow/Snakefile | 4 ++-- workflow/scripts/verif_single_init.py | 11 +++++++++-- 3 files changed, 14 insertions(+), 4 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index 612265b9..a9cd0b3d 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -54,6 +54,7 @@ def _fallback(): # RMSE: "U_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "V_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, + "WS_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "TD_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "PMSL.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, @@ -65,6 +66,7 @@ def _fallback(): # MAE: "U_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "V_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, + "WS_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "TD_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "PMSL.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, @@ -77,6 +79,7 @@ def _fallback(): # Red-Blue colour scheme for all variables except precipitation, where a Brown-Green scheme is more suggestive. "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, + "WS_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, diff --git a/workflow/Snakefile b/workflow/Snakefile index 51d13b41..a2989eb5 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -57,7 +57,7 @@ rule experiment_all: run_id=collect_all_candidates(), leadtime=[6], metric=["BIAS"], - param=["T_2M"], + param=["WS_10M"], region=["centraleurope", "switzerland"], season=["all", "DJF", "MAM", "JJA", "SON"], experiment=EXPERIMENT_HASH @@ -67,7 +67,7 @@ rule experiment_all: baseline_id=collect_all_baselines(), leadtime=[6], metric=["BIAS"], - param=["T_2M"], + param=["WS_10M"], region=["centraleurope", "switzerland"], season=["all", "DJF", "MAM", "JJA", "SON"], experiment=EXPERIMENT_HASH diff --git a/workflow/scripts/verif_single_init.py b/workflow/scripts/verif_single_init.py index 52d9e08f..44cb2ddb 100644 --- a/workflow/scripts/verif_single_init.py +++ b/workflow/scripts/verif_single_init.py @@ -111,8 +111,15 @@ def main(args: ScriptConfig): analysis, ) - # compute metrics and statistics + # before verifying, calculate wind speed: + for ds in [fcst, analysis]: + if "U_10M" in ds and "V_10M" in ds: + LOG.info("Calculating Wind Speed (WS_10M)...") + ds["WS_10M"] = (ds["U_10M"]**2 + ds["V_10M"]**2)**0.5 + # Optional: Add metadata for the netCDF output + ds["WS_10M"].attrs = {"units": "m/s", "long_name": "10m Wind Speed"} + # compute metrics and statistics results = verify(fcst, analysis, args.label, args.analysis_label, args.regions) LOG.info("Verification results:\n%s", results) @@ -150,7 +157,7 @@ def main(args: ScriptConfig): parser.add_argument( "--params", type=lambda x: x.split(","), - default=["T_2M", "TD_2M", "U_10M", "V_10M", "PS", "PMSL", "TOT_PREC"], + default=["T_2M", "TD_2M", "U_10M", "V_10M", "WS_10M", "PS", "PMSL", "TOT_PREC"], ) parser.add_argument( "--steps", From ccea8dbca6de9b20c978d8ed6f75fbe0b625a1e2 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 19 Feb 2026 16:55:34 +0100 Subject: [PATCH 49/58] Fix for Wind speed calculation Now works! --- workflow/scripts/verif_single_init.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/workflow/scripts/verif_single_init.py b/workflow/scripts/verif_single_init.py index 44cb2ddb..27b9cbf6 100644 --- a/workflow/scripts/verif_single_init.py +++ b/workflow/scripts/verif_single_init.py @@ -157,7 +157,7 @@ def main(args: ScriptConfig): parser.add_argument( "--params", type=lambda x: x.split(","), - default=["T_2M", "TD_2M", "U_10M", "V_10M", "WS_10M", "PS", "PMSL", "TOT_PREC"], + default=["T_2M", "TD_2M", "U_10M", "V_10M", "PS", "PMSL", "TOT_PREC"], ) parser.add_argument( "--steps", From ce4242fd5cc14b1b8e0426b4cd773fc240ff675c Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Mon, 23 Feb 2026 14:39:12 +0100 Subject: [PATCH 50/58] Name Wind Speed consistent with the previous def. --- src/plotting/colormap_defaults.py | 6 +++--- workflow/Snakefile | 4 ++-- workflow/scripts/verif_single_init.py | 6 +++--- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index a9cd0b3d..06790652 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -54,7 +54,7 @@ def _fallback(): # RMSE: "U_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "V_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, - "WS_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, + "SP_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "TD_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "PMSL.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, @@ -66,7 +66,7 @@ def _fallback(): # MAE: "U_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "V_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, - "WS_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, + "SP_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "TD_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "PMSL.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, @@ -79,7 +79,7 @@ def _fallback(): # Red-Blue colour scheme for all variables except precipitation, where a Brown-Green scheme is more suggestive. "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, - "WS_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, + "SP_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, diff --git a/workflow/Snakefile b/workflow/Snakefile index a2989eb5..aa1ab24a 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -57,7 +57,7 @@ rule experiment_all: run_id=collect_all_candidates(), leadtime=[6], metric=["BIAS"], - param=["WS_10M"], + param=["SP_10M"], region=["centraleurope", "switzerland"], season=["all", "DJF", "MAM", "JJA", "SON"], experiment=EXPERIMENT_HASH @@ -67,7 +67,7 @@ rule experiment_all: baseline_id=collect_all_baselines(), leadtime=[6], metric=["BIAS"], - param=["WS_10M"], + param=["SP_10M"], region=["centraleurope", "switzerland"], season=["all", "DJF", "MAM", "JJA", "SON"], experiment=EXPERIMENT_HASH diff --git a/workflow/scripts/verif_single_init.py b/workflow/scripts/verif_single_init.py index 27b9cbf6..dae6020a 100644 --- a/workflow/scripts/verif_single_init.py +++ b/workflow/scripts/verif_single_init.py @@ -114,10 +114,10 @@ def main(args: ScriptConfig): # before verifying, calculate wind speed: for ds in [fcst, analysis]: if "U_10M" in ds and "V_10M" in ds: - LOG.info("Calculating Wind Speed (WS_10M)...") - ds["WS_10M"] = (ds["U_10M"]**2 + ds["V_10M"]**2)**0.5 + LOG.info("Calculating Wind Speed (SP_10M)...") + ds["SP_10M"] = (ds["U_10M"]**2 + ds["V_10M"]**2)**0.5 # Optional: Add metadata for the netCDF output - ds["WS_10M"].attrs = {"units": "m/s", "long_name": "10m Wind Speed"} + ds["SP_10M"].attrs = {"units": "m/s", "long_name": "10m Wind Speed"} # compute metrics and statistics results = verify(fcst, analysis, args.label, args.analysis_label, args.regions) From 9628fae72db3f287faa9ce3b5bb06ee37fa0a4f7 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Tue, 24 Feb 2026 10:09:27 +0100 Subject: [PATCH 51/58] Verification files not temporary any more Some SLURM optimizations. --- workflow/rules/verif.smk | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/workflow/rules/verif.smk b/workflow/rules/verif.smk index 0566ccf2..4d60045d 100644 --- a/workflow/rules/verif.smk +++ b/workflow/rules/verif.smk @@ -27,13 +27,14 @@ rule verif_metrics_baseline: analysis_label=config["analysis"].get("label"), regions=REGION_TXT, output: - temp(OUT_ROOT / "data/baselines/{baseline_id}/{init_time}/verif.nc"), + OUT_ROOT / "data/baselines/{baseline_id}/{init_time}/verif.nc", log: OUT_ROOT / "logs/verif_metrics_baseline/{baseline_id}-{init_time}.log", resources: cpus_per_task=24, mem_mb=50_000, runtime="60m", + slurm_extra="--exclude=nid001229,nid001225,nid001226,nid001227,nid001230" shell: """ uv run {input.script} \ @@ -63,7 +64,7 @@ rule verif_metrics: inference_okfile=rules.execute_inference.output.okfile, analysis_zarr=config["analysis"].get("analysis_zarr"), output: - temp(OUT_ROOT / "data/runs/{run_id}/{init_time}/verif.nc"), + OUT_ROOT / "data/runs/{run_id}/{init_time}/verif.nc", # wildcard_constraints: # run_id="^" # to avoid ambiguitiy with run_baseline_verif # TODO: implement logic to use experiment name instead of run_id as wildcard @@ -81,6 +82,7 @@ rule verif_metrics: cpus_per_task=24, mem_mb=50_000, runtime="60m", + slurm_extra="--exclude=nid001229,nid001225,nid001226,nid001227,nid001230" shell: """ uv run {input.script} \ @@ -117,7 +119,8 @@ rule verif_metrics_aggregation: resources: cpus_per_task=24, mem_mb=250_000, - runtime="2h", + runtime="24h", + slurm_extra="--exclude=nid001229,nid001225,nid001226,nid001227,nid001230" shell: """ uv run {input.script} {input.verif_nc} --output {output} > {log} 2>&1 From 59d4b12a1369d33ed764550e242ad6f78d7f2464 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 26 Feb 2026 13:53:10 +0100 Subject: [PATCH 52/58] Precipitation plotting: meters to millimeters. --- src/plotting/colormap_defaults.py | 6 +- workflow/scripts/plot_summary_stat_maps.mo.py | 78 ++++++++++++++++++- 2 files changed, 80 insertions(+), 4 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index 06790652..4fc652e4 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -59,7 +59,7 @@ def _fallback(): "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "PMSL.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "PS.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, - "TOT_PREC.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.001, 0.0015, 0.002, 0.003, 0.004]} | {"units": "mm"}, + "TOT_PREC.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 1, 1.5, 2, 3, 4]} | {"units": "mm"}, # would ideally want a 6th colour on the high end of the colour scale, but somehow # matplotlib does not do that -> ? @@ -71,7 +71,7 @@ def _fallback(): "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "PMSL.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, "PS.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, - "TOT_PREC.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.001, 0.0015, 0.002, 0.003, 0.004]} | {"units": "mm"}, + "TOT_PREC.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 1, 1.5, 2, 3, 4]} | {"units": "mm"}, # the levels for precip are a bit on the bright side, but still worth keeping consistent with RMSE. # Bias: @@ -84,7 +84,7 @@ def _fallback(): "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "PS.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, - "TOT_PREC.BIAS.spatial": {"cmap": plt.get_cmap("BrBG", 9), "levels": [-0.001, -0.0005, -0.00025, -0.0001, 0.0001, 0.00025, 0.0005, 0.001]} | {"units": "mm"} + "TOT_PREC.BIAS.spatial": {"cmap": plt.get_cmap("BrBG", 9), "levels": [-1, -0.5, -0.25, -0.1, 0.1, 0.25, 0.5, 1]} | {"units": "mm"} } CMAP_DEFAULTS = defaultdict(_fallback, _CMAP_DEFAULTS) diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index 4ff57749..0db54b3f 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -125,6 +125,76 @@ def get_style(param, metric, units_override=None): } return (get_style,) +# @app.cell +# def _(LOG, np): +# """Preprocess fields with pint-based unit conversion and derived quantities.""" +# try: +# import pint # type: ignore + +# _ureg = pint.UnitRegistry() + +# def _k_to_c(arr): +# # robust conversion with pint, fallback if dtype unsupported +# try: +# return (_ureg.Quantity(arr, _ureg.kelvin).to(_ureg.degC)).magnitude +# except Exception: +# return arr - 273.15 + +# def _ms_to_knots(arr): +# # robust conversion with pint, fallback if dtype unsupported +# try: +# return ( +# _ureg.Quantity(arr, _ureg.meter / _ureg.second).to(_ureg.knot) +# ).magnitude +# except Exception: +# return arr * 1.943844 + +# def _m_to_mm(arr): +# # robust conversion with pint, fallback if dtype unsupported +# try: +# return (_ureg.Quantity(arr, _ureg.meter).to(_ureg.millimeter)).magnitude +# except Exception: +# return arr * 1000 + +# except Exception: +# LOG.warning("pint not available; falling back hardcoded conversions") + +# def _k_to_c(arr): +# return arr - 273.15 + +# def _ms_to_knots(arr): +# return arr * 1.943844 + +# def _m_to_mm(arr): +# return arr * 1000 + +# 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"] +# # temperature variables +# if param in ("T_2M", "TD_2M", "T", "TD"): +# return _k_to_c(fields[param]), "°C" +# # 10m wind speed (allow legacy 'uv' alias) +# if param == "SP_10M": +# u = fields["U_10M"] +# v = fields["V_10M"] +# return np.sqrt(u**2 + v**2), "m/s" +# # wind speed from standard-level components +# if param == "SP": +# u = fields["U"] +# v = fields["V"] +# return np.sqrt(u**2 + v**2), "m/s" +# if param == "TOT_PREC": +# return _m_to_mm(fields[param]), "mm" +# # default: passthrough +# return fields[param], None + +# return (preprocess_field,) @app.cell def _( @@ -162,7 +232,13 @@ def _( # # preprocess field (unit conversion, derived quantities) # field, units_override = preprocess_field(param, state) - plotter.plot_field(subplot, ds.values.ravel(), **get_style(var, metric)) + # Quick fix for precipitation (might have to use preprocess_field in the end (see above)) + if param == "TOT_PREC": + plot_vals = ds.values.ravel()*1000 + else: + plot_vals = ds.values.ravel() + + plotter.plot_field(subplot, plot_vals, **get_style(var, metric)) # subplot.ax.add_geometries( # state["lam_envelope"], # edgecolor="black", From 28ad4c5c343014109b4bd546e1d6ceeb04e59de4 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 26 Feb 2026 13:54:31 +0100 Subject: [PATCH 53/58] Map plotting not on default busy nodes. Change back (clean up in the end) --- workflow/rules/plot.smk | 1 + 1 file changed, 1 insertion(+) diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 10f85c47..23b049d7 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -82,6 +82,7 @@ rule plot_summary_stat_maps: slurm_partition="postproc", cpus_per_task=1, runtime="10m", + slurm_extra="--exclude=nid001229,nid001225,nid001226,nid001227,nid001230" params: nc_out_dir=lambda wc: ( Path(OUT_ROOT) / f"data/runs/{wc.run_id}/verif_aggregated.nc" From db60db5feded23a055c7257806c7fde374c1ce4b Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 26 Feb 2026 15:19:20 +0100 Subject: [PATCH 54/58] Preliminary Colour Levels for Wind Bias. Might have to slightly adjust for ICON (emulator). --- src/plotting/colormap_defaults.py | 10 +++++----- workflow/scripts/plot_summary_stat_maps.mo.py | 2 ++ 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index 4fc652e4..b46fde95 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -77,11 +77,11 @@ def _fallback(): # Bias: # diverging colour scheme for the Bias to reflect the nature of the data (can be positive or negative, symmetric). # Red-Blue colour scheme for all variables except precipitation, where a Brown-Green scheme is more suggestive. - "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, - "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, - "SP_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "m/s"}, - "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, - "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels" : np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, + "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -2.25, stop = 2.26, step = 0.5)} | {"units": "m/s"}, + "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -2.25, stop = 2.26, step = 0.5)} | {"units": "m/s"}, + "SP_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -2.25, stop = 2.26, step = 0.5)} | {"units": "m/s"}, + "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, + "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "PS.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, "TOT_PREC.BIAS.spatial": {"cmap": plt.get_cmap("BrBG", 9), "levels": [-1, -0.5, -0.25, -0.1, 0.1, 0.25, 0.5, 1]} | {"units": "mm"} diff --git a/workflow/scripts/plot_summary_stat_maps.mo.py b/workflow/scripts/plot_summary_stat_maps.mo.py index 0db54b3f..728f49db 100644 --- a/workflow/scripts/plot_summary_stat_maps.mo.py +++ b/workflow/scripts/plot_summary_stat_maps.mo.py @@ -233,6 +233,8 @@ def _( # field, units_override = preprocess_field(param, state) # Quick fix for precipitation (might have to use preprocess_field in the end (see above)) + # for wind speed, preprocess_field only has conversion from m/s to knots + # (not the other way around), so I assume the values are in m/s if param == "TOT_PREC": plot_vals = ds.values.ravel()*1000 else: From 2fcfd8a339050a69f02d93c5dc44cdc729064bbb Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Thu, 26 Feb 2026 15:56:52 +0100 Subject: [PATCH 55/58] Replace EXPERIMENT_HASH by EXPERIMENT_NAME Solves a Problem that occurred due to the merge. --- workflow/Snakefile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/workflow/Snakefile b/workflow/Snakefile index 695abc77..2c11ecb9 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -134,7 +134,7 @@ rule experiment_all: param=["SP_10M"], region=["centraleurope", "switzerland"], season=["all", "DJF", "MAM", "JJA", "SON"], - experiment=EXPERIMENT_HASH + experiment=EXPERIMENT_NAME ), expand( OUT_ROOT / "results/{experiment}/metrics/spatial/baselines/{baseline_id}/{param}_{metric}_{region}_{season}_{leadtime}.png", @@ -144,7 +144,7 @@ rule experiment_all: param=["SP_10M"], region=["centraleurope", "switzerland"], season=["all", "DJF", "MAM", "JJA", "SON"], - experiment=EXPERIMENT_HASH + experiment=EXPERIMENT_NAME ) From c00c5a074444fa29201090ddeb615c0b51b16491 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Mon, 2 Mar 2026 11:34:47 +0100 Subject: [PATCH 56/58] Colour Levels for Pressure Bias. --- src/plotting/colormap_defaults.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index b46fde95..21bde983 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -82,8 +82,8 @@ def _fallback(): "SP_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -2.25, stop = 2.26, step = 0.5)} | {"units": "m/s"}, "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, - "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, - "PS.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11)} | {"units": "Pa"}, + "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -110, stop = 111, step = 20)} | {"units": "Pa"}, + "PS.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -110, stop = 111, step = 20)} | {"units": "Pa"}, "TOT_PREC.BIAS.spatial": {"cmap": plt.get_cmap("BrBG", 9), "levels": [-1, -0.5, -0.25, -0.1, 0.1, 0.25, 0.5, 1]} | {"units": "mm"} } From 23d8c9c55a8bf4d68dc0057b25050bbc0a11f399 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Mon, 9 Mar 2026 16:56:31 +0100 Subject: [PATCH 57/58] Colour Levels for Pressure MAE and RMSE --- src/plotting/colormap_defaults.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index 21bde983..a555d1aa 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -57,8 +57,8 @@ def _fallback(): "SP_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "TD_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, - "PMSL.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, - "PS.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, + "PMSL.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 7), "levels": [0, 50, 100, 150, 200, 250, 300, 350]} | {"units": "Pa"}, + "PS.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 7), "levels": [0, 50, 100, 150, 200, 250, 300, 350]} | {"units": "Pa"}, "TOT_PREC.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 1, 1.5, 2, 3, 4]} | {"units": "mm"}, # would ideally want a 6th colour on the high end of the colour scale, but somehow # matplotlib does not do that -> ? @@ -69,17 +69,17 @@ def _fallback(): "SP_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, "TD_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, - "PMSL.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, - "PS.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "Pa"}, + "PMSL.MAE.spatial": {"cmap": plt.get_cmap("Reds", 7), "levels": [0, 50, 100, 150, 200, 250, 300, 350]} | {"units": "Pa"}, + "PS.MAE.spatial": {"cmap": plt.get_cmap("Reds", 7), "levels": [0, 50, 100, 150, 200, 250, 300, 350]} | {"units": "Pa"}, "TOT_PREC.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 1, 1.5, 2, 3, 4]} | {"units": "mm"}, # the levels for precip are a bit on the bright side, but still worth keeping consistent with RMSE. # Bias: # diverging colour scheme for the Bias to reflect the nature of the data (can be positive or negative, symmetric). # Red-Blue colour scheme for all variables except precipitation, where a Brown-Green scheme is more suggestive. - "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -2.25, stop = 2.26, step = 0.5)} | {"units": "m/s"}, - "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -2.25, stop = 2.26, step = 0.5)} | {"units": "m/s"}, - "SP_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -2.25, stop = 2.26, step = 0.5)} | {"units": "m/s"}, + "U_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 9), "levels": np.arange(start = -2.25, stop = 2.26, step = 0.5)} | {"units": "m/s"}, + "V_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 9), "levels": np.arange(start = -2.25, stop = 2.26, step = 0.5)} | {"units": "m/s"}, + "SP_10M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 9), "levels": np.arange(start = -2.25, stop = 2.26, step = 0.5)} | {"units": "m/s"}, "TD_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, "T_2M.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -2.75, stop = 2.76, step = 0.5)} | {"units": "°C"}, "PMSL.BIAS.spatial": {"cmap": plt.get_cmap("RdBu_r", 11), "levels": np.arange(start = -110, stop = 111, step = 20)} | {"units": "Pa"}, From 84edd34a6c4a994ecd123ca07881655fc53f8b49 Mon Sep 17 00:00:00 2001 From: Louis-Frey Date: Tue, 10 Mar 2026 17:18:26 +0100 Subject: [PATCH 58/58] Colour Levels for 10m Wind MAE and RMSE --- src/plotting/colormap_defaults.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/plotting/colormap_defaults.py b/src/plotting/colormap_defaults.py index a555d1aa..7ad4878a 100644 --- a/src/plotting/colormap_defaults.py +++ b/src/plotting/colormap_defaults.py @@ -52,9 +52,9 @@ def _fallback(): # always start at 0 so that the saturation of the colour corresponds to the error magnitude. # RMSE: - "U_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, - "V_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, - "SP_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, + "U_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "m/s"}, + "V_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "m/s"}, + "SP_10M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "m/s"}, "TD_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "T_2M.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "PMSL.RMSE.spatial": {"cmap": plt.get_cmap("Reds", 7), "levels": [0, 50, 100, 150, 200, 250, 300, 350]} | {"units": "Pa"}, @@ -64,9 +64,9 @@ def _fallback(): # matplotlib does not do that -> ? # MAE: - "U_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, - "V_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, - "SP_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "vmin": 0} | {"units": "m/s"}, + "U_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "m/s"}, + "V_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "m/s"}, + "SP_10M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "m/s"}, "TD_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "T_2M.MAE.spatial": {"cmap": plt.get_cmap("Reds", 6), "levels": [0, 0.5, 1, 1.5, 2, 2.5, 3]} | {"units": "°C"}, "PMSL.MAE.spatial": {"cmap": plt.get_cmap("Reds", 7), "levels": [0, 50, 100, 150, 200, 250, 300, 350]} | {"units": "Pa"},