Regional power spectra using Discrete Cosine Transform (DCT)

pyproject.toml

@@ -30,6 +30,7 @@ dependencies = [
     "markdown-it-py >= 3.0",
     "nc-time-axis",
     "iris-grib",
+    "scipy",
     "scikit-image",
 ]
 

requirements/environment.yml

@@ -13,6 +13,7 @@ dependencies:
   - markdown-it-py >= 3.0
   - nc-time-axis
   - iris-grib
+  - scipy
   - importlib_resources # For python < 3.12
   - scikit-image # For image processing techniques.
 

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

@@ -212,6 +212,24 @@ type=python_boolean
 compulsory=true
 sort-key=0surface7b
 
+[template variables=SPECTRUM_SURFACE_FIELD]
+ns=Diagnostics/Fields
+description=Create power spectra plots of specified surface fields.
+            Use SPECTRUM_SURFACE_FIELD_SEQUENCE to set plotting mode.
+type=python_boolean
+compulsory=true
+sort-key=0surface8
+
+[template variables=SPECTRUM_SURFACE_FIELD_SEQUENCE]
+ns=Diagnostics/Fields
+description=Select analysis method for output surface field spectra.
+            Set to True for spectrum at each model output time.
+            Set to False for a spectrum of each ensemble member.
+type=python_boolean
+compulsory=true
+sort-key=0surface8a
+
+
 #######################################################################
 # Pressure level fields.
 [Diagnostics/Pressure]
@@ -393,6 +411,23 @@ type=real,real
 compulsory=true
 sort-key=1pressure8b
 
+[template variables=SPECTRUM_PLEVEL_FIELD]
+ns=Diagnostics/Pressure
+description=Create spectrum of specified pressure level fields.
+            Use SPECTRUM_PLEVEL_FIELD_SEQUENCE to set plotting mode.
+type=python_boolean
+compulsory=true
+sort-key=1pressure9
+
+[template variables=SPECTRUM_PLEVEL_FIELD_SEQUENCE]
+ns=Diagnostics/Pressure
+description=Select analysis method for output pressure level spectra.
+            Set to True for spectrum at each model output time.
+            Set to False for a spectrum of each ensemble member.
+type=python_boolean
+compulsory=true
+sort-key=1pressure9a
+
 [template variables=SPATIAL_STRUCTURAL_SIMILARITY_PLEVEL_FIELD]
 ns=Diagnostics/Pressure
 description=Create spatially mapped structural similarity plots for
@@ -402,7 +437,7 @@ description=Create spatially mapped structural similarity plots for
 help=Requires a field with a time coordinate of "time" or "hour".
 type=python_boolean
 compulsory=true
-sort-key=1pressure9
+sort-key=1pressure9e
 
 [template variables=MEAN_STRUCTURAL_SIMILARITY_PLEVEL_FIELD]
 ns=Diagnostics/Pressure
@@ -412,7 +447,8 @@ description=Create time series plots of the mean structural similarity
 help=Requires a field with a time coordinate of "time" or "hour".
 type=python_boolean
 compulsory=true
-sort-key=1pressure9b
+sort-key=1pressure9f
+
 
 ########################################################################
 # Model-level fields.
@@ -595,6 +631,23 @@ type=real,real
 compulsory=true
 sort-key=2modellevel8c
 
+[template variables=SPECTRUM_MLEVEL_FIELD]
+ns=Diagnostics/ModelLevel
+description=Create spectrum of specified model level fields.
+            Use SPECTRUM_MLEVEL_FIELD_SEQUENCE to set plotting mode.
+type=python_boolean
+compulsory=true
+sort-key=2modellevel9
+
+[template variables=SPECTRUM_MLEVEL_FIELD_SEQUENCE]
+ns=Diagnostics/ModelLevel
+description=Select analysis method for output model level spectrum.
+            Set to True for spectrum at each model output time.
+            Set to False for a spectrum of each ensemble member.
+type=python_boolean
+compulsory=true
+sort-key=2modellevel9a
+
 [template variables=SPATIAL_STRUCTURAL_SIMILARITY_MLEVEL]
 ns=Diagnostics/ModelLevel
 description=Create spatially mapped structural similarity plots for
@@ -603,7 +656,7 @@ description=Create spatially mapped structural similarity plots for
 help=Requires a field with a time coordinate of "time" or "hour".
 type=python_boolean
 compulsory=true
-sort-key=2modellevel9
+sort-key=2modellevel9e
 
 [template variables=MEAN_STRUCTURAL_SIMILARITY_MLEVEL]
 ns=Diagnostics/ModelLevel
@@ -613,7 +666,8 @@ description=Create time series plots of the mean structural similarity
 help=Requires a field with a time coordinate of "time" or "hour".
 type=python_boolean
 compulsory=true
-sort-key=2modellevel9b
+sort-key=2modellevel9f
+
 
 #######################################################
 # Verification

src/CSET/loaders/__init__.py

@@ -22,6 +22,7 @@
 from CSET.loaders import (
     aoa,
     histograms,
+    power_spectrum,
     profiles,
     spatial_difference_field,
     spatial_field,
@@ -33,6 +34,7 @@
 __all__ = [
     "aoa",
     "histograms",
+    "power_spectrum",
     "profiles",
     "spatial_difference_field",
     "spatial_field",

src/CSET/loaders/power_spectrum.py

@@ -0,0 +1,95 @@
+# © Crown copyright, Met Office (2022-2025) and CSET contributors.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Load power spectrum recipes."""
+
+import itertools
+
+from CSET.recipes import Config, RawRecipe, get_models
+
+
+def load(conf: Config):
+    """Yield recipes from the given workflow configuration."""
+    # Load a list of model detail dictionaries.
+    models = get_models(conf.asdict())
+
+    # Surface (2D) fields.
+    if conf.SPECTRUM_SURFACE_FIELD:
+        for field in conf.SURFACE_FIELDS:
+            yield RawRecipe(
+                recipe="generic_surface_power_spectrum_series.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,
+                },
+                model_ids=[model["id"] for model in models],
+                aggregation=False,
+            )
+
+    # Pressure level fields.
+    if conf.SPECTRUM_PLEVEL_FIELD:
+        for field, plevel in itertools.product(
+            conf.PRESSURE_LEVEL_FIELDS,
+            conf.PRESSURE_LEVELS,
+        ):
+            yield RawRecipe(
+                recipe="generic_level_power_spectrum_series.yaml",
+                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,
+                },
+                model_ids=[model["id"] for model in models],
+                aggregation=False,
+            )
+
+    # Model level fields
+    if conf.SPECTRUM_MLEVEL_FIELD:
+        for field, mlevel in itertools.product(
+            conf.MODEL_LEVEL_FIELDS,
+            conf.MODEL_LEVELS,
+        ):
+            yield RawRecipe(
+                recipe="generic_level_power_spectrum_series.yaml",
+                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,
+                },
+                model_ids=[model["id"] for model in models],
+                aggregation=False,
+            )