From e134137690fb18c6aada77584e4e2e78a72f5052 Mon Sep 17 00:00:00 2001 From: Carol Halliwell Date: Mon, 29 Jun 2026 14:35:49 +0100 Subject: [PATCH 1/5] Add functionality to calculate and plot mean and rolling mean of spectra --- .../meta/diagnostics/rose-meta.conf | 58 ++++++++++++++ src/CSET/loaders/power_spectrum.py | 79 +++++++++++++++++++ src/CSET/operators/collapse.py | 35 +++++++- ...mlevel_power_spectrum_series_mean_all.yaml | 59 ++++++++++++++ ...el_power_spectrum_series_mean_rolling.yaml | 60 ++++++++++++++ ...plevel_power_spectrum_series_mean_all.yaml | 59 ++++++++++++++ ...el_power_spectrum_series_mean_rolling.yaml | 60 ++++++++++++++ ...urface_power_spectrum_series_mean_all.yaml | 59 ++++++++++++++ ...ce_power_spectrum_series_mean_rolling.yaml | 60 ++++++++++++++ 9 files changed, 527 insertions(+), 2 deletions(-) create mode 100644 src/CSET/recipes/level_fields/generic_mlevel_power_spectrum_series_mean_all.yaml create mode 100644 src/CSET/recipes/level_fields/generic_mlevel_power_spectrum_series_mean_rolling.yaml create mode 100644 src/CSET/recipes/level_fields/generic_plevel_power_spectrum_series_mean_all.yaml create mode 100644 src/CSET/recipes/level_fields/generic_plevel_power_spectrum_series_mean_rolling.yaml create mode 100644 src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_all.yaml create mode 100644 src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_rolling.yaml diff --git a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf index 77483988b1..0e277dae45 100644 --- a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf +++ b/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=0surface8b + +[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=0surface8c + ####################################################################### # 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=1pressure9b + +[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=1pressure9c + + + [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=2modellevel9b + +[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=2modellevel9c + + + [template variables=SPATIAL_STRUCTURAL_SIMILARITY_MLEVEL] ns=Diagnostics/ModelLevel description=Create spatially mapped structural similarity plots for diff --git a/src/CSET/loaders/power_spectrum.py b/src/CSET/loaders/power_spectrum.py index faf4ad7e4c..2a175e0089 100644 --- a/src/CSET/loaders/power_spectrum.py +++ b/src/CSET/loaders/power_spectrum.py @@ -99,3 +99,82 @@ 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"] + + # Surface (2D) fields. + for atype, field in itertools.product(AGGREGATION_TYPES, conf.SURFACE_FIELDS): + if conf.SPECTRUM_SURFACE_FIELD_AGGREGATION[AGGREGATION_TYPES.index(atype)]: + yield RawRecipe( + recipe=f"generic_surface_power_spectrum_series_mean_{atype}.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, + "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], + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, + } + + # 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], + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, + } + + # 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, + ) diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index dcd2994844..ac18adc4eb 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -36,6 +36,7 @@ def collapse( coordinate: str | list[str], method: str, additional_percent: float = None, + window_hours: float = None, **kwargs, ) -> iris.cube.Cube | iris.cube.CubeList: """Collapse coordinate(s) of a single cube or of every cube in a cube list. @@ -54,10 +55,14 @@ def collapse( be given. method: str Type of collapse i.e. method: 'MEAN', 'MAX', 'MIN', 'MEDIAN', - 'PERCENTILE' getattr creates iris.analysis.MEAN, etc. For PERCENTILE YAML - file requires i.e. method: 'PERCENTILE' additional_percent: 90. + 'PERCENTILE', 'ROLLING_MEAN' getattr creates iris.analysis.MEAN, etc. + For PERCENTILE YAML file requires i.e. method: 'PERCENTILE' additional_percent: 90. + For ROLLING_MEAN YAML file requires method: 'ROLLING_MEAN' and window_hours specified + or set in gui through $WINDOW_LEN_SURFACE additional_percent: float, optional Required for the PERCENTILE method. This is a number between 0 and 100. + window_hours: float + number of hours used for rolling mean (ROLLING_MEAN) Returns ------- @@ -112,6 +117,32 @@ def collapse( cube_max = cube.collapsed(coordinate, iris.analysis.MAX) cube_min = cube.collapsed(coordinate, iris.analysis.MIN) collapsed_cubes.append(cube_max - cube_min) + elif method == "ROLLING_MEAN": + # Rolling mean over window length based on specified window_hours + if window_hours is None: + raise ValueError("ROLLING_MEAN requires kwarg: window_hours=") + # 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}." + ) + # Create rolling window along time + rolled = cube.rolling_window_time_aggregation( + coord="time", + window=window_len, + aggregator=iris.analysis.MEAN, + ) + # Collapse over the window dimension + collapsed_cubes.append(rolled) else: collapsed_cubes.append( cube.collapsed(coordinate, getattr(iris.analysis, method)) diff --git a/src/CSET/recipes/level_fields/generic_mlevel_power_spectrum_series_mean_all.yaml b/src/CSET/recipes/level_fields/generic_mlevel_power_spectrum_series_mean_all.yaml new file mode 100644 index 0000000000..451125ad54 --- /dev/null +++ b/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 diff --git a/src/CSET/recipes/level_fields/generic_mlevel_power_spectrum_series_mean_rolling.yaml b/src/CSET/recipes/level_fields/generic_mlevel_power_spectrum_series_mean_rolling.yaml new file mode 100644 index 0000000000..54dc632c51 --- /dev/null +++ b/src/CSET/recipes/level_fields/generic_mlevel_power_spectrum_series_mean_rolling.yaml @@ -0,0 +1,60 @@ +category: Power Spectrum +title: Power Spectrum $VARNAME $LEVELTYPE $LEVEL $SUBAREA_NAME +description: | + Extracts and plots the rolling mean power spectral density of $VARNAME from a + file at $LEVELTYPE level $LEVEL of $MODEL_NAME using specified window length. + + 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: ROLLING_MEAN + coordinate: "time" + window_hours: $WINDOW_LEN_MLEVEL + + - 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 diff --git a/src/CSET/recipes/level_fields/generic_plevel_power_spectrum_series_mean_all.yaml b/src/CSET/recipes/level_fields/generic_plevel_power_spectrum_series_mean_all.yaml new file mode 100644 index 0000000000..13359f1f17 --- /dev/null +++ b/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 diff --git a/src/CSET/recipes/level_fields/generic_plevel_power_spectrum_series_mean_rolling.yaml b/src/CSET/recipes/level_fields/generic_plevel_power_spectrum_series_mean_rolling.yaml new file mode 100644 index 0000000000..9986ea78d7 --- /dev/null +++ b/src/CSET/recipes/level_fields/generic_plevel_power_spectrum_series_mean_rolling.yaml @@ -0,0 +1,60 @@ +category: Power Spectrum +title: Power Spectrum $VARNAME $LEVELTYPE $LEVEL +description: | + Extracts and plots the rolling mean power spectral density of $VARNAME from a + file at $LEVELTYPE level $LEVEL of $MODEL_NAME using specified window length. + + 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: ROLLING_MEAN + coordinate: "time" + window_hours: $WINDOW_LEN_PLEVEL + + - 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 diff --git a/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_all.yaml b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_all.yaml new file mode 100644 index 0000000000..40e9f53a60 --- /dev/null +++ b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_all.yaml @@ -0,0 +1,59 @@ +category: Power Spectrum +title: "Power Spectrum $VARNAME $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" + stamp_coordinate: "realization" + single_plot: $SPECTRUM_SURFACE_FIELD_SEQUENCE + + - operator: write.write_cube_to_nc + overwrite: True diff --git a/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_rolling.yaml b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_rolling.yaml new file mode 100644 index 0000000000..fa4010cb83 --- /dev/null +++ b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_rolling.yaml @@ -0,0 +1,60 @@ +category: Power Spectrum +title: "Power Spectrum $VARNAME $SUBAREA_NAME" +description: | + Extracts and plots the rolling mean power spectral density of surface $VARNAME over specified window length + + 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: ROLLING_MEAN + coordinate: "time" + window_hours: $WINDOW_LEN_SURFACE + + - operator: plot.plot_line_series + series_coordinate: "physical_wavenumber" + stamp_coordinate: "realization" + single_plot: $SPECTRUM_SURFACE_FIELD_SEQUENCE + + - operator: write.write_cube_to_nc + overwrite: True From fea02bc2158029332990152ca818f8c7e694111d Mon Sep 17 00:00:00 2001 From: Carol Halliwell Date: Mon, 29 Jun 2026 15:12:59 +0100 Subject: [PATCH 2/5] bug fixes in collapse operator when calling rolling mean --- src/CSET/operators/collapse.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index ac18adc4eb..b83132f6d7 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -28,7 +28,10 @@ import numpy as np from CSET._common import iter_maybe -from CSET.operators.aggregate import add_hour_coordinate +from CSET.operators.aggregate import ( + add_hour_coordinate, + rolling_window_time_aggregation, +) def collapse( @@ -136,10 +139,8 @@ def collapse( f"Window of {window_hours} hours is too small for timestep {dt_hours}." ) # Create rolling window along time - rolled = cube.rolling_window_time_aggregation( - coord="time", - window=window_len, - aggregator=iris.analysis.MEAN, + rolled = rolling_window_time_aggregation( + cube, "MEAN", window=window_len ) # Collapse over the window dimension collapsed_cubes.append(rolled) From b32d2e1e34ffe9ca1d16337dfdde83ddff3c546e Mon Sep 17 00:00:00 2001 From: Carol Halliwell Date: Thu, 2 Jul 2026 14:10:02 +0100 Subject: [PATCH 3/5] Bug fix for multi-model plots of mean spectra --- .../meta/diagnostics/rose-meta.conf | 12 ++-- src/CSET/loaders/power_spectrum.py | 15 ++++ src/CSET/operators/collapse.py | 68 ++++++++----------- ...urface_power_spectrum_series_mean_all.yaml | 1 - ...ce_power_spectrum_series_mean_rolling.yaml | 6 +- 5 files changed, 51 insertions(+), 51 deletions(-) diff --git a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf index 0e277dae45..a77f645df5 100644 --- a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf @@ -237,7 +237,7 @@ description=Create case aggregated power spectrum plots for surface. Option2: Rolling time mean. type=python_boolean,python_boolean compulsory=true -sort-key=0surface8b +sort-key=c1surface8b [template variables=WINDOW_LEN_SURFACE] ns=Diagnostics/Fields @@ -245,7 +245,7 @@ description= The length of the window (in hours) that the rolling mean uses. help= A single value type=integer compulsory=true -sort-key=0surface8c +sort-key=c1surface8c ####################################################################### @@ -489,7 +489,7 @@ description=Create case aggregated power spectrum plots for specified pressure l Option2: Rolling time mean. type=python_boolean,python_boolean compulsory=true -sort-key=1pressure9b +sort-key=c2pressure9b [template variables=WINDOW_LEN_PLEVEL] ns=Diagnostics/Pressure @@ -497,7 +497,7 @@ description= The length of the window (in hours) that the rolling mean uses. help= A single value type=integer compulsory=true -sort-key=1pressure9c +sort-key=c2pressure9c @@ -729,7 +729,7 @@ description=Create case aggregated power spectrum plots for specified model leve Option2: Rolling time mean. type=python_boolean,python_boolean compulsory=true -sort-key=2modellevel9b +sort-key=c3modellevel9b [template variables=WINDOW_LEN_MLEVEL] ns=Diagnostics/ModelLevel @@ -737,7 +737,7 @@ description= The length of the window (in hours) that the rolling mean uses. help= A single value type=integer compulsory=true -sort-key=2modellevel9c +sort-key=c3modellevel9c diff --git a/src/CSET/loaders/power_spectrum.py b/src/CSET/loaders/power_spectrum.py index 2a175e0089..7d161f5c1c 100644 --- a/src/CSET/loaders/power_spectrum.py +++ b/src/CSET/loaders/power_spectrum.py @@ -111,10 +111,15 @@ def load(conf: Config): 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, @@ -134,8 +139,13 @@ def load(conf: Config): "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 @@ -162,8 +172,13 @@ def load(conf: Config): "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 diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index b83132f6d7..ee2ef7112f 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -28,10 +28,7 @@ import numpy as np from CSET._common import iter_maybe -from CSET.operators.aggregate import ( - add_hour_coordinate, - rolling_window_time_aggregation, -) +from CSET.operators.aggregate import add_hour_coordinate def collapse( @@ -39,7 +36,6 @@ def collapse( coordinate: str | list[str], method: str, additional_percent: float = None, - window_hours: float = None, **kwargs, ) -> iris.cube.Cube | iris.cube.CubeList: """Collapse coordinate(s) of a single cube or of every cube in a cube list. @@ -58,14 +54,10 @@ def collapse( be given. method: str Type of collapse i.e. method: 'MEAN', 'MAX', 'MIN', 'MEDIAN', - 'PERCENTILE', 'ROLLING_MEAN' getattr creates iris.analysis.MEAN, etc. - For PERCENTILE YAML file requires i.e. method: 'PERCENTILE' additional_percent: 90. - For ROLLING_MEAN YAML file requires method: 'ROLLING_MEAN' and window_hours specified - or set in gui through $WINDOW_LEN_SURFACE + 'PERCENTILE' getattr creates iris.analysis.MEAN, etc. For PERCENTILE YAML + file requires i.e. method: 'PERCENTILE' additional_percent: 90. additional_percent: float, optional Required for the PERCENTILE method. This is a number between 0 and 100. - window_hours: float - number of hours used for rolling mean (ROLLING_MEAN) Returns ------- @@ -87,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.") @@ -104,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. @@ -120,34 +131,11 @@ def collapse( cube_max = cube.collapsed(coordinate, iris.analysis.MAX) cube_min = cube.collapsed(coordinate, iris.analysis.MIN) collapsed_cubes.append(cube_max - cube_min) - elif method == "ROLLING_MEAN": - # Rolling mean over window length based on specified window_hours - if window_hours is None: - raise ValueError("ROLLING_MEAN requires kwarg: window_hours=") - # 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}." - ) - # Create rolling window along time - rolled = rolling_window_time_aggregation( - cube, "MEAN", window=window_len - ) - # Collapse over the window dimension - collapsed_cubes.append(rolled) else: collapsed_cubes.append( cube.collapsed(coordinate, getattr(iris.analysis, method)) ) + if len(collapsed_cubes) == 1: return collapsed_cubes[0] else: diff --git a/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_all.yaml b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_all.yaml index 40e9f53a60..8fb46b553d 100644 --- a/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_all.yaml +++ b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_all.yaml @@ -52,7 +52,6 @@ steps: - operator: plot.plot_line_series series_coordinate: "physical_wavenumber" - stamp_coordinate: "realization" single_plot: $SPECTRUM_SURFACE_FIELD_SEQUENCE - operator: write.write_cube_to_nc diff --git a/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_rolling.yaml b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_rolling.yaml index fa4010cb83..2ff7dac596 100644 --- a/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_rolling.yaml +++ b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_rolling.yaml @@ -46,14 +46,12 @@ steps: - operator: power_spectrum.calculate_power_spectrum - - operator: collapse.collapse + - operator: aggregate.rolling_window_time_aggregation method: ROLLING_MEAN - coordinate: "time" - window_hours: $WINDOW_LEN_SURFACE + window: $WINDOW_LEN_SURFACE - operator: plot.plot_line_series series_coordinate: "physical_wavenumber" - stamp_coordinate: "realization" single_plot: $SPECTRUM_SURFACE_FIELD_SEQUENCE - operator: write.write_cube_to_nc From 6eaf101d0a343a9e811fb236d8eca7cb69876ba0 Mon Sep 17 00:00:00 2001 From: Carol Halliwell Date: Thu, 2 Jul 2026 14:42:19 +0100 Subject: [PATCH 4/5] changes to rolling_window_time_aggregation to change window hour into number of points --- src/CSET/operators/aggregate.py | 24 +++++++++++++++++-- ...ce_power_spectrum_series_mean_rolling.yaml | 4 ++-- 2 files changed, 24 insertions(+), 4 deletions(-) diff --git a/src/CSET/operators/aggregate.py b/src/CSET/operators/aggregate.py index 0fa3f2e725..f3c0f9fe0a 100644 --- a/src/CSET/operators/aggregate.py +++ b/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=") + # 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) diff --git a/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_rolling.yaml b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_rolling.yaml index 2ff7dac596..76a3306e06 100644 --- a/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_rolling.yaml +++ b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_rolling.yaml @@ -47,8 +47,8 @@ steps: - operator: power_spectrum.calculate_power_spectrum - operator: aggregate.rolling_window_time_aggregation - method: ROLLING_MEAN - window: $WINDOW_LEN_SURFACE + method: MEAN + window_hours: $WINDOW_LEN_SURFACE - operator: plot.plot_line_series series_coordinate: "physical_wavenumber" From 7a8cc3587361478d3f16e7dcdf177ea8fbaf127d Mon Sep 17 00:00:00 2001 From: Carol Halliwell Date: Tue, 7 Jul 2026 16:01:19 +0100 Subject: [PATCH 5/5] Changes to power_spectrum loader --- src/CSET/loaders/power_spectrum.py | 13 +++++++++---- src/CSET/operators/plot.py | 3 +++ ...eric_surface_power_spectrum_series_mean_all.yaml | 2 +- 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/src/CSET/loaders/power_spectrum.py b/src/CSET/loaders/power_spectrum.py index 7d161f5c1c..5afbfb13d5 100644 --- a/src/CSET/loaders/power_spectrum.py +++ b/src/CSET/loaders/power_spectrum.py @@ -102,18 +102,23 @@ def load(conf: Config): # 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)]: + # 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", + # "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 diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index 051e438bd0..6cb3db4d83 100644 --- a/src/CSET/operators/plot.py +++ b/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, diff --git a/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_all.yaml b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_all.yaml index 8fb46b553d..a22fc8033a 100644 --- a/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_all.yaml +++ b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_all.yaml @@ -1,5 +1,5 @@ category: Power Spectrum -title: "Power Spectrum $VARNAME $SUBAREA_NAME" +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