diff --git a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf index 77483988b..a77f645df 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=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 diff --git a/src/CSET/loaders/power_spectrum.py b/src/CSET/loaders/power_spectrum.py index faf4ad7e4..5afbfb13d 100644 --- a/src/CSET/loaders/power_spectrum.py +++ b/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, + ) diff --git a/src/CSET/operators/aggregate.py b/src/CSET/operators/aggregate.py index 0fa3f2e72..f3c0f9fe0 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/operators/collapse.py b/src/CSET/operators/collapse.py index dcd299484..ee2ef7112 100644 --- a/src/CSET/operators/collapse.py +++ b/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: diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index 051e438bd..6cb3db4d8 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/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 000000000..451125ad5 --- /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 000000000..54dc632c5 --- /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 000000000..13359f1f1 --- /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 000000000..9986ea78d --- /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 000000000..a22fc8033 --- /dev/null +++ b/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 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 000000000..76a3306e0 --- /dev/null +++ b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series_mean_rolling.yaml @@ -0,0 +1,58 @@ +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: aggregate.rolling_window_time_aggregation + method: MEAN + window_hours: $WINDOW_LEN_SURFACE + + - operator: plot.plot_line_series + series_coordinate: "physical_wavenumber" + single_plot: $SPECTRUM_SURFACE_FIELD_SEQUENCE + + - operator: write.write_cube_to_nc + overwrite: True