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

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

@@ -229,6 +229,15 @@ 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.
+type=python_boolean
+compulsory=true
+sort-key=c1surface8b
+
 
 #######################################################################
 # Pressure level fields.
@@ -463,6 +472,17 @@ 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.
+type=python_boolean
+compulsory=true
+sort-key=c2pressure9b
+
+
+
 [template variables=SPATIAL_STRUCTURAL_SIMILARITY_PLEVEL_FIELD]
 ns=Diagnostics/Pressure
 description=Create spatially mapped structural similarity plots for
@@ -683,6 +703,16 @@ 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.
+type=python_boolean
+compulsory=true
+sort-key=c3modellevel9b
+
+
 [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,86 @@ def load(conf: Config):
                 model_ids=[model["id"] for model in models],
                 aggregation=False,
             )
+
+    # Create a list of case aggregation types.
+    AGGREGATION_TYPES = ["all"]
+    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)
+        print("ALL INFO ", atype, index, AGGREGATION_TYPES)
+        if len(AGGREGATION_TYPES) > index and AGGREGATION_TYPES[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,
+                },
+                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,
+            }
+
+            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,
+            }
+
+            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/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

@@ -2050,6 +2050,7 @@ def plot_line_series(
     # Ensure we have a name for the plot file.
     recipe_title = get_recipe_metadata().get("title", iter_maybe(cube)[0].name())
 
+    print("RECIPE TITLE ", recipe_title)
     num_models = get_num_models(cube)
 
     validate_cube_shape(cube, num_models)
@@ -2165,6 +2166,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

src/CSET/recipes/level_fields/generic_plevel_power_spectrum_series_mean_all.yaml

@@ -0,0 +1,59 @@
+category: Power Spectrum
+title: Power Spectrum $VARNAME $LEVELTYPE $LEVEL
+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_PLEVEL_FIELD_SEQUENCE
+
+  - operator: write.write_cube_to_nc
+    overwrite: True

src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_all.yaml

@@ -0,0 +1,58 @@
+category: Power Spectrum
+title: "Power Spectrum $VARNAME\n Aggregation over all cases $SUBAREA_NAME"
+description: |
+  Extracts and plots the mean power spectral density of surface $VARNAME over whole run
+
+  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.
+
+  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.
+  Alternatively the power spectral density can be plotted against wavelength (needs to be selected in recipe),
+  where low wavelengths represent small-scale features and high wavelengths represent large-scale features.
+
+
+  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
+    file_paths: $INPUT_PATHS
+    model_names: $MODEL_NAME
+    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
+      pressure_level_constraint:
+        operator: constraints.generate_level_constraint
+        coordinate: pressure
+        levels: []
+    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"
+    single_plot: $SPECTRUM_SURFACE_FIELD_SEQUENCE
+
+  - operator: write.write_cube_to_nc
+    overwrite: True