2229 scores rmse vertical profile pressure levels add ensemble capability

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

@@ -73,3 +73,25 @@ description=Create a timeseries of the RMSE for the specified surface fields.
 type=python_boolean
 compulsory=true
 sort-key=scores2
+
+[template variables=SCORES_RMSE_VERTICAL_PROFILES]
+ns=Verification/Scores
+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=scores3
+
+[template variables=SCORES_RMSE_VERTICAL_PROFILES_SEQUENCE]
+ns=Verification/Scores
+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=scores4

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

@@ -127,6 +127,8 @@ RAIN_PRESENCE_SPATIAL_DIFFERENCE=False
 RAIN_PRESENCE_SPATIAL_PLOT=False
 SCORES_RMSE_SPATIAL=False
 SCORES_RMSE_TIMESERIES=False
+SCORES_RMSE_VERTICAL_PROFILES=True
+SCORES_RMSE_VERTICAL_PROFILES_SEQUENCE=False
 SCREEN_LEVEL_TEMPERATURE_PROBABILITIES=False
 !!SCREEN_LEVEL_TEMPERATURE_SPATIAL_PROBABILITY_WITHOUT_CONTROL_MEMBER=False
 SELECT_SUBAREA=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

@@ -24,9 +24,9 @@ def load(conf: Config):
     # Load a list of model detail dictionaries.
     models = get_models(conf.asdict())
     # Models are listed in order, so model 1 is the first element.
+    base_model = models[0]
 
     if conf.SCORES_RMSE_SPATIAL:
-        base_model = models[0]
         for model, field, method in itertools.product(
             models[1:], conf.SURFACE_FIELDS, conf.SPATIAL_SURFACE_FIELD_METHOD
         ):
@@ -47,7 +47,6 @@ def load(conf: Config):
             )
 
     if conf.SCORES_RMSE_TIMESERIES:
-        base_model = models[0]
         for model, field in itertools.product(models[1:], conf.SURFACE_FIELDS):
             yield RawRecipe(
                 recipe="timeseries_surface_rmse_scores.yaml",
@@ -63,3 +62,41 @@ def load(conf: Config):
                 model_ids=[base_model["id"], model["id"]],
                 aggregation=False,
             )
+
+    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", "realization"],
+                    "AGGREGATION_MODE": "pressure",
+                    "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_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", "realization"],
+                    "AGGREGATION_MODE": "timeseries",
+                    "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

@@ -2181,28 +2181,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

@@ -140,6 +140,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.
 
@@ -154,7 +201,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
@@ -181,14 +228,48 @@ 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

src/CSET/recipes/verification/generic_level_rmse_scores_profile.yaml

@@ -0,0 +1,40 @@
+category: Scores
+title: "$VARNAME\nRMSE vertical profile ($AGGREGATION_MODE) between $OTHER_MODEL and $BASE_MODEL"
+description: |
+  Extracts and plots the Root Mean Square Error (RMSE) of $VARNAME
+  as a vertical profile on pressure levels. The RMSE is calculated
+  pairwise between $BASE_MODEL and $OTHER_MODEL at each grid point
+  and vertical level, then aggregated over horizontal grid points.
+
+  Aggregation mode:
+  - pressure: Aggregates over the entire time period, producing one RMSE profile
+  - timeseries: Preserves time, producing one RMSE profile per time step
+
+  A larger RMSE implies a greater error than a smaller RMSE. An
+  RMSE of zero implies the two fields match exactly.
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    model_names: [$BASE_MODEL, $OTHER_MODEL]
+    constraint:
+      operator: constraints.combine_constraints
+      variable_constraint:
+        operator: constraints.generate_var_constraint
+        varname: $VARNAME
+      level_constraint:
+        operator: constraints.generate_level_constraint
+        coordinate: pressure
+        levels: "*"
+    subarea_type: $SUBAREA_TYPE
+    subarea_extent: $SUBAREA_EXTENT
+
+  - operator: scoreswrappers.scores_rmse
+    preserved_coordinates: $PRESERVED_COORDS
+
+  - operator: plot.plot_vertical_line_series
+    series_coordinate: pressure
+    sequence_coordinate: time
+
+  - operator: write.write_cube_to_nc
+    overwrite: True