Add RMSE scores based vertical pressure profile recipe

src/CSET/cset_workflow/includes/metplus_point_stat.cylc

@@ -103,7 +103,7 @@
     # Runs VerPy metloader utility to create the VerPy databases
         inherit = VERPY_METLOADER
         [[[environment]]]
-            VERPY_DIR = {{VERPY_DIR}}
+            VERPY_DIR = {{VERPY_DIR|default("")}}
             VER_METHOD = area
             STAT_TYPE = cnt
             STREAM = point_stat

src/CSET/cset_workflow/meta/verification/rose-meta.conf

@@ -54,36 +54,71 @@ sort-key=met5
 
 ###################################
 # Scores
-[Verification/Scores]
-ns=Verification/Scores
+[Verification/Scores spatial]
+ns=Verification/Scores spatial
 sort-key=scores
-title=Scores
+title=Scores spatial
 
 [template variables=SCORES_ALL]
-ns=Verification/Scores
+ns=Verification/Scores spatial
 description=If selected, all scores shown below will be calculated.
             NOTE scores currently only available for model-model intercomparisons.
 type=python_boolean
 trigger=template variables=SCORES_RMSE_SPATIAL: False;
         template variables=SCORES_RMSE_TIMESERIES: False;
 compulsory=true
-sort-key=scores1
+sort-key=scoresspatial1
 
 [template variables=SCORES_RMSE_SPATIAL]
-ns=Verification/Scores
+ns=Verification/Scores spatial
 description=Create spatially mapped RMSE plots for the specified surface fields.
             The analysis method(s) set in SPATIAL_SURFACE_FIELD_METHOD will be used.
             NOTE currently only available for model-model intercomparisons.
 type=python_boolean
 compulsory=true
-sort-key=scores2
+sort-key=scoresspatial2
+
+[Verification/Scores timeseries]
+ns=Verification/Scores timeseries
+sort-key=scorestimeseries
+title=Scores timeseries
 
+[Verification/Scores timeseries]
 [template variables=SCORES_RMSE_TIMESERIES]
-ns=Verification/Scores
+ns=Verification/Scores timeseries
 description=Create a timeseries of the RMSE for the specified surface fields.
             NOTE currently only available for model-model intercomparisons.
 type=python_boolean
 compulsory=true
+sort-key=scorestimeseries1
+
+[Verification/Scores vertical profiles]
+ns=Verification/Scores vertical profiles
+sort-key=scoresvertical
+title=Scores vertical profiles
+
+[Verification/Scores vertical profiles]
+[template variables=SCORES_RMSE_VERTICAL_PROFILES]
+ns=Verification/Scores vertical profiles
+description=Create vertical profile RMSE plots for the specified level fields.
+    RMSE profiles are calculated at each grid point for all timesteps within a case study at each pressure level.
+type=python_boolean
+compulsory=true
+sort-key=scoresprofile1
+
+[template variables=SCORES_RMSE_VERTICAL_PROFILES_SEQUENCE]
+ns=Verification/Scores vertical profiles
+description=Also produce one RMSE vertical profile per time step.
+help=When True, the RMSE is calculated at each grid point and level for every
+    time step, collapsing only the spatial (horizontal) dimensions. This produces
+    a sequence of vertical RMSE profile plots, one per time step, in addition to
+    the profile that is aggregated across the whole case study.
+
+    When False (default), only the aggregated vertical RMSE profile is produced.
+type=python_boolean
+compulsory=true
+#trigger=template variables=SCORES_RMSE_VERTICAL_PROFILES: True;
+sort-key=scoresprofile2
 sort-key=scores3
 
 [template variables=SCORES_CRPS_FOR_ENSEMBLE]

src/CSET/cset_workflow/rose-suite.conf.example

@@ -125,6 +125,8 @@ RUN_METPLUS_POINT_STAT=False
 SCORES_ALL=False
 SCORES_RMSE_SPATIAL=False
 SCORES_RMSE_TIMESERIES=False
+SCORES_RMSE_VERTICAL_PROFILES=False
+SCORES_RMSE_VERTICAL_PROFILES_SEQUENCE=False
 SCORES_CRPS_FOR_ENSEMBLE=False
 SATURATION_FRACTION=False
 SCREEN_LEVEL_TEMPERATURE_PROBABILITIES=False

src/CSET/cset_workflow/site/nci-gadi.cylc

@@ -23,7 +23,7 @@
             """
         [[[ environment ]]]
             PROJECT = {{ PROJECT }}
-            PYTHONPATH = "{{VERPY_DIR}}"
+            PYTHONPATH = "{{VERPY_DIR|default("")}}"
 
 {% if RUN_METPLUS_GRID_STAT|default(False) or RUN_METPLUS_POINT_STAT|default(False) %}
     [[METPLUS]]

src/CSET/loaders/verification.py

@@ -23,7 +23,10 @@ def load(conf: Config):
     """Yield recipes from the given workflow configuration."""
     # Load a list of model detail dictionaries.
     models = get_models(conf.asdict())
+    if not models:
+        return
     # Models are listed in order, so model 1 is the first element.
+    base_model = models[0]
 
     if conf.SCORES_RMSE_SPATIAL or conf.SCORES_ALL:
         base_model = models[0]
@@ -63,6 +66,25 @@ def load(conf: Config):
                 model_ids=[base_model["id"], model["id"]],
                 aggregation=False,
             )
+    # including AGGREGATION_MODE in the variables dictionary to allow for clearer labeling of plots.
+    if conf.SCORES_RMSE_VERTICAL_PROFILES:
+        for model, field in itertools.product(models[1:], conf.PRESSURE_LEVEL_FIELDS):
+            yield RawRecipe(
+                recipe="generic_level_rmse_scores_profile.yaml",
+                variables={
+                    "VARNAME": field,
+                    "BASE_MODEL": base_model["name"],
+                    "OTHER_MODEL": model["name"],
+                    "PRESERVED_COORDS": ["pressure"],
+                    "AGGREGATION_MODE": "Case-study RMSE",
+                    "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
+                    "SUBAREA_EXTENT": conf.SUBAREA_EXTENT
+                    if conf.SELECT_SUBAREA
+                    else None,
+                },
+                model_ids=[base_model["id"], model["id"]],
+                aggregation=False,
+            )
 
     if conf.SCORES_CRPS_FOR_ENSEMBLE:
         for model, field in itertools.product(models, conf.SURFACE_FIELDS):
@@ -80,3 +102,22 @@ def load(conf: Config):
                 model_ids=[model["id"]],
                 aggregation=False,
             )
+    # including AGGREGATION_MODE in the variables dictionary to allow for clearer labeling of plots.
+    if conf.SCORES_RMSE_VERTICAL_PROFILES_SEQUENCE:
+        for model, field in itertools.product(models[1:], conf.PRESSURE_LEVEL_FIELDS):
+            yield RawRecipe(
+                recipe="generic_level_rmse_scores_profile.yaml",
+                variables={
+                    "VARNAME": field,
+                    "BASE_MODEL": base_model["name"],
+                    "OTHER_MODEL": model["name"],
+                    "PRESERVED_COORDS": ["time", "pressure"],
+                    "AGGREGATION_MODE": "Time-step RMSE",
+                    "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
+                    "SUBAREA_EXTENT": conf.SUBAREA_EXTENT
+                    if conf.SELECT_SUBAREA
+                    else None,
+                },
+                model_ids=[base_model["id"], model["id"]],
+                aggregation=False,
+            )

src/CSET/operators/plot.py

@@ -2270,28 +2270,74 @@ def plot_vertical_line_series(
         vmin = min(x_levels)
         vmax = max(x_levels)
 
-    # Matching the slices (matching by seq coord point; it may happen that
-    # evaluated models do not cover the same seq coord range, hence matching
-    # necessary)
-    cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
+    # Check if the cube has a sequence coordinate (e.g. time). If not, plot
+    # a single profile directly without iterating over a sequence.
+    sequence_coords = [
+        cube.coord(sequence_coordinate)
+        for cube in cubes
+        if cube.coords(sequence_coordinate)
+    ]
+    has_sequence_coord = len(sequence_coords) == len(cubes) and all(
+        np.size(coord.points) > 1 for coord in sequence_coords
+    )
+    has_scalar_sequence_coord = len(sequence_coords) == len(cubes) and all(
+        np.size(coord.points) == 1 for coord in sequence_coords
+    )
 
-    # Create a plot for each value of the sequence coordinate.
-    # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
-    # similar sequence values. Passing a CubeList into the internal plotting function
-    # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
-    # or an iris.cube.CubeList.
     plot_index = []
-    nplot = np.size(cubes[0].coord(sequence_coordinate).points)
-    for cubes_slice in cube_iterables:
-        # Format the coordinate value in a unit appropriate way.
-        seq_coord = cubes_slice[0].coord(sequence_coordinate)
+    if has_sequence_coord:
+        # Matching the slices (matching by seq coord point; it may happen that
+        # evaluated models do not cover the same seq coord range, hence matching
+        # necessary)
+        cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
+        nplot = np.size(cubes[0].coord(sequence_coordinate).points)
+        for cubes_slice in cube_iterables:
+            # Format the coordinate value in a unit appropriate way.
+            seq_coord = cubes_slice[0].coord(sequence_coordinate)
+            plot_title, plot_filename = _set_title_and_filename(
+                seq_coord, nplot, recipe_title, filename
+            )
+
+            # Do the actual plotting.
+            _plot_and_save_vertical_line_series(
+                cubes_slice,
+                coords,
+                "realization",
+                plot_filename,
+                series_coordinate,
+                title=plot_title,
+                vmin=vmin,
+                vmax=vmax,
+            )
+            plot_index.append(plot_filename)
+    elif has_scalar_sequence_coord:
+        # Scalar sequence coordinate (typically aggregated time bounds):
+        # make one plot and include sequence period in title/filename.
         plot_title, plot_filename = _set_title_and_filename(
-            seq_coord, nplot, recipe_title, filename
+            sequence_coords[0], 1, recipe_title, filename
         )
 
-        # Do the actual plotting.
         _plot_and_save_vertical_line_series(
-            cubes_slice,
+            cubes,
+            coords,
+            "realization",
+            plot_filename,
+            series_coordinate,
+            title=plot_title,
+            vmin=vmin,
+            vmax=vmax,
+        )
+        plot_index.append(plot_filename)
+    else:
+        # 1D case: no sequence coordinate, plot a single profile.
+        plot_title = recipe_title
+        if filename:
+            plot_filename = filename
+        else:
+            plot_filename = f"{slugify(plot_title)}.png"
+
+        _plot_and_save_vertical_line_series(
+            cubes,
             coords,
             "realization",
             plot_filename,

src/CSET/operators/scoreswrappers.py

@@ -146,6 +146,53 @@ def _sort_cubes_for_verification(cubes: CubeList):
     return base, other
 
 
+def _resolve_preserve_dims(
+    cube: Cube,
+    data_array: xr.DataArray,
+    preserved_coordinates: list[str] | str | None,
+) -> list[str] | None:
+    """Resolve preserve coordinates to xarray dimension names.
+
+    The ``scores`` package expects preserve dimensions to match xarray
+    dimension names. In Iris data, commonly used coordinates such as ``time``
+    may be auxiliary coordinates attached to a differently named dimension
+    (e.g. ``dim0``). This helper maps coordinate names to their underlying
+    dimension names.
+    """
+    if preserved_coordinates is None:
+        return None
+
+    coord_names = (
+        [preserved_coordinates]
+        if isinstance(preserved_coordinates, str)
+        else preserved_coordinates
+    )
+    preserve_dims: list[str] = []
+
+    for coord_name in coord_names:
+        # Already an xarray dimension name.
+        if coord_name in data_array.dims:
+            if coord_name not in preserve_dims:
+                preserve_dims.append(coord_name)
+            continue
+
+        # Otherwise, map coordinate name to dimension index/indices.
+        try:
+            dim_indices = cube.coord_dims(coord_name)
+        except iris.exceptions.CoordinateNotFoundError:
+            # Keep original name so scores raises a clear error for unknown keys.
+            if coord_name not in preserve_dims:
+                preserve_dims.append(coord_name)
+            continue
+
+        for dim_index in dim_indices:
+            dim_name = data_array.dims[dim_index]
+            if dim_name not in preserve_dims:
+                preserve_dims.append(dim_name)
+
+    return preserve_dims
+
+
 def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
     r"""Calculate the Root Mean Square Error (RMSE) using scores.
 
@@ -160,7 +207,7 @@ def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None =
         A CubeList containing exactly two cubes: a base and an "other" model,
         this can be an analysis and the model.
     preserved_coordinates: list[str] | str | None, default is None.
-        The coordinates that you wish to preserve in the calculaiton of the
+        The coordinates (or xarray dimension names) that you wish to preserve in the calculaiton of the
         RMSE. For example if you want a map of each time you can preserve
         ["time","grid_latitude", "grid_longitude"] or if you want a time series
         you can preserve ["time"], if you want to collapse to a single value
@@ -187,15 +234,49 @@ def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None =
         forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
     """
     base, other = _sort_cubes_for_verification(cubes)
+
+    other_xr = xr.DataArray.from_iris(other)
+    base_xr = xr.DataArray.from_iris(base)
+    preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
+
     # Scores operators on xarray data arrays, so we transform the iris cube into an array,
     # apply scores, and then transform it back.
     RMSE = xr.DataArray.to_iris(
         scores.continuous.rmse(
-            xr.DataArray.from_iris(other),
-            xr.DataArray.from_iris(base),
-            preserve_dims=preserved_coordinates,
+            other_xr,
+            base_xr,
+            preserve_dims=preserve_dims,
         )
     )
+
+    # If time is aggregated out, attach a scalar time coordinate with bounds
+    # so plotting can display the aggregated period in the title.
+    try:
+        if not RMSE.coords("time"):
+            base_time = base.coord("time")
+            time_vals = (
+                base_time.bounds.flatten()
+                if base_time.has_bounds()
+                else base_time.points
+            )
+            t_start = float(time_vals[0])
+            t_end = float(time_vals[-1])
+            t_mid = 0.5 * (t_start + t_end)
+
+            RMSE.add_aux_coord(
+                iris.coords.AuxCoord(
+                    t_mid,
+                    standard_name=base_time.standard_name,
+                    long_name=base_time.long_name,
+                    var_name=base_time.var_name,
+                    units=base_time.units,
+                    bounds=np.array([t_start, t_end]),
+                    attributes=base_time.attributes.copy(),
+                )
+            )
+    except iris.exceptions.CoordinateNotFoundError:
+        pass
+
     RMSE.rename(f"RMSE_of_{base.name()}")
     return RMSE