Add RMSE scores based vertical pressure profile recipe

src/CSET/cset_workflow/includes/metplus_point_stat.cylc

@@ -80,7 +80,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

@@ -77,7 +77,7 @@ trigger=template variables=SCORES_SPATIAL_RMSE: False;
         template variables=SCORES_TIMESERIES_AB: False;
         template variables=SCORES_TIMESERIES_PC: False;
 compulsory=true
-sort-key=scores1
+sort-key=scoresspatial1
 
 
 # Scores spatial plots
@@ -226,3 +226,31 @@ help=Usually this should be set to 0, but in some case it may be 1.
 type=integer
 compulsory=true
 sort-key=scores4b
+
+[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

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

@@ -150,6 +150,8 @@ SCORES_ALL=False
 SCORES_CRPS_FOR_ENSEMBLE=False
 SCORES_RMSE_SPATIAL=False
 SCORES_RMSE_TIMESERIES=False
+SCORES_RMSE_VERTICAL_PROFILES_SEQUENCE=False
+SCORES_RMSE_VERTICAL_PROFILES=False
 SCORES_SPATIAL_AB=False
 SCORES_SPATIAL_MAE=False
 SCORES_SPATIAL_RMSE=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

@@ -49,12 +49,14 @@ 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]
 
     scores_spatial_methods = _get_scores_spatial_methods(conf)
     if scores_spatial_methods:
         # Produce 2D spatial plots of scores metrics.
-        base_model = models[0]
         for model, field, method, scores_method in itertools.product(
             models[1:],
             conf.SURFACE_FIELDS,
@@ -101,7 +103,6 @@ def load(conf: Config):
     scores_timeseries_methods = _get_scores_timeseries_methods(conf)
     if scores_timeseries_methods:
         # Produce timeseries plots of scores metrics averaged over the domain for each case study.
-        base_model = models[0]
         for model, field, scores_method in itertools.product(
             models[1:], conf.SURFACE_FIELDS, scores_timeseries_methods
         ):
@@ -137,3 +138,43 @@ 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:
+        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,
+            )
+
+    # 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

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

@@ -208,7 +208,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

src/CSET/recipes/verification/generic_level_rmse_scores_profile.yaml

@@ -0,0 +1,48 @@
+category: Scores
+title: "$VARNAME\nRMSE vertical profile between $OTHER_MODEL and $BASE_MODEL ($AGGREGATION_MODE)"
+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 for each pressure
+  level field utilising the preserved coordinates in the scores package.
+
+  By adding a preserved coordinate to the recipe, the RMSE can be calculated for each
+  time step in a case study or for all time steps combined.
+
+  References:
+  https://scores.readthedocs.io/en/stable/
+  https://scores.readthedocs.io/en/stable/api.html#scores.continuous.rmse
+
+
+  Aggregation mode:
+  - Case-study RMSE: Single RMSE profile calculated using all time points in the case study.
+  - Time-step RMSE: RMSE profiles calculated at each time point in a case study.
+
+  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