diff --git a/pyproject.toml b/pyproject.toml index bdfab840f..186aa5a65 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,6 +30,7 @@ dependencies = [ "markdown-it-py >= 3.0", "nc-time-axis", "iris-grib", + "scipy", "scikit-image", ] diff --git a/requirements/environment.yml b/requirements/environment.yml index f49e965f3..501cb590b 100644 --- a/requirements/environment.yml +++ b/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. diff --git a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf index c304a0c3f..9dadb2211 100644 --- a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf +++ b/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 diff --git a/src/CSET/loaders/__init__.py b/src/CSET/loaders/__init__.py index 385ec97ea..172596700 100644 --- a/src/CSET/loaders/__init__.py +++ b/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", diff --git a/src/CSET/loaders/power_spectrum.py b/src/CSET/loaders/power_spectrum.py new file mode 100644 index 000000000..6b3256243 --- /dev/null +++ b/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, + ) diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index a4b08eab9..c8d3ccb27 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -35,6 +35,7 @@ import matplotlib.colors as mcolors import matplotlib.pyplot as plt import numpy as np +import scipy.fft as fft from markdown_it import MarkdownIt from CSET._common import ( @@ -1269,6 +1270,192 @@ def _plot_and_save_postage_stamps_in_single_plot_histogram_series( plt.close(fig) +def _plot_and_save_power_spectrum_series( + cubes: iris.cube.Cube | iris.cube.CubeList, + filename: str, + title: str, + **kwargs, +): + """Plot and save a power spectrum series. + + Parameters + ---------- + cubes: Cube or CubeList + 2 dimensional Cube or CubeList of the data to plot as power spectrum. + filename: str + Filename of the plot to write. + title: str + Plot title. + """ + fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k") + ax = plt.gca() + + model_colors_map = _get_model_colors_map(cubes) + + for cube in iter_maybe(cubes): + # Calculate power spectrum + + # Extract time coordinate and convert to datetime + time_coord = cube.coord("time") + time_points = time_coord.units.num2date(time_coord.points) + + # Choose one time point (e.g., the first one) + target_time = time_points[0] + + # Bind target_time inside the lambda using a default argument + time_constraint = iris.Constraint( + time=lambda cell, target_time=target_time: cell.point == target_time + ) + + cube = cube.extract(time_constraint) + + if cube.ndim == 2: + cube_3d = cube.data[np.newaxis, :, :] + logging.debug("Adding in new axis for a 2 dimensional cube.") + elif cube.ndim == 3: + cube_3d = cube.data + else: + raise ValueError("Cube dimensions unsuitable for power spectra code") + raise ValueError( + f"Cube is {cube.ndim} dimensional. Cube should be 2 or 3 dimensional." + ) + + # Calculate spectra + ps_array = _DCT_ps(cube_3d) + + ps_cube = iris.cube.Cube( + ps_array, + long_name="power_spectra", + ) + + ps_cube.attributes["model_name"] = cube.attributes.get("model_name") + + # Create a frequency/wavelength array for coordinate + ps_len = ps_cube.data.shape[1] + freqs = np.arange(1, ps_len + 1) + freq_coord = iris.coords.DimCoord(freqs, long_name="frequency", units="1") + + # Convert datetime to numeric time using original units + numeric_time = time_coord.units.date2num(time_points) + # Create a new DimCoord with numeric time + new_time_coord = iris.coords.DimCoord( + numeric_time, standard_name="time", units=time_coord.units + ) + + # Add time and frequency coordinate to spectra cube. + ps_cube.add_dim_coord(new_time_coord.copy(), 0) + ps_cube.add_dim_coord(freq_coord.copy(), 1) + + # Extract data from the cube + frequency = ps_cube.coord("frequency").points + power_spectrum = ps_cube.data + + label = None + color = "black" + if model_colors_map: + label = ps_cube.attributes.get("model_name") + color = model_colors_map[label] + ax.plot(frequency, power_spectrum[0], color=color, label=label) + + # Add some labels and tweak the style. + ax.set_title(title, fontsize=16) + ax.set_xlabel("Wavenumber", fontsize=14) + ax.set_ylabel("Power", fontsize=14) + ax.tick_params(axis="both", labelsize=12) + + # Set log-log scale + ax.set_xscale("log") + ax.set_yscale("log") + + # Overlay grid-lines onto power spectrum plot. + ax.grid(linestyle="--", color="grey", linewidth=1) + if model_colors_map: + ax.legend(loc="best", ncol=1, frameon=False, fontsize=16) + + # Save plot. + fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution()) + logging.info("Saved line plot to %s", filename) + plt.close(fig) + + +def _plot_and_save_postage_stamp_power_spectrum_series( + cube: iris.cube.Cube, + filename: str, + title: str, + stamp_coordinate: str, + **kwargs, +): + """Plot and save postage (ensemble members) stamps for a power spectrum series. + + Parameters + ---------- + cube: Cube + 2 dimensional Cube of the data to plot as power spectrum. + filename: str + Filename of the plot to write. + title: str + Plot title. + stamp_coordinate: str + Coordinate that becomes different plots. + """ + # Use the smallest square grid that will fit the members. + grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points)))) + + fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k") + # Make a subplot for each member. + for member, subplot in zip( + cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False + ): + # Implicit interface is much easier here, due to needing to have the + # cartopy GeoAxes generated. + plt.subplot(grid_size, grid_size, subplot) + + frequency = member.coord("frequency").points + power_spectrum = member.data + + ax = plt.gca() + ax.plot(frequency, power_spectrum[0]) + ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}") + + # Overall figure title. + fig.suptitle(title, fontsize=16) + + fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution()) + logging.info("Saved power spectra postage stamp plot to %s", filename) + plt.close(fig) + + +def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series( + cube: iris.cube.Cube, + filename: str, + title: str, + stamp_coordinate: str, + **kwargs, +): + fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k") + ax.set_title(title, fontsize=16) + ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14) + ax.set_ylabel("Power", fontsize=14) + # Loop over all slices along the stamp_coordinate + for member in cube.slices_over(stamp_coordinate): + frequency = member.coord("frequency").points + power_spectrum = member.data + ax.plot( + frequency, + power_spectrum[0], + label=f"Member #{member.coord(stamp_coordinate).points[0]}", + ) + + # Add a legend + ax.legend(fontsize=16) + + # Save the figure to a file + plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution()) + + # Close the figure + plt.close(fig) + + def _spatial_plot( method: Literal["contourf", "pcolormesh"], cube: iris.cube.Cube, @@ -2412,3 +2599,251 @@ def plot_histogram_series( _make_plot_html_page(complete_plot_index) return cubes + + +def plot_power_spectrum_series( + cubes: iris.cube.Cube | iris.cube.CubeList, + filename: str = None, + sequence_coordinate: str = "time", + stamp_coordinate: str = "realization", + single_plot: bool = False, + **kwargs, +) -> iris.cube.Cube | iris.cube.CubeList: + """Plot a power spectrum plot for each vertical level provided. + + A power spectrum plot can be plotted, but if the sequence_coordinate (i.e. time) + is present then a sequence of plots will be produced using the time slider + functionality to scroll through power spectrum against time. If a + stamp_coordinate is present then postage stamp plots will be produced. If + stamp_coordinate and single_plot is True, all postage stamp plots will be + plotted in a single plot instead of separate postage stamp plots. + + Parameters + ---------- + cubes: Cube | iris.cube.CubeList + Iris cube or CubeList of the data to plot. It should have a single dimension other + than the stamp coordinate. + The cubes should cover the same phenomenon i.e. all cubes contain temperature data. + We do not support different data such as temperature and humidity in the same CubeList for plotting. + filename: str, optional + Name of the plot to write, used as a prefix for plot sequences. Defaults + to the recipe name. + sequence_coordinate: str, optional + Coordinate about which to make a plot sequence. Defaults to ``"time"``. + This coordinate must exist in the cube and will be used for the time + slider. + stamp_coordinate: str, optional + Coordinate about which to plot postage stamp plots. Defaults to + ``"realization"``. + single_plot: bool, optional + If True, all postage stamp plots will be plotted in a single plot. If + False, each postage stamp plot will be plotted separately. Is only valid + if stamp_coordinate exists and has more than a single point. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + The original Cube or CubeList (so further operations can be applied). + Plotted data. + + Raises + ------ + ValueError + If the cube doesn't have the right dimensions. + TypeError + If the cube isn't a Cube or CubeList. + """ + recipe_title = get_recipe_metadata().get("title", "Untitled") + + cubes = iter_maybe(cubes) + # Ensure we have a name for the plot file. + if filename is None: + filename = slugify(recipe_title) + + # Internal plotting function. + plotting_func = _plot_and_save_power_spectrum_series + + num_models = _get_num_models(cubes) + + _validate_cube_shape(cubes, num_models) + + # If several power spectra are plotted with time as sequence_coordinate for the + # time slider option. + for cube in cubes: + try: + cube.coord(sequence_coordinate) + except iris.exceptions.CoordinateNotFoundError as err: + raise ValueError( + f"Cube must have a {sequence_coordinate} coordinate." + ) from err + + # Make postage stamp plots if stamp_coordinate exists and has more than a + # single point. If single_plot is True: + # -- all postage stamp plots will be plotted in a single plot instead of + # separate postage stamp plots. + # -- model names (hidden in cube attrs) are ignored, that is stamp plots are + # produced per single model only + if num_models == 1: + if ( + stamp_coordinate in [c.name() for c in cubes[0].coords()] + and cubes[0].coord(stamp_coordinate).shape[0] > 1 + ): + if single_plot: + plotting_func = ( + _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series + ) + else: + plotting_func = _plot_and_save_postage_stamp_power_spectrum_series + cube_iterables = cubes[0].slices_over(sequence_coordinate) + else: + all_points = sorted( + set( + itertools.chain.from_iterable( + cb.coord(sequence_coordinate).points for cb in cubes + ) + ) + ) + all_slices = list( + itertools.chain.from_iterable( + cb.slices_over(sequence_coordinate) for cb in cubes + ) + ) + # Matched slices (matched by seq coord point; it may happen that + # evaluated models do not cover the same seq coord range, hence matching + # necessary) + cube_iterables = [ + iris.cube.CubeList( + s for s in all_slices if s.coord(sequence_coordinate).points[0] == point + ) + for point in all_points + ] + + plot_index = [] + nplot = np.size(cube.coord(sequence_coordinate).points) + # Create a plot for each value of the sequence coordinate. Allowing for + # multiple cubes in a CubeList to be plotted in the same plot for similar + # sequence values. Passing a CubeList into the internal plotting function + # for similar values of the sequence coordinate. cube_slice can be an + # iris.cube.Cube or an iris.cube.CubeList. + for cube_slice in cube_iterables: + single_cube = cube_slice + if isinstance(cube_slice, iris.cube.CubeList): + single_cube = cube_slice[0] + + # Use sequence value so multiple sequences can merge. + sequence_value = single_cube.coord(sequence_coordinate).points[0] + plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png" + coord = single_cube.coord(sequence_coordinate) + # Format the coordinate value in a unit appropriate way. + title = f"{recipe_title}\n [{coord.units.title(coord.points[0])}]" + # Use sequence (e.g. time) bounds if plotting single non-sequence outputs + if nplot == 1 and coord.has_bounds: + if np.size(coord.bounds) > 1: + title = f"{recipe_title}\n [{coord.units.title(coord.bounds[0][0])} to {coord.units.title(coord.bounds[0][1])}]" + # Do the actual plotting. + plotting_func( + cube_slice, + filename=plot_filename, + stamp_coordinate=stamp_coordinate, + title=title, + ) + plot_index.append(plot_filename) + + # Add list of plots to plot metadata. + complete_plot_index = _append_to_plot_index(plot_index) + + # Make a page to display the plots. + _make_plot_html_page(complete_plot_index) + + return cubes + + +def _DCT_ps(y_3d): + """Calculate power spectra for regional domains. + + Parameters + ---------- + y_3d: 3D array + 3 dimensional array to calculate spectrum for. + (2D field data with 3rd dimension of time) + + Returns: ps_array + Array of power spectra values calculated for input field (for each time) + + Method for regional domains: + Calculate power spectra over limited area domain using Discrete Cosine Transform (DCT) + as described in Denis et al 2002 [Denis_etal_2002]_. + + References + ---------- + .. [Denis_etal_2002] Bertrand Denis, Jean Côté and René Laprise (2002) + "Spectral Decomposition of Two-Dimensional Atmospheric Fields on + Limited-Area Domains Using the Discrete Cosine Transform (DCT)" + Monthly Weather Review, Vol. 130, 1812-1828 + doi: https://doi.org/10.1175/1520-0493(2002)130<1812:SDOTDA>2.0.CO;2 + """ + Nt, Ny, Nx = y_3d.shape + + # Max coefficient + Nmin = min(Nx - 1, Ny - 1) + + # Create alpha matrix (of wavenumbers) + alpha_matrix = _create_alpha_matrix(Ny, Nx) + + # Prepare output array + ps_array = np.zeros((Nt, Nmin)) + + # Loop over time to get spectrum for each time. + for t in range(Nt): + y_2d = y_3d[t] + + # Apply 2D DCT to transform y_3d[t] from physical space to spectral space. + # fkk is a 2D array of DCT coefficients, representing the amplitudes of + # cosine basis functions at different spatial frequencies. + + # normalise spectrum to allow comparison between models. + fkk = fft.dctn(y_2d, norm="ortho") + + # Normalise fkk + fkk = fkk / np.sqrt(Ny * Nx) + + # calculate variance of spectral coefficient + sigma_2 = fkk**2 / Nx / Ny + + # Group ellipses of alphas into the same wavenumber k/Nmin + for k in range(1, Nmin + 1): + alpha = k / Nmin + alpha_p1 = (k + 1) / Nmin + + # Sum up elements matching k + mask_k = np.where((alpha_matrix >= alpha) & (alpha_matrix < alpha_p1)) + ps_array[t, k - 1] = np.sum(sigma_2[mask_k]) + + return ps_array + + +def _create_alpha_matrix(Ny, Nx): + """Construct an array of 2D wavenumbers from 2D wavenumber pair. + + Parameters + ---------- + Ny, Nx: + Dimensions of the 2D field for which the power spectra is calculated. Used to + create the array of 2D wavenumbers. Each Ny, Nx pair is associated with a + single-scale parameter. + + Returns: alpha_matrix + normalisation of 2D wavenumber axes, transforming the spectral domain into + an elliptic coordinate system. + + """ + # Create x_indices: each row is [1, 2, ..., Nx] + x_indices = np.tile(np.arange(1, Nx + 1), (Ny, 1)) + + # Create y_indices: each column is [1, 2, ..., Ny] + y_indices = np.tile(np.arange(1, Ny + 1).reshape(Ny, 1), (1, Nx)) + + # Compute alpha_matrix + alpha_matrix = np.sqrt((x_indices**2) / Nx**2 + (y_indices**2) / Ny**2) + + return alpha_matrix diff --git a/src/CSET/recipes/level_fields/generic_level_power_spectrum_series.yaml b/src/CSET/recipes/level_fields/generic_level_power_spectrum_series.yaml new file mode 100644 index 000000000..96a0abd0d --- /dev/null +++ b/src/CSET/recipes/level_fields/generic_level_power_spectrum_series.yaml @@ -0,0 +1,51 @@ +category: Power Spectrum +title: Power Spectrum $VARNAME $LEVELTYPE $LEVEL +description: | + Extracts and plots the power spectrum 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) + 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: plot.plot_power_spectrum_series + sequence_coordinate: $SEQUENCE + # stamp_coordinate and single_plot optional and only required for ensemble data + stamp_coordinate: "realization" + single_plot: False + + - operator: write.write_cube_to_nc + overwrite: True diff --git a/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series.yaml b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series.yaml new file mode 100644 index 000000000..596c20176 --- /dev/null +++ b/src/CSET/recipes/surface_fields/generic_surface_power_spectrum_series.yaml @@ -0,0 +1,45 @@ +category: Power Spectrum +title: Power Spectrum $VARNAME +description: | + Extracts and plots the power spectrum of surface $VARNAME + + 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) + + 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 + 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: write.write_cube_to_nc + overwrite: True + + - operator: plot.plot_power_spectrum_series + sequence_coordinate: $SEQUENCE diff --git a/tests/conftest.py b/tests/conftest.py index e98cf2186..7bc1c40bb 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -131,6 +131,36 @@ def histogram_cube(histogram_cube_readonly): return histogram_cube_readonly.copy() +@pytest.fixture(scope="session") +def field2d_cube_readonly(): + """Get a 2D Cube for testing power spectrum code. It is NOT safe to modify.""" + from CSET.operators import read + + return read.read_cube("tests/test_data/air_temperature_lat_lon.nc") + + +@pytest.fixture() +def field2d_cube(field2d_cube_readonly): + """Get a 2D cube for testing power spectrum code.""" + return field2d_cube_readonly.copy() + + +@pytest.fixture(scope="session") +def power_spectrum_cube_readonly(): + """Get a Cube for testing power spectrum code. It is NOT safe to modify.""" + from CSET.operators import read + + return read.read_cube( + "tests/test_data/power_spectrum_temperature_at_pressure_levels_pressure_250_1time.nc" + ) + + +@pytest.fixture() +def power_spectrum_cube(power_spectrum_cube_readonly): + """Get a Cube for testing power spectrum code.""" + return power_spectrum_cube_readonly.copy() + + @pytest.fixture(scope="session") def regrid_rectilinear_cube_readonly(): """Get a cube to test with. It is NOT safe to modify.""" diff --git a/tests/operators/test_plot.py b/tests/operators/test_plot.py index 6496c1c6c..2a86a0616 100644 --- a/tests/operators/test_plot.py +++ b/tests/operators/test_plot.py @@ -891,6 +891,155 @@ def test_plot_and_save_postage_stamps_in_single_plot_histogram_series( assert Path("test.png").is_file() +def test_plot_power_spectrum_with_filename(field2d_cube, tmp_working_dir): + """Testing power spectrum code produces file.""" + plot.plot_power_spectrum_series( + field2d_cube, filename="test", sequence_coordinate="time" + ) + assert Path("test_464569.0.png").is_file() + + +def test_plot_and_save_postage_stamp_power_spectrum_series( + power_spectrum_cube, tmp_working_dir +): + """Test plotting a postage stamp power spectrum.""" + plot._plot_and_save_postage_stamp_power_spectrum_series( + cube=power_spectrum_cube, + filename="test.png", + title="Test", + stamp_coordinate="realization", + histtype="step", + ) + assert Path("test.png").is_file() + + +def test_plot_and_save_postage_stamps_in_single_plot_power_spectrum_series( + power_spectrum_cube, tmp_working_dir +): + """Test plotting a multiline power spectrum for multiple ensemble members.""" + plot._plot_and_save_postage_stamps_in_single_plot_power_spectrum_series( + cube=power_spectrum_cube, + filename="test.png", + title="Test", + stamp_coordinate="realization", + histtype="step", + ) + assert Path("test.png").is_file() + + +def test_create_alpha_matrix_shape(): + """Test shape of alpha matrix used in power spectrum calculation.""" + Ny, Nx = 10, 15 + alpha = plot._create_alpha_matrix(Ny, Nx) + assert alpha.shape == (Ny, Nx), "Alpha matrix shape mismatch" + + +def test_create_alpha_matrix_values(): + """Test alpha matrix contains only positive values.""" + Ny, Nx = 4, 4 + alpha = plot._create_alpha_matrix(Ny, Nx) + assert np.all(alpha >= 0), "Alpha matrix contains negative values" + assert np.isclose(alpha[0, 0], np.sqrt((1 / Nx) ** 2 + (1 / Ny) ** 2)), ( + "Alpha matrix value incorrect" + ) + + +def test_dct_ps_output_shape(): + """Test shape of power spectrum output from _DCT_ps.""" + Nt, Ny, Nx = 5, 10, 10 + y_3d = np.random.rand(Nt, Ny, Nx) + ps = plot._DCT_ps(y_3d) + expected_shape = (Nt, min(Nx - 1, Ny - 1)) + assert ps.shape == expected_shape, "Power spectrum output shape mismatch" + + +def test_dct_ps_non_negative(): + """Test power spectrum only contains positive values.""" + Nt, Ny, Nx = 3, 8, 8 + y_3d = np.random.rand(Nt, Ny, Nx) + ps = plot._DCT_ps(y_3d) + assert np.all(ps >= 0), "Power spectrum contains negative values" + + +def test_dct_ps_known_input(): + """Test _DCT_ps produces non-zero spectrum for constant input.""" + # Use a constant field to test expected behavior + Nt, Ny, Nx = 2, 4, 4 + y_3d = np.ones((Nt, Ny, Nx)) + ps = plot._DCT_ps(y_3d) + assert np.allclose(ps[:, 1:], 0, atol=1e-6), "Non-zero spectrum for constant input" + + +def test_plot_power_spectrum_no_sequence_coordinate(field2d_cube, tmp_working_dir): + """Error when cube is missing sequence coordinate (time).""" + field2d_cube.remove_coord("time") + with pytest.raises(ValueError, match="Cube must have a time coordinate."): + plot.plot_power_spectrum_series(field2d_cube, series_coordinate="pressure") + + +def make_test_cube_power_spectrum(shape=(1, 10, 10), time_points=None): + """Create test cube for use with the power spectrum tests.""" + data = np.random.rand(*shape) + if time_points is None: + time_points = [0] + time_coord = iris.coords.DimCoord( + time_points, standard_name="time", units="hours since 1970-01-01 00:00:00" + ) + y_coord = iris.coords.DimCoord(np.arange(shape[1]), long_name="y", units="1") + x_coord = iris.coords.DimCoord(np.arange(shape[2]), long_name="x", units="1") + cube = iris.cube.Cube( + data, + dim_coords_and_dims=[(time_coord, 0), (y_coord, 1), (x_coord, 2)], + long_name="test_data", + ) + return cube + + +def test_calculate_power_spectrum_raises_for_bad_dim(): + """Check error is raised if the cube has too many dimensions.""" + cube_3d = make_test_cube_power_spectrum() + + # Add 2 new dimensions to cube_3d to make 5D + new_data = cube_3d.data[np.newaxis, np.newaxis, :, :, :] + + # Create dummy coordinates for the new dimensions + coord_0 = iris.coords.DimCoord([0], long_name="extra_dim_0") + coord_1 = iris.coords.DimCoord([0], long_name="extra_dim_1") + + # Build dim_coords_and_dims manually + dim_coords_and_dims = [(coord_0, 0), (coord_1, 1)] + for i, coord in enumerate(cube_3d.dim_coords): + dim_coords_and_dims.append((coord, i + 2)) # shift by 2 for new axes + + # Create the new 4D cube + cube_5d = iris.cube.Cube(new_data, dim_coords_and_dims=dim_coords_and_dims) + + if isinstance(cube_5d, iris.cube.CubeList): + cube_5d = cube_5d[0] + + with pytest.raises( + ValueError, match="Cube dimensions unsuitable for power spectra code" + ): + plot.plot_power_spectrum_series(cubes=cube_5d) + + +def test_calculate_power_spectrum_raises_for_bad_dim_1D(): + """Check error is raised if the cube has too few dimensions.""" + cube_3d = make_test_cube_power_spectrum() + + # Make a 1D field + + cube_1d = cube_3d.collapsed(["x", "y"], iris.analysis.MEAN) + + if isinstance(cube_1d, iris.cube.CubeList): + cube_1d = cube_1d[0] + + with pytest.raises( + ValueError, match="Cube dimensions unsuitable for power spectra code" + ): + plot.plot_power_spectrum_series(cubes=cube_1d) + + def test_scatter_plot(cube, vertical_profile_cube, tmp_working_dir): """Save a scatter plot.""" cube_y = collapse.collapse(cube, ["time", "grid_longitude"], "MEAN")[0:4] diff --git a/tests/test_data/power_spectrum_temperature_at_pressure_levels_pressure_250_1time.nc b/tests/test_data/power_spectrum_temperature_at_pressure_levels_pressure_250_1time.nc new file mode 100644 index 000000000..4d61db49e Binary files /dev/null and b/tests/test_data/power_spectrum_temperature_at_pressure_levels_pressure_250_1time.nc differ