Add verification support for surface variable added with post-processors

config/varda-single-1.0.yaml

@@ -5,9 +5,12 @@ description: |
 config_label: varda-single-1.0
 
 dates:
-  start: 2025-03-01T00:00
-  end: 2025-03-03T00:00
-  frequency: 24h
+  # start: 2025-01-01T06:00
+  # end: 2025-12-31T06:00
+  # frequency: 294h
+  - 2025-01-01T06:00
+  - 2025-01-02T12:00
+  - 2025-01-03T18:00
 
 runs:
   - temporal_downscaler:
@@ -17,14 +20,17 @@ runs:
       config: resources/inference/configs/sgm-temporal-downscaler-global_trimedge_multi.yaml
       extra_requirements:
         - anemoi-datasets==0.5.35
-        # NOTE this pins the merge commit that introduced 'copy-prognostic-from-forecaster'
-        # TODO pin a release tag once the plugins are versioned. Do not pin the main branch to avoid stale code (changes do to main to not change the evalml hash)
-        - git+https://github.com/MeteoSwiss/anemoi-plugins-meteoswiss.git@b4835d76346e1e5ec0181884520da00e317fbfbc
-        # - anemoi-inference==0.11.0
+        - git+https://github.com/MeteoSwiss/anemoi-plugins-meteoswiss@1dc3eae81182e9c49ff07cf9d02cab8907f7e6dd
+        - git+https://github.com/ecmwf/anemoi-inference@88c1ec13632c570f246cc2580ce324e45f36a03b
       forecaster:
         checkpoint: https://service.meteoswiss.ch/mlstore#/models/varda-forecaster-det-sgm/versions/1
         config: resources/inference/configs/sgm-multidataset-forecaster-global-ich1-oper.yaml
         steps: 0/120/6
+        extra_requirements:
+          - anemoi-datasets==0.5.35
+          - git+https://github.com/MeteoSwiss/anemoi-plugins-meteoswiss@1dc3eae81182e9c49ff07cf9d02cab8907f7e6dd
+          - git+https://github.com/ecmwf/anemoi-inference@88c1ec13632c570f246cc2580ce324e45f36a03b
+
   - baseline:
       label: INCA
       root: /store_new/mch/msclim/INCA
@@ -41,20 +47,19 @@ runs:
 truth:
   label: SwissMetNet
   root: jretrievedwh:1,2
-  # To verify against SwissMetNet observations from the DWH via jretrievedwh,
-  # set instead (requires jretrievedwh.py on $PATH and $OPR_HOME set):
-  # Other selectors: root: jretrievedwh:locations=ARO,KLO,LUG
-  #                  root: jretrievedwh:bbox=45.8,47.8,5.9,10.5
-  #                  append ;stage=devt to target a non-prod DWH stage
 
 experiment:
   params:
     - T_2M
     - TD_2M
+    - U_10M
+    - V_10M
     - SP_10M
+    - DD_10M
+    - RELHUM_2M
     - TOT_PREC1
-    - TOT_PREC6
     - PMSL
+    - PS
   stratification:
     regions:
       - icon
@@ -66,11 +71,7 @@ experiment:
     root: /store_new/mch/msopr/ml/regions/Prognoseregionen_LV95_20220517
   thresholds:
     TOT_PREC1:
-      gt: [0.0, 0.1, 1, 5, 10]
-    TOT_PREC6:
-      gt: [0.0, 0.1, 5, 10, 20, 50]
-    SP_10M:
-      gt: [2.5, 5.0, 10.0]
+      gt: [0.0, 1, 5]
     T_2M:
       lt: [273.15]
       gt: [288.15, 298.15]
@@ -80,7 +81,7 @@ experiment:
       # - init_hour
       - season
   scorecards:
-    enabled: true
+    enabled: false
     sections:
       nowcasting:
         baseline: INCA
@@ -90,7 +91,7 @@ experiment:
           - "U_10M:RMSE,R2,ETS"
           - "V_10M:RMSE,R2,ETS"
           - "T_2M:RMSE,R2,ETS"
-          - "TOT_PREC1:RMSE,R2,ETS"
+          - "TOT_PREC:RMSE,R2,ETS"
       short_range:
         baseline: ICON-CH1-CTRL
         lead_times: "6/33/6"
@@ -116,6 +117,7 @@ showcase:
   params:
     - T_2M
     - SP_10M
+    - DD_10M
     - TOT_PREC1
   meteograms:
     enabled: false

resources/inference/configs/sgm-temporal-downscaler-global_trimedge_multi.yaml

@@ -57,6 +57,16 @@ constant_forcings:
 
 patch_metadata: resources/sgm-temporal-downscaler-ich1-oper-patch.yaml
 
+# The checkpoint never predicts 10-metre wind speed/direction or 2-metre
+# relative humidity, so anemoi-inference has no metadata of its own for them.
+typed_variables:
+  SP_10M:
+    mars: {param: SP_10M, levtype: sfc}
+  DD_10M:
+    mars: {param: DD_10M, levtype: sfc}
+  RELHUM_2M:
+    mars: {param: RELHUM_2M, levtype: sfc}
+
 post_processors:
   - accumulate_from_start_of_forecast: # accumulate tp from start of forecast
       accumulations:
@@ -79,15 +89,19 @@ output:
           - extract_mask: # removes global points
               mask: "lam_0/cutout_mask"
               as_slice: true
-          # here, the trimedge mask can be specified when available
-          # at overlap steps (multiples of the forecaster stride) replace the
-          # re-predicted prognostics with the forecaster's boundary values
           - forward_transform_filter:
               copy-prognostic-from-forecaster:
                 forecaster_path: forecaster/20*
                 common_leadtime: 6h
                 params_to_keep: [tp] # NOTE, if more diagnostics are added, update accordingly.
                 namer: *namer
+          - forward_transform_filter:
+              surface-diagnostics:
+                u_component: U_10M
+                v_component: V_10M
+                temperature: T_2M
+                dewpoint: TD_2M
+                variables: [SP_10M, DD_10M, RELHUM_2M]
 
     - grib:
         path: grib/ifs-{dateTime}_{step:03}.grib
@@ -108,6 +122,13 @@ output:
                 common_leadtime: 6h
                 params_to_keep: [tp] # NOTE, if more diagnostics are added, update accordingly.
                 namer: *namer
+          - forward_transform_filter:
+              surface-diagnostics:
+                u_component: U_10M
+                v_component: V_10M
+                temperature: T_2M
+                dewpoint: TD_2M
+                variables: [SP_10M, DD_10M, RELHUM_2M]
         modifiers:
           - patches:
             - variable:
@@ -132,6 +153,9 @@ output:
                 TOT_PREC: {"param": 228, "shortName": "tp"}
                 tp: {"param": 228, "shortName": "tp"}
                 z: {"param": 129, "shortName": "z"}
+                SP_10M: {"param": 207, "shortName": "10si"}
+                DD_10M: {"param": 260260, "shortName": "10wdir"}
+                RELHUM_2M: {"param": 260242, "shortName": "2r"}
                 "^q_(\\d+)$": {"param": 133, "shortName": "q"}
                 "^t_(\\d+)$": {"param": 130, "shortName": "t"}
                 "^u_(\\d+)$": {"param": 131, "shortName": "u"}

resources/inference/metadata/sgm-temporal-downscaler-ich1-oper-patch.yaml

@@ -49,6 +49,35 @@ dataset:
           param: T_2M
           step: 12
           time: 0
+      DD_10M:
+        # See the SP_10M comment below
+          date: 20050101
+          levtype: sfc
+          param: DD_10M
+          step: 12
+          time: 0
+      RELHUM_2M:
+        # See the SP_10M comment below
+        mars:
+          date: 20050101
+          levtype: sfc
+          param: RELHUM_2M
+          step: 12
+          time: 0
+      SP_10M:
+        # Not predicted by the checkpoint -- created mid-pipeline by the
+        # surface-diagnostics post-processor filter. Every
+        # field a forward_transform_filter chain touches (including ones a
+        # prior filter just created) must be a known variable here: anemoi-inference's
+        # TransformFilter.process() calls wrap_state(state, metadata.typed_variables)
+        # before EVERY filter in the chain, and that lookup is a hard KeyError
+        # for any field name not in this variables_metadata map
+        mars:
+          date: 20050101
+          levtype: sfc
+          param: SP_10M
+          step: 12
+          time: 0
       cos_julian_day:
         computed_forcing: true
         constant_in_time: false

resources/inference/templates/templates_index_icon.yaml

@@ -3,7 +3,7 @@
 - - {levtype: pl, param: [T, QV, U, V, OMEGA, FI]}
   - resources/icon-ch1-typeOfLevel=isobaricInhPa.grib
 
-- - {levtype: sfc, param: [T_2M, TD_2M, U_10M, V_10M]}
+- - {levtype: sfc, param: [T_2M, TD_2M, U_10M, V_10M, SP_10M, DD_10M, RELHUM_2M]}
   - resources/icon-ch1-typeOfLevel=heightAboveGround.grib
 
 - - {levtype: sfc, param: [FR_LAND, PS, FIS, T_G, SKT, SSO_SIGMA,SSO_STDH]}

src/data_input/__init__.py

@@ -7,6 +7,7 @@
 from typing import Callable, Literal, Any
 
 import earthkit.data as ekd
+import earthkit.meteo.thermo as ekdt
 import earthkit.meteo.vertical as ekdv
 import numpy as np
 import pandas as pd
@@ -54,7 +55,8 @@
 # Spatial/multi-field derivations: output param → required input params.
 DERIVED_PARAMS: dict[str, tuple[str, ...]] = {
     "SP_10M": ("U_10M", "V_10M"),
-    "SP": ("U", "V"),
+    "DD_10M": ("U_10M", "V_10M"),
+    "RELHUM_2M": ("T_2M", "TD_2M"),
 }
 
 
@@ -88,12 +90,26 @@ def compute_derived(ds: xr.Dataset, param: str) -> xr.DataArray:
             "name": "10m wind speed",
         }
         return da
-    if param == "SP":
-        da = (ds["U"] ** 2 + ds["V"] ** 2) ** 0.5
+    if param == "DD_10M":
+        # Meteorological convention: direction the wind is blowing FROM,
+        # clockwise from North. Matches the U/V <-> DD/FF convention used in
+        # load_obs_data_from_jretrieve (U = -FF*sin(DD), V = -FF*cos(DD)).
+        da = np.mod(np.degrees(np.arctan2(-ds["U_10M"], -ds["V_10M"])), 360.0)
         da.attrs["parameter"] = {
-            "shortName": "SP",
-            "units": "m/s",
-            "name": "Wind speed",
+            "shortName": "DD_10M",
+            "units": "degrees",
+            "name": "10m wind direction",
+        }
+        return da
+    if param == "RELHUM_2M":
+        # Same formula as the anemoi-inference surface-diagnostics post-processor
+        # (anemoi_plugins_meteoswiss.transform.filters.surface_diagnostics), so
+        # baselines are verified against the ML forecaster on a like-for-like basis.
+        da = ekdt.relative_humidity_from_dewpoint(ds["T_2M"], ds["TD_2M"])
+        da.attrs["parameter"] = {
+            "shortName": "RELHUM_2M",
+            "units": "%",
+            "name": "2m relative humidity",
         }
         return da
     raise ValueError(f"No recipe for derived variable '{param}'")
@@ -713,6 +729,7 @@ def load_obs_data_from_jretrieve(
         "SP_10M": "fkl010z0",
         "DD_10M": "dkl010z0",
         "VMAX_10M": "fkl010z1",
+        "RELHUM_2M": "ure200h0",
     }
     DWH_WIND_SPEED = "fkl010z0"
     DWH_WIND_DIR = "dkl010z0"
@@ -786,7 +803,7 @@ def load_truth_data(
     """Load truth data from an analysis Zarr dataset or DWH observations via jretrieve.
 
     Handles derived and aggregated params transparently (same contract as
-    load_forecast_data): SP_10M is computed from U_10M/V_10M; TOT_PREC6 is
+    load_forecast_data): SP_10M/DD_10M are computed from U_10M/V_10M; TOT_PREC6 is
     disaggregated from cumulative TOT_PREC; plain TOT_PREC is returned as-is.
     Returns a dataset with 'time' dimension (valid datetimes).
 
@@ -867,6 +884,7 @@ def _load_INCA_baseline_from_netcdf(
                    CLCT      total cloud cover         %        1h/10min   CT        10min    %         2022
                    U_10M     10 m zonal wind           m/s      1h/10min   derived from DD_10M, FF_10M
                    V_10M     10 m meridional wind      m/s      1h/10min   derived from DD_10M, FF_10M
+                   RELHUM_2M 2 m relative humidity     %        1h/10min   derived from T_2M, TD_2M
 
                  U_10M and V_10M use the meteorological convention: DD is
                  the direction the wind blows FROM, clockwise from North.
@@ -1042,7 +1060,11 @@ def _load_shifted(param: str, prefix: str) -> xr.DataArray:
             "TOT_PREC": "RP",
         },
     }
-    DERIVED_DEPS = {"U_10M": ["DD_10M", "FF_10M"], "V_10M": ["DD_10M", "FF_10M"]}
+    DERIVED_DEPS = {
+        "U_10M": ["DD_10M", "FF_10M"],
+        "V_10M": ["DD_10M", "FF_10M"],
+        "RELHUM_2M": ["T_2M", "TD_2M"],
+    }
     PARAM_UNITS = {
         "T_2M": "K",
         "TD_2M": "K",
@@ -1054,6 +1076,7 @@ def _load_shifted(param: str, prefix: str) -> xr.DataArray:
         "VMAX_10M": "m/s",
         "U_10M": "m/s",
         "V_10M": "m/s",
+        "RELHUM_2M": "%",
     }
     FREQ_TO_TD = {
         "1h": np.timedelta64(1, "h"),
@@ -1211,6 +1234,14 @@ def _load_cumul(
         if "V_10M" in params:
             merged["V_10M"] = (-ff * np.cos(dd_rad)).assign_attrs(units="m/s")
 
+    if "RELHUM_2M" in params:
+        # Same formula as the anemoi-inference surface-diagnostics post-processor
+        # (anemoi_plugins_meteoswiss.transform.filters.surface_diagnostics), so
+        # INCA is verified against the ML forecaster on a like-for-like basis.
+        merged["RELHUM_2M"] = ekdt.relative_humidity_from_dewpoint(
+            merged["T_2M"], merged["TD_2M"]
+        ).assign_attrs(units="%")
+
     # Restructure to match the earthkit GRIB engine profile: `step` is the
     # lead-time dimension, `valid_time` and `forecast_reference_time` are coords.
     ref_time_np = np.datetime64(reftime, "ns")
@@ -1355,7 +1386,13 @@ def load_forecast_data(
     """Load forecast data from GRIB files or an ICON archive.
 
     Handles derived and aggregated params transparently:
-    - ``SP_10M`` is computed from ``U_10M`` / ``V_10M``
+    - For the ML inference GRIB output: ``SP_10M``/``DD_10M``/``RELHUM_2M`` are
+      expected to be produced directly by the anemoi-inference
+      surface-diagnostics post-processor filter and are
+      loaded as native fields, not recomputed here.
+    - For baselines (ICON archive, INCA): these archives have no such
+      filters, so ``SP_10M``/``DD_10M`` are computed from ``U_10M``/``V_10M``
+      (or, for INCA, read from its native FF/DD product) as before.
     - ``TOT_PREC6`` is disaggregated from the cumulative ``TOT_PREC`` field
     - Plain ``TOT_PREC`` is returned as cumulative-from-start without disaggregation
 
@@ -1366,13 +1403,22 @@ def load_forecast_data(
     3. Otherwise → ICON operational archive (via :func:`_load_icon_baseline_from_grib`)
     """
     root = Path(root)
-    load_params = get_base_params(params)
     load_steps = get_steps(steps, params)
     if any(root.glob("*.grib")):
         LOG.info("Loading forecasts from GRIB files...")
+        # ML inference output: request params as-is (only stripping
+        # aggregation suffixes, e.g. TOT_PREC6 -> TOT_PREC) so fields the
+        # anemoi-inference filters already produced (SP_10M, DD_10M,
+        # RELHUM_2M) are read natively instead of being recomputed from
+        # their components. _disaggregated_and_derived_params below only
+        # falls back to compute_derived for a DERIVED_PARAMS entry that
+        # isn't already present in the loaded dataset.
+        ml_load_params = list(
+            dict.fromkeys(parse_aggregated_param(p)[0] for p in params)
+        )
         ds = _load_forecast_data_from_grib(
             files=_collect_ml_grib_files(root, load_steps),
-            params=load_params,
+            params=ml_load_params,
         )
         # Try to derive elevation from surface geopotential (FIS/z at step 0)
         # before falling back to ICON grid constants lookup.
@@ -1388,7 +1434,8 @@ def load_forecast_data(
             ds = _try_assign_elevation(ds)
     elif "INCA" in root.parts:
         LOG.info("Loading INCA baseline from NetCDF files...")
-        # INCA provides wind speed natively (FF), so don't expand SP_10M → U_10M/V_10M.
+        # INCA provides wind speed/direction natively (FF/DD), so don't expand
+        # SP_10M/DD_10M → U_10M/V_10M.
         # Only expand aggregated params (e.g. TOT_PREC6 → TOT_PREC).
         inca_load_params = list(
             dict.fromkeys(parse_aggregated_param(p)[0] for p in params)
@@ -1398,7 +1445,10 @@ def load_forecast_data(
         )
     else:
         LOG.info("Loading baseline forecasts from ICON GRIB archive...")
+        # ICON's own archive has no post-processor filters, so SP_10M/DD_10M
+        # are expanded to U_10M/V_10M here and computed by
+        # _disaggregated_and_derived_params below, same as before.
         ds = _load_icon_baseline_from_grib(
-            root, reftime, load_steps, load_params, member=member
+            root, reftime, load_steps, get_base_params(params), member=member
         )
     return _disaggregated_and_derived_params(ds, steps, params)