Add summary statistic (mean, max, min,...) plots and introduce additional aggregation options

src/CSET/operators/collapse.py

@@ -69,6 +69,8 @@ def collapse(
     ValueError
         If additional_percent wasn't supplied while using PERCENTILE method.
     """
+    if method == "SEQ" or method == "" or method is None:
+        return cubes
     if method == "PERCENTILE" and additional_percent is None:
         raise ValueError("Must specify additional_percent")
     collapsed_cubes = iris.cube.CubeList([])
@@ -88,6 +90,10 @@ def collapse(
                         percent=additional_percent,
                     )
                 )
+            elif method == "RANGE":
+                cube_max = cube.collapsed(coordinate, iris.analysis.MAX)
+                cube_min = cube.collapsed(coordinate, iris.analysis.MIN)
+                collapsed_cubes.append(cube_max - cube_min)
             else:
                 collapsed_cubes.append(
                     cube.collapsed(coordinate, getattr(iris.analysis, method))

src/CSET/operators/plot.py

@@ -325,7 +325,6 @@ def _colorbar_map_levels(cube: iris.cube.Cube, axis: Literal["x", "y"] | None =
             norm = mpl.colors.BoundaryNorm(levels, ncolors=cmap.N)
             logging.debug("Using levels for %s colorbar.", varname)
             logging.info("Using levels: %s", levels)
-
         except KeyError:
             # Get the range for this variable.
             vmin, vmax = var_colorbar["min"], var_colorbar["max"]
@@ -529,6 +528,11 @@ def _plot_and_save_postage_stamp_spatial_plot(
             else:
                 raise TypeError("Unknown vmin and vmax range.")
                 vmin, vmax = None, None
+            # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
+            # if levels are defined.
+            if norm is not None:
+                vmin = None
+                vmax = None
             # pcolormesh plot of the field.
             plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
         else:
@@ -597,6 +601,9 @@ def _plot_and_save_line_series(
     # Store min/max ranges.
     y_levels = []
 
+    # Check match-up across sequence coords gives consistent sizes
+    _validate_cubes_coords(cubes, coords)
+
     for cube, coord in zip(cubes, coords, strict=True):
         label = None
         color = "black"
@@ -681,6 +688,9 @@ def _plot_and_save_vertical_line_series(
 
     model_colors_map = _get_model_colors_map(cubes)
 
+    # Check match-up across sequence coords gives consistent sizes
+    _validate_cubes_coords(cubes, coords)
+
     for cube, coord in zip(cubes, coords, strict=True):
         label = None
         color = "black"
@@ -1076,13 +1086,18 @@ def _spatial_plot(
 
     # Create a plot for each value of the sequence coordinate.
     plot_index = []
+    nplot = np.size(cube.coord(sequence_coordinate).points)
     for cube_slice in cube.slices_over(sequence_coordinate):
         # Use sequence value so multiple sequences can merge.
         sequence_value = cube_slice.coord(sequence_coordinate).points[0]
         plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
         coord = cube_slice.coord(sequence_coordinate)
         # Format the coordinate value in a unit appropriate way.
-        title = f"{recipe_title}\n{coord.units.title(coord.points[0])}"
+        title = f"{recipe_title}\n [{coord.units.title(coord.points[0])}]"
+        # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
+        if nplot == 1 and coord.has_bounds:
+            if np.size(coord.bounds) > 1:
+                title = f"{recipe_title}\n [{coord.units.title(coord.bounds[0][0])} to {coord.units.title(coord.bounds[0][1])}]"
         # Do the actual plotting.
         plotting_func(
             cube_slice,
@@ -1201,6 +1216,19 @@ def _validate_cube_shape(
         )
 
 
+def _validate_cubes_coords(
+    cubes: iris.cube.CubeList, coords: list[iris.coords.Coord]
+) -> None:
+    """Check same number of cubes as sequence coordinate for zip functions."""
+    if len(cubes) != len(coords):
+        raise ValueError(
+            f"The number of CubeList entries ({len(cubes)}) should equal the number "
+            f"of sequence coordinates ({len(coords)})."
+            f"Check that number of time entries in input data are consistent if "
+            f"performing time-averaging steps prior to plotting outputs."
+        )
+
+
 ####################
 # Public functions #
 ####################
@@ -1514,13 +1542,18 @@ def filter_cube_iterables(cube_iterables) -> bool:
     # 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:
         # Use sequence value so multiple sequences can merge.
         seq_coord = cubes_slice[0].coord(sequence_coordinate)
         sequence_value = seq_coord.points[0]
         plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
         # Format the coordinate value in a unit appropriate way.
-        title = f"{recipe_title}\n{seq_coord.units.title(sequence_value)}"
+        title = f"{recipe_title}\n [{seq_coord.units.title(sequence_value)}]"
+        # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
+        if nplot == 1 and seq_coord.has_bounds:
+            if np.size(seq_coord.bounds) > 1:
+                title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
         # Do the actual plotting.
         _plot_and_save_vertical_line_series(
             cubes_slice,
@@ -1778,6 +1811,7 @@ def plot_histogram_series(
         ]
 
     plot_index = []
+    nplot = np.size(cube.coord(sequence_coordinate).points)
     # 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
@@ -1793,7 +1827,11 @@ def plot_histogram_series(
         plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
         coord = single_cube.coord(sequence_coordinate)
         # Format the coordinate value in a unit appropriate way.
-        title = f"{recipe_title}\n{coord.units.title(coord.points[0])}"
+        title = f"{recipe_title}\n [{coord.units.title(coord.points[0])}]"
+        # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
+        if nplot == 1 and coord.has_bounds:
+            if np.size(coord.bounds) > 1:
+                title = f"{recipe_title}\n [{coord.units.title(coord.bounds[0][0])} to {coord.units.title(coord.bounds[0][1])}]"
         # Do the actual plotting.
         plotting_func(
             cube_slice,

src/CSET/operators/read.py

@@ -71,7 +71,7 @@ def read_cube(
         Constraints to filter data by. Defaults to unconstrained.
     model_names: str | list[str], optional
         Names of the models that correspond to respective paths in file_paths.
-    subarea_type: "realworld" | "modelrelative", optional
+    subarea_type: "gridcells" | "modelrelative" | "realworld", optional
         Whether to constrain data by model relative coordinates or real world
         coordinates.
     subarea_extent: list, optional
@@ -187,7 +187,6 @@ def read_cubes(
 
     # Select sub region.
     cubes = _cutout_cubes(cubes, subarea_type, subarea_extent)
-
     # Merge and concatenate cubes now metadata has been fixed.
     cubes = cubes.merge()
     cubes = cubes.concatenate()
@@ -540,7 +539,7 @@ def _longitude_fix_callback(cube: iris.cube.Cube):
         long_points -= 180.0
     long_coord.points = long_points
 
-    # Update coord bounds to be consistent with wrapping
+    # Update coord bounds to be consistent with wrapping.
     if long_coord.has_bounds() and np.size(long_coord) > 1:
         long_coord.bounds = None
         long_coord.guess_bounds()
@@ -685,6 +684,7 @@ def _fix_um_radtime(cube: iris.cube.Cube):
     """
     try:
         if cube.attributes["STASH"] in [
+            "m01s01i207",
             "m01s01i208",
             "m01s02i205",
             "m01s02i201",
@@ -697,7 +697,6 @@ def _fix_um_radtime(cube: iris.cube.Cube):
             # Convert time points to datetime objects
             time_unit = time_coord.units
             time_points = time_unit.num2date(time_coord.points)
-
             # Skip if times don't need fixing.
             if time_points[0].minute == 0 and time_points[0].second == 0:
                 return
@@ -869,6 +868,7 @@ def _lfric_time_callback(cube: iris.cube.Cube):
     # Construct forecast_reference time if it doesn't exist.
     try:
         tcoord = cube.coord("time")
+        # Set time coordinate to common basis "hours since 1970"
         try:
             tcoord.convert_units("hours since 1970-01-01 00:00:00")
         except ValueError:
@@ -941,3 +941,30 @@ def _lfric_forecast_period_standard_name_callback(cube: iris.cube.Cube):
             coord.standard_name = "forecast_period"
     except iris.exceptions.CoordinateNotFoundError:
         pass
+
+
+def _remove_time0(cubes: iris.cube.CubeList):
+    """Remove T0 from UM inputs to allow time-averaged comparison with LFRic.
+
+    A number of UM outputs contain T=0 initial time diagnostic fields, while
+    LFRic diagnostics begin from T=1 output step. This does not cause issues
+    for comparing UM with LFRic for hour-by-hour comparisons, for which times
+    are matched up in code. However, for any recipes requiring collapse by
+    time ahead of comparison, computing averages over e.g. 24h vs 25h window
+    results in different timestamps in each collapsed cube, breaking subsequent
+    CSET time-checking steps to compare like-with-like.
+
+    This function removes any outputs at T=0 to support time-processed comparisons.
+    """
+    valid_cubes = iris.cube.CubeList()
+    for cube in cubes:
+        if cube.coords("forecast_period"):
+            valid_cube = cube.extract(
+                iris.Constraint(forecast_period=lambda cell: cell >= 0.25)
+            )
+            if valid_cube:
+                valid_cubes.append(valid_cube)
+        else:
+            valid_cubes.append(cube)
+
+    return valid_cubes

src/CSET/recipes/generic_level_domain_mean_time_series_case_aggregation_all.yaml

@@ -0,0 +1,34 @@
+category: $LEVELTYPE Level Time Series
+title: "$VARNAME $LEVEL domain mean time series\n Aggregation over all cases"
+description: Plots a time series of the domain mean $VARNAME at $LEVELTYPE $LEVEL by using all available times.
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    model_names: $MODEL_NAME
+    constraint:
+      operator: constraints.combine_constraints
+      varname_constraint:
+        operator: constraints.generate_var_constraint
+        varname: $VARNAME
+      level_constraint:
+        operator: constraints.generate_level_constraint
+        coordinate: $LEVELTYPE
+        levels: $LEVEL
+
+  - operator: collapse.collapse
+    coordinate: [grid_latitude, grid_longitude]
+    method: MEAN
+
+  - operator: aggregate.ensure_aggregatable_across_cases
+
+  - operator: collapse.collapse_by_validity_time
+    method: SEQ
+
+  # Make a single NetCDF with all the data inside it.
+  - operator: write.write_cube_to_nc
+    overwrite: True
+
+  # Plot the data.
+  - operator: plot.plot_line_series
+    series_coordinate: time

src/CSET/recipes/generic_level_domain_mean_vertical_profile_series.yaml

@@ -1,6 +1,7 @@
 category: Profiles
-title: Domain mean $VARNAME vertical profile
-description: Plots a time series of vertical profiles for the domain horizontal mean $VARNAME using a $LEVELTYPE coordinate.
+title: Domain mean $VARNAME vertical profile $METHOD
+description: Plots the $METHOD of vertical profiles for the domain horizontal mean
+  $VARNAME using a $LEVELTYPE coordinate.
 
 steps:
   - operator: read.read_cubes
@@ -22,6 +23,11 @@ steps:
     coordinate: [grid_latitude, grid_longitude]
     method: MEAN
 
+  # If required, collapse data using specified time analysis method (or sequence)
+  - operator: collapse.collapse
+    coordinate: [time]
+    method: $METHOD
+
   # Plot the vertical line series.
   - operator: plot.plot_vertical_line_series
     series_coordinate: $LEVELTYPE

src/CSET/recipes/generic_level_domain_mean_vertical_profile_series_case_aggregation_all.yaml

@@ -0,0 +1,38 @@
+category: Profiles
+title: "Domain mean $VARNAME vertical profile\n Aggregation over all cases"
+description: |
+  Plots a vertical profile for the domain horizontal mean $VARNAME using a
+  $LEVELTYPE coordinate. Multiple cases are aggregated together for all available times.
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    model_names: $MODEL_NAME
+    constraint:
+      operator: constraints.combine_constraints
+      variable_constraint:
+        operator: constraints.generate_var_constraint
+        varname: $VARNAME
+      level_constraint:
+        operator: constraints.generate_level_constraint
+        coordinate: $LEVELTYPE
+        levels: "*"
+
+  - operator: collapse.collapse
+    coordinate: [grid_latitude, grid_longitude]
+    method: MEAN
+
+  - operator: aggregate.ensure_aggregatable_across_cases
+
+  - operator: collapse.collapse
+    coordinate: [time]
+    method: MEAN
+
+  # Plot the vertical line series.
+  - operator: plot.plot_vertical_line_series
+    series_coordinate: $LEVELTYPE
+    sequence_coordinate: time
+
+  # Make a single NetCDF with all the data inside it.
+  - operator: write.write_cube_to_nc
+    overwrite: True

src/CSET/recipes/generic_level_histogram_series.yaml

@@ -30,7 +30,7 @@ steps:
     subarea_extent: $SUBAREA_EXTENT
 
   - operator: plot.plot_histogram_series
-    sequence_coordinate: time
+    sequence_coordinate: $SEQUENCE
     # stamp_coordinate and single_plot optional and only required for ensemble data
     stamp_coordinate: "realization"
     single_plot: False

src/CSET/recipes/generic_level_histogram_series_case_aggregation_all.yaml

@@ -1,5 +1,5 @@
 category: Histogram
-title: $VARNAME $LEVELTYPE $LEVEL histogram case aggregation
+title: "Histogram $VARNAME $LEVELTYPE $LEVEL\n Aggregation over all cases"
 description: |
   Extracts and plots the probability density of $LEVEL level $VARNAME from a
   file at $LEVELTYPE level $LEVEL of $MODEL_NAME. It uses

src/CSET/recipes/generic_level_spatial_plot_sequence.yaml

@@ -1,7 +1,7 @@
 category: $LEVELTYPE Level Spatial
-title: $MODEL_NAME $VARNAME $LEVELTYPE$LEVEL
+title: $MODEL_NAME $VARNAME $LEVELTYPE$LEVEL $METHOD
 description: |
-  Extracts ands plots the $VARNAME from $MODEL_NAME at $LEVELTYPE level $LEVEL.
+  Extracts and plots the $METHOD of $VARNAME from $MODEL_NAME at $LEVELTYPE level $LEVEL.
 
 steps:
   - operator: read.read_cube
@@ -18,6 +18,10 @@ steps:
     subarea_type: $SUBAREA_TYPE
     subarea_extent: $SUBAREA_EXTENT
 
+  - operator: collapse.collapse
+    coordinate: [time]
+    method: $METHOD
+
   - operator: plot.spatial_contour_plot
     sequence_coordinate: time
 

src/CSET/recipes/generic_level_spatial_plot_sequence_case_aggregation_mean_all.yaml

@@ -0,0 +1,30 @@
+category: $LEVELTYPE Level Spatial
+title: "$MODEL_NAME $VARNAME $LEVELTYPE$LEVEL\n Aggregation over all cases"
+description: |
+  Extracts ands plots the $VARNAME from $MODEL_NAME at $LEVELTYPE level $LEVEL.
+  Multiple cases are aggregated together for all available times.
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    constraint:
+      operator: constraints.combine_constraints
+      variable_constraint:
+        operator: constraints.generate_var_constraint
+        varname: $VARNAME
+      level_constraint:
+        operator: constraints.generate_level_constraint
+        coordinate: $LEVELTYPE
+        levels: $LEVEL
+
+  - operator: aggregate.ensure_aggregatable_across_cases
+
+  - operator: collapse.collapse
+    coordinate: time
+    method: MEAN
+
+  - operator: plot.spatial_contour_plot
+    sequence_coordinate: time
+
+  - operator: write.write_cube_to_nc
+    overwrite: True

src/CSET/recipes/generic_surface_domain_mean_time_series.yaml

@@ -1,6 +1,6 @@
 category: Surface Time Series
 title: Domain mean $VARNAME time series
-description: Plots a time series of the domain mean $VARNAME.
+description: Plots a time series of the domain horizontal mean $VARNAME.
 
 steps:
   - operator: read.read_cubes