Add functionality to calculate and plot mean and rolling mean of spectra

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

@@ -229,6 +229,24 @@ type=python_boolean
 compulsory=true
 sort-key=c1surface8a
 
+[template variables=SPECTRUM_SURFACE_FIELD_AGGREGATION]
+ns=Diagnostics/Fields
+description=Create case aggregated power spectrum plots for surface.
+            Select all options required.
+       Option1: Mean over all times.
+       Option2: Rolling time mean.
+type=python_boolean,python_boolean
+compulsory=true
+sort-key=c1surface8b
+
+[template variables=WINDOW_LEN_SURFACE]
+ns=Diagnostics/Fields
+description= The length of the window (in hours) that the rolling mean uses.
+help= A single value
+type=integer
+compulsory=true
+sort-key=c1surface8c
+
 
 #######################################################################
 # Pressure level fields.
@@ -463,6 +481,26 @@ type=python_boolean
 compulsory=true
 sort-key=c2pressure9a
 
+[template variables=SPECTRUM_PLEVEL_FIELD_AGGREGATION]
+ns=Diagnostics/Pressure
+description=Create case aggregated power spectrum plots for specified pressure levels.
+            Select all options required.
+       Option1: Mean over all times.
+       Option2: Rolling time mean.
+type=python_boolean,python_boolean
+compulsory=true
+sort-key=c2pressure9b
+
+[template variables=WINDOW_LEN_PLEVEL]
+ns=Diagnostics/Pressure
+description= The length of the window (in hours) that the rolling mean uses.
+help= A single value
+type=integer
+compulsory=true
+sort-key=c2pressure9c
+
+
+
 [template variables=SPATIAL_STRUCTURAL_SIMILARITY_PLEVEL_FIELD]
 ns=Diagnostics/Pressure
 description=Create spatially mapped structural similarity plots for
@@ -683,6 +721,26 @@ type=python_boolean
 compulsory=true
 sort-key=c3modellevel9a
 
+[template variables=SPECTRUM_MLEVEL_FIELD_AGGREGATION]
+ns=Diagnostics/ModelLevel
+description=Create case aggregated power spectrum plots for specified model levels.
+            Select all options required.
+       Option1: Mean over all times.
+       Option2: Rolling time mean.
+type=python_boolean,python_boolean
+compulsory=true
+sort-key=c3modellevel9b
+
+[template variables=WINDOW_LEN_MLEVEL]
+ns=Diagnostics/ModelLevel
+description= The length of the window (in hours) that the rolling mean uses.
+help= A single value
+type=integer
+compulsory=true
+sort-key=c3modellevel9c
+
+
+
 [template variables=SPATIAL_STRUCTURAL_SIMILARITY_MLEVEL]
 ns=Diagnostics/ModelLevel
 description=Create spatially mapped structural similarity plots for

src/CSET/loaders/power_spectrum.py

@@ -99,3 +99,102 @@ def load(conf: Config):
                 model_ids=[model["id"] for model in models],
                 aggregation=False,
             )
+
+    # Create a list of case aggregation types.
+    AGGREGATION_TYPES = ["all", "rolling"]
+    print("In ps loader ", AGGREGATION_TYPES, conf.SURFACE_FIELDS)
+
+    # Surface (2D) fields.
+    for atype, field in itertools.product(AGGREGATION_TYPES, conf.SURFACE_FIELDS):
+        #        if conf.SPECTRUM_SURFACE_FIELD_AGGREGATION[AGGREGATION_TYPES.index(atype)]:
+        index = AGGREGATION_TYPES.index(atype)
+        aggregations = conf.SPECTRUM_SURFACE_FIELD_AGGREGATION
+        print("ALL INFO ", atype, index, aggregations)
+        if len(aggregations) > index and aggregations[index]:
+            yield RawRecipe(
+                recipe=f"generic_surface_power_spectrum_series_mean_{atype}.yaml",
+                variables={
+                    "VARNAME": field,
+                    "MODEL_NAME": [model["name"] for model in models],
+                    #                    "SEQUENCE": "time"
+                    #                    if conf.SPECTRUM_SURFACE_FIELD_SEQUENCE
+                    #                    else "realization",
+                    "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
+                    "SUBAREA_EXTENT": conf.SUBAREA_EXTENT
+                    if conf.SELECT_SUBAREA
+                    else None,
+                    "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "",
+                    "SPECTRUM_SURFACE_FIELD_SEQUENCE": conf.SPECTRUM_SURFACE_FIELD_SEQUENCE,
+                    "WINDOW_LEN_SURFACE": conf.WINDOW_LEN_SURFACE
+                    if atype == "rolling"
+                    else None,
+                },
+                model_ids=[model["id"] for model in models],
+                aggregation=True,
+            )
+
+    # Pressure level fields.
+    for atype, field, plevel in itertools.product(
+        AGGREGATION_TYPES, conf.PRESSURE_LEVEL_FIELDS, conf.PRESSURE_LEVELS
+    ):
+        if conf.SPECTRUM_PLEVEL_FIELD_AGGREGATION[AGGREGATION_TYPES.index(atype)]:
+            # Build the variables dict *without* WINDOW_LEN_PLEVEL first
+            variables = {
+                "VARNAME": field,
+                "LEVELTYPE": "pressure",
+                "LEVEL": [plevel],
+                "MODEL_NAME": [model["name"] for model in models],
+                "SEQUENCE": "time"
+                if conf.SPECTRUM_PLEVEL_FIELD_SEQUENCE
+                else "realization",
+                "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
+                "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None,
+                "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "",
+                "SPECTRUM_PLEVEL_FIELD_SEQUENCE": conf.SPECTRUM_PLEVEL_FIELD_SEQUENCE,
+            }
+
+            # Add WINDOW_LEN_PLEVEL *only if* rolling
+            if atype == "rolling":
+                variables["WINDOW_LEN_PLEVEL"] = conf.WINDOW_LEN_PLEVEL
+            else:
+                variables["WINDOW_LEN_PLEVEL"] = None
+
+            yield RawRecipe(
+                recipe=f"generic_level_power_spectrum_series_plevel_mean_{atype}.yaml",
+                variables=variables,
+                model_ids=[model["id"] for model in models],
+                aggregation=True,
+            )
+
+    # Model level fields.
+    for atype, field, mlevel in itertools.product(
+        AGGREGATION_TYPES, conf.MODEL_LEVEL_FIELDS, conf.MODEL_LEVELS
+    ):
+        if conf.POWER_SPECTRUM_MLEVEL_FIELD_AGGREGATION[AGGREGATION_TYPES.index(atype)]:
+            # Build the variables dict *without* WINDOW_LEN_MLEVEL first
+            variables = {
+                "VARNAME": field,
+                "LEVELTYPE": "model_level_number",
+                "LEVEL": [mlevel],
+                "MODEL_NAME": [model["name"] for model in models],
+                "SEQUENCE": "time"
+                if conf.SPECTRUM_MLEVEL_FIELD_SEQUENCE
+                else "realization",
+                "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
+                "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None,
+                "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "",
+                "SPECTRUM_MLEVEL_FIELD_SEQUENCE": conf.SPECTRUM_MLEVEL_FIELD_SEQUENCE,
+            }
+
+            # Add WINDOW_LEN_MLEVEL *only if* rolling
+            if atype == "rolling":
+                variables["WINDOW_LEN_MLEVEL"] = conf.WINDOW_LEN_MLEVEL
+            else:
+                variables["WINDOW_LEN_MLEVEL"] = None
+
+            yield RawRecipe(
+                recipe=f"generic_level_power_spectrum_series_mlevel_mean_{atype}.yaml",
+                variables=variables,
+                model_ids=[model["id"] for model in models],
+                aggregation=True,
+            )

src/CSET/operators/aggregate.py

@@ -242,7 +242,7 @@ def add_hour_coordinate(
 
 
 def rolling_window_time_aggregation(
-    cubes: iris.cube.Cube | iris.cube.CubeList, method: str, window: int
+    cubes: iris.cube.Cube | iris.cube.CubeList, method: str, window_hours: int
 ) -> iris.cube.Cube | iris.cube.CubeList:
     """Aggregate a cube along the time dimension using a rolling window.
 
@@ -273,8 +273,28 @@ def rolling_window_time_aggregation(
     for cube in iter_maybe(cubes):
         # Use a rolling window in time to applied specified aggregation method
         # over a specified window length.
+        # Input window length in hours and convert to number of time points in
+        # the window.
+        # Add catches
+        if window_hours is None:
+            raise ValueError("ROLLING_MEAN requires kwarg: window_hours=<int>")
+        # Determine number of steps based on time spacing
+        time_coord = cube.coord("time")
+        points = time_coord.units.num2date(time_coord.points)
+        if len(points) < 2:
+            raise ValueError("Not enough time points for rolling window.")
+        # Time step in hours:
+        dt_hours = (points[1] - points[0]).total_seconds() / 3600.0
+        if dt_hours <= 0:
+            raise ValueError("Time coordinate must be increasing.")
+        window_len = int(round(window_hours / dt_hours))
+        if window_len < 1:
+            raise ValueError(
+                f"Window of {window_hours} hours is too small for timestep {dt_hours}."
+            )
+
         window_cube = cube.rolling_window(
-            "time", getattr(iris.analysis, method), window
+            "time", getattr(iris.analysis, method), window_len
         )
         new_cubelist.append(window_cube)
 

src/CSET/operators/collapse.py

@@ -79,12 +79,30 @@ def collapse(
         logging.debug(
             "Extracting common time points as multiple model inputs detected."
         )
-        for cube in cubes:
-            cube.coord("forecast_reference_time").bounds = None
-            cube.coord("forecast_period").bounds = None
-        cubes = cubes.extract_overlapping(
-            ["forecast_reference_time", "forecast_period"]
+        is_power_spectrum = any(
+            cubes[0].coords(coord)
+            for coord in ["frequency", "physical_wavenumber", "wavelength"]
         )
+        if is_power_spectrum:
+            for cube in cubes:
+                cube.coord("time").bounds = None
+            cubes = cubes.extract_overlapping(["time"])
+
+            for cube in cubes:
+                t = cube.coord("time")
+                t.points = t.points.astype(np.float64)
+
+                if t.bounds is not None:
+                    t.bounds = t.bounds.astype(np.float64)
+
+        else:
+            for cube in cubes:
+                cube.coord("forecast_reference_time").bounds = None
+                cube.coord("forecast_period").bounds = None
+            cubes = cubes.extract_overlapping(
+                ["forecast_reference_time", "forecast_period"]
+            )
+
         if len(cubes) == 0:
             raise ValueError("No overlapping times detected in input cubes.")
 
@@ -96,6 +114,7 @@ def collapse(
         warnings.filterwarnings(
             "ignore", "Collapsing spatial coordinate.+without weighting", UserWarning
         )
+
         for cube in iter_maybe(cubes):
             # Apply a mask to check for invalid data, this will allow NaNs to
             # be ignored.
@@ -116,6 +135,7 @@ def collapse(
                 collapsed_cubes.append(
                     cube.collapsed(coordinate, getattr(iris.analysis, method))
                 )
+
     if len(collapsed_cubes) == 1:
         return collapsed_cubes[0]
     else:

src/CSET/operators/plot.py

@@ -2005,6 +2005,7 @@ def plot_line_series(
     # Ensure we have a name for the plot file.
     recipe_title = get_recipe_metadata().get("title", "Untitled")
 
+    print("RECIPE TITLE ", recipe_title)
     num_models = get_num_models(cube)
 
     validate_cube_shape(cube, num_models)
@@ -2120,6 +2121,8 @@ def plot_line_series(
                 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])}]"
 
+            print("TITLE ", title)
+
             # Do the actual plotting.
             plotting_func(
                 cube_slice,

src/CSET/recipes/level_fields/generic_mlevel_power_spectrum_series_mean_all.yaml

@@ -0,0 +1,59 @@
+category: Power Spectrum
+title: Power Spectrum $VARNAME $LEVELTYPE $LEVEL $SUBAREA_NAME
+description: |
+  Extracts and plots the mean power spectral density of $VARNAME from a
+  file at $LEVELTYPE level $LEVEL of $MODEL_NAME.
+
+  Power Spectra calculated on regional domains using Discrete Cosine Transform (DCT)
+  as described in
+  [Denis_etal_2002](https://doi.org/10.1175/1520-0493(2002)130<1812:SDOTDA>2.0.CO;2).
+  The power spectra is then converted into a power spectral density to enable fair comparison
+  between different model resolutions, comparison with observations with different number of data points
+  or comparison of models between different domain sizes, cases, run lengths.
+  In case of ensemble data choose from postage stamp plot or
+  single plot via the single_plot option in the recipe directly.
+
+  Power spectra can be used to identify how the variance (energy) of a field is distributed across scales.
+  The power (energy) is plotted against wavenumber, where low wavenumbers represent large-scale features
+  and high numbers represent small-scale features.  Using a log-log scale for plotting allows for ease of
+  interpretation and the spectra are normalised to allow for comparison between different models.
+
+  The dominant scales of motion can be identified by peaks in the spectrum and the slope of the curve shows
+  how the energy cascades across scales; the slope can be used to identify, for example, the inertial subrange
+  (slope of gradient -5/3).  Caution is advised when interpreting the spectra at high wavenumbers due to,
+  for example, aliasing which can occur when the grid is too coarse to resolve the smallest-scale features
+  and often manifests as spurious spikes or increases in the power at the high wavenumbers.
+
+steps:
+  - operator: read.read_cubes
+    model_names: $MODEL_NAME
+    file_paths: $INPUT_PATHS
+    constraint:
+      operator: constraints.combine_constraints
+      variable_constraint:
+        operator: constraints.generate_var_constraint
+        varname: $VARNAME
+      cell_methods_constraint:
+        operator: constraints.generate_cell_methods_constraint
+        cell_methods: []
+        varname: $VARNAME
+      level_constraint:
+        operator: constraints.generate_level_constraint
+        coordinate: $LEVELTYPE
+        levels: $LEVEL
+    subarea_type: $SUBAREA_TYPE
+    subarea_extent: $SUBAREA_EXTENT
+
+  - operator: power_spectrum.calculate_power_spectrum
+
+  - operator: collapse.collapse
+    method: MEAN
+    coordinate: "time"
+
+  - operator: plot.plot_line_series
+    series_coordinate: "physical_wavenumber"
+    stamp_coordinate: "realization"
+    single_plot: $SPECTRUM_MLEVEL_FIELD_SEQUENCE
+
+  - operator: write.write_cube_to_nc
+    overwrite: True