Calculate mean spectra over the case for surface fields

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 in case.
+type=python_boolean
+compulsory=true
+sort-key=c1surface8b
+
 
 #######################################################################
 # Pressure level fields.
@@ -463,6 +472,15 @@ 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 in case.
+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 +701,15 @@ 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 in case.
+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/cset_workflow/rose-suite.conf.example

@@ -219,10 +219,13 @@ SPATIAL_SURFACE_FIELD_AGGREGATION=False,False,False,False
 SPATIAL_SURFACE_FIELD_METHOD=[""]
 SPECTRUM_MLEVEL_FIELD=False
 SPECTRUM_MLEVEL_FIELD_SEQUENCE=False
+SPECTRUM_MLEVEL_FIELD_AGGREGATION=False
 SPECTRUM_PLEVEL_FIELD=False
 SPECTRUM_PLEVEL_FIELD_SEQUENCE=False
+SPECTRUM_PLEVEL_FIELD_AGGREGATION=False
 SPECTRUM_SURFACE_FIELD=False
 SPECTRUM_SURFACE_FIELD_SEQUENCE=False
+SPECTRUM_SURFACE_FIELD_AGGREGATION=False
 !!SUBAREA_EXTENT=0,0,0,0
 !!SUBAREA_NAME=""
 !!SUBAREA_TYPE="realworld"

src/CSET/loaders/power_spectrum.py

@@ -52,6 +52,7 @@ def load(conf: Config):
             conf.PRESSURE_LEVEL_FIELDS,
             conf.PRESSURE_LEVELS,
         ):
+            print("PRESSURE FIELDS ", field, plevel)
             yield RawRecipe(
                 recipe="generic_plevel_power_spectrum_series.yaml",
                 variables={
@@ -99,3 +100,74 @@ def load(conf: Config):
                 model_ids=[model["id"] for model in models],
                 aggregation=False,
             )
+
+    # Surface (2D) fields.
+
+    #    if conf.SPECTRUM_SURFACE_FIELD_AGGREGATION:
+    if conf.SPECTRUM_SURFACE_FIELD and conf.SPECTRUM_SURFACE_FIELD_AGGREGATION:
+        for field in conf.SURFACE_FIELDS:
+            print("LOADER AGGREGATION ")
+            yield RawRecipe(
+                recipe="generic_surface_power_spectrum_series_mean_case.yaml",
+                variables={
+                    "VARNAME": field,
+                    "MODEL_NAME": [model["name"] for model in models],
+                    "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.
+    if conf.SPECTRUM_PLEVEL_FIELD and conf.SPECTRUM_PLEVEL_FIELD_AGGREGATION:
+        for field in conf.PRESSURE_LEVEL_FIELDS:
+            print("PRESSURE FIELDS AGG", field, plevel)
+            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="generic_plevel_power_spectrum_series_mean_case.yaml",
+                variables=variables,
+                model_ids=[model["id"] for model in models],
+                aggregation=True,
+            )
+
+    # Model level fields.
+    if conf.SPECTRUM_MLEVEL_FIELD and conf.POWER_SPECTRUM_MLEVEL_FIELD_AGGREGATION:
+        for field in conf.MODEL_LEVEL_FIELDS:
+            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="generic_mlevel_power_spectrum_series_mean_case.yaml",
+                variables=variables,
+                model_ids=[model["id"] for model in models],
+                aggregation=True,
+            )

src/CSET/operators/collapse.py

@@ -78,13 +78,30 @@ def collapse(
 
     # Retain only common time points between different models if multiple model inputs.
     if isinstance(cubes, iris.cube.CubeList) and len(cubes) > 1:
-        logger.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.")
 

src/CSET/operators/plot.py

@@ -2249,6 +2249,9 @@ def plot_line_series(
 
     # Iterate over all cubes and extract coordinate to plot.
     cubes = iris.cube.CubeList(iter_maybe(cube))
+
+    print("CUBES in plot_line_series ", cubes)
+
     coords = []
     for model_cube in cubes:
         try:

src/CSET/operators/power_spectrum.py

@@ -310,6 +310,8 @@ def _power_spectrum(cube: iris.cube.Cube) -> iris.cube.Cube:
     time_coord = cube.coord("time")
     time_points = time_coord.units.num2date(time_coord.points)
 
+    print("CUBE in _power_spectrum ", cube)
+
     if cube.ndim == 2:
         cube_3d = cube.data[np.newaxis, :, :]
         logger.debug("Adding in new axis for a 2 dimensional cube.")

src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_case.yaml

@@ -0,0 +1,58 @@
+category: Power Spectrum
+title: "Power Spectrum $VARNAME\n $SUBAREA_NAME Mean over the case"
+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