Regional power spectra using Discrete Cosine Transform (DCT) #1765

cehalliwell · 2025-10-15T13:49:44Z

Follows on from #215

Calculate power spectra of 2D field in regional model using Discrete Fourier Transform (DCT)

Contribution checklist

Aim to have all relevant checks ticked off before merging. See the developer's guide for more detail.

Documentation has been updated to reflect change.
New code has tests, and affected old tests have been updated.
All tests and CI checks pass.
Ensured the pull request title is descriptive.
Conda lock files have been updated if dependencies have changed.
Attributed any Generative AI, such as GitHub Copilot, used in this PR.
Marked the PR as ready to review.

github-actions · 2025-10-15T13:51:45Z

cehalliwell · 2025-10-20T09:53:35Z

Power spectra calculated over regional (limited area) domain using Discrete Cosine Transform (DCT). This code is based on Denis et al 2002.

The code calculates power spectrum on a single (surface, pressure or model) level for each specified time/level but does not yet allow aggregation or consider ensembles.

Mean spectra over a specified time period for each level and/or mean spectra over a range of levels at a specified time may be future options.

M365 Copilot was used to optimise the code and give code suggestions.

cehalliwell · 2025-10-22T11:41:39Z

Power Spectrum code changed so that the power spectrum is now calculated in _plot_and_save_power_spectrum_series rather than a separate operator. This is so results from multiple models can be plotted on the same plot.

The unit tests for the power spectrum code have now been updated to reflect these changes.
Microsoft Copilot was used to help design the unit tests.

There are currently only 2 recipes:
level_fields/generic_level_power_spectrum_series.yaml
surface_fields/generic_surface_power_spectrum_series.yaml

which allow plotting of power spectrum on model, pressure and surface levels. There is no method of aggregation at the moment. This might involve taking mean of spectra over times or levels but needs further investigation.

cehalliwell · 2025-10-23T11:46:40Z

This is an example plot from CSET using this method. Currently, models must be on the same grid.

power_spectrum_vertical_wind_at_pressure_levels_pressure_250_459459 0

daflack

This PR is looking good. I've got a few comments here, and some more specific ones in the code.

It is worth considering tests for the following: 3D data, a cubelist, no coordinate named time present (if allowing it, see comments in code), a multi-model run, a test to cover the other options of L2649, L2692, and L2702of plot.py if those if statements are needed.

It would be nice to include somewhere (either in the function description, or the recipe description) a quick guide to interpretation of power spectra.

You also need to update the conda lock files too: Actions -> Update conda lock files -> Run workflow, select your branch, and then merge in the PR that the bot creates into your PR.

Finally, it might be worth getting someone else to do a technical review of this PR as well, or discuss with @jfrost-mo whether the power spectra code should be separate from the plotting code. My only thought behind this is that for aggregation purposes it might be easier to have the calculation of the power spectra separate.

It's really great to see progress on this, and the output looks fantastic.

daflack · 2025-10-27T11:24:47Z

src/CSET/loaders/power_spectrum.py

+                variables={
+                    "VARNAME": field,
+                    "MODEL_NAME": [model["name"] for model in models],
+                    "SEQUENCE": "time"


Applies to all three loaders here. Should the SEQUENCE be set to conf.SPECTRUM_SURFACE_FIELD_SEQUENCE? Otherwise won't this setting overwrite what is in there recipe and always do a sequence rather than the FALSE option? Or is this a temporary thing related to lack of aggregation?

This is a copy of the loader for histograms, basically replacing "HISTOGRAM" with "SPECTRUM" and removing anything relating to aggregation at the moment. Maybe I misunderstand what the code but I thought conf.SPECTRUM_SURFACE_FIELD_SEQUENCE was a logical and "SEQUENCE" was specifying the coordinate, which would always be time?

Do we actually need that logical then? I don't see it being used anywhere so does the logical come in with aggregation (i.e. aggregation across a model run) and so should be included later?

I'm not sure. Maybe a question for whoever wrote the loader for the histograms? I don't see it in the aggregation section of the histogram loader.

I suppose one thing to try is running it with and without that switch to False and see what happens. It may be more of a question for who put that functionality into CSET rather than the loader, as the loader was based on the original includes file.

It ran with that set to False when using cset bake. (?!)

Ir ran with the switch set to False. ?!?

daflack · 2025-10-27T11:25:08Z

src/CSET/loaders/power_spectrum.py

+                    "LEVELTYPE": "pressure",
+                    "LEVEL": plevel,
+                    "MODEL_NAME": [model["name"] for model in models],
+                    "SEQUENCE": "time"


See question above.

See comment above

daflack · 2025-10-27T11:25:20Z

src/CSET/loaders/power_spectrum.py

+                    "LEVELTYPE": "model_level_number",
+                    "LEVEL": mlevel,
+                    "MODEL_NAME": [model["name"] for model in models],
+                    "SEQUENCE": "time"


See question above.

See comment above

daflack · 2025-10-27T12:06:58Z

src/CSET/operators/plot.py

+        # If cube is a CubeList, extract the first cube
+        if isinstance(cube, iris.cube.CubeList):
+            cube = cube[0]


Suggested change

# If cube is a CubeList, extract the first cube

if isinstance(cube, iris.cube.CubeList):

cube = cube[0]

I don't think this is needed any more? cubes will be a cube/cubelist and iter_maybe will loop over each cube within the cubelist, so will produce this result. Or is it aiming at something else?

I think you are right - it was necessary in an earlier version of the code but I moved it into the for loop. I have removed that section now.

daflack · 2025-10-27T12:16:20Z

src/CSET/operators/plot.py

+        elif cube.ndim == 3:
+            cube_3d = cube.data
+        else:
+            raise ValueError("Cube dimensions unsuitable for power spectra code")


Suggested change

raise ValueError("Cube dimensions unsuitable for power spectra code")

raise ValueError(f"Cube is {cube.ndim} dimensional. Cube should be 2 or 3 dimensional.")

This error message is a little more useful as it tells the user how the cube is unsuitable and what needs to be done to fix the error.

I agree. I have added the extra message. Thanks.

daflack · 2025-10-27T12:34:15Z

src/CSET/operators/plot.py

+        if isinstance(cube_slice, iris.cube.CubeList):
+            single_cube = cube_slice[0]


A little confused around this section? The cubes may not always be read in the same order so you could end up changing cubes by selecting the first cube_slice. In what instances will the slice be a cubelist?

Again, this is copied from the histogram code. When running the code, this is loop is accessed (I checked with a print statement)

The code does not work without this and, checking with a print statement, this if statement is accessed, implying that it can be a cubelist.

daflack · 2025-10-27T12:36:40Z

src/CSET/operators/plot.py

+    return cubes
+
+
+def DCT_ps(y_3d):


Suggested change

def DCT_ps(y_3d):

def _DCT_ps(y_3d):

If we are not using this function within other recipes (i.e. only using it as part of a plotting function) we could define this as a non-public function.

Yes, that makes sense. Changing the function to non-public.

You are right. Function changed to non-public. Thanks.

daflack · 2025-10-27T12:37:09Z

src/CSET/operators/plot.py

+    # Regional domains:
+    # Calculate power spectra over limited are domain using Discrete Cosine Transform (DCT)
+    # as described in Denis et al 2002 [Denis_etal_2002]_.
+


Useful to include a description of y_3d.

Good point. I've added description of input (y_3d) and output (ps_array).

I have expanded the description in the _DCT_ps function.

daflack · 2025-10-27T12:39:03Z

src/CSET/operators/plot.py

+    return ps_array
+
+
+def create_alpha_matrix(Ny, Nx):


Suggested change

def create_alpha_matrix(Ny, Nx):

def _create_alpha_matrix(Ny, Nx):

Again I don't think this needs to be a public function as not used elsewhere.

Changed to non-public

I have made this function non-public.

daflack · 2025-10-27T12:39:31Z

src/CSET/operators/plot.py

+
+    Each pair is associated with a single-scale parameter. alpha_matrix is normalisation of
+    2D wavenumber axes, transforming the spectral domain into an elliptic coordinate system.
+    """


Useful to have a description of the Nx and Ny parameters.

Description expanded in _create_alpha_matrix

add functionality to calculate and plot spectra

324798a

cehalliwell self-assigned this Oct 15, 2025

cehalliwell added 2 commits October 15, 2025 14:52

correct line in operators.rst

12ecba1

unit tests for power spectrum code

c81877a

cehalliwell added 6 commits October 20, 2025 15:22

move power spectra operator into plot operator

8cbf4fd

Tidy up spurious comments

0db1864

remove any occurrence of power_spectrum operator (now removed)

9a8de34

Update unit tests for power spectra code and tidy up

93ab85c

Remove occurrences of deleted power_spectrum operator

e53f5d4

Remove redundant comments

c31fd48

Tidy up comments

ab55d05

cehalliwell requested a review from daflack October 23, 2025 11:42

Correct x-axis label for power spectra

8cc454c

daflack requested changes Oct 27, 2025

View reviewed changes

daflack added this to the CSET next milestone Oct 28, 2025

jfrost-mo self-requested a review October 31, 2025 16:18

	# If cube is a CubeList, extract the first cube
	if isinstance(cube, iris.cube.CubeList):
	cube = cube[0]

	raise ValueError("Cube dimensions unsuitable for power spectra code")
	raise ValueError(f"Cube is {cube.ndim} dimensional. Cube should be 2 or 3 dimensional.")

		if isinstance(cube_slice, iris.cube.CubeList):
		single_cube = cube_slice[0]

	def create_alpha_matrix(Ny, Nx):
	def _create_alpha_matrix(Ny, Nx):

Regional power spectra using Discrete Cosine Transform (DCT) #1765

Are you sure you want to change the base?

Regional power spectra using Discrete Cosine Transform (DCT) #1765

Uh oh!

Conversation

cehalliwell commented Oct 15, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Contribution checklist

Uh oh!

github-actions bot commented Oct 15, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

cehalliwell commented Oct 20, 2025

Uh oh!

cehalliwell commented Oct 22, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

cehalliwell commented Oct 23, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

daflack left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

cehalliwell commented Oct 15, 2025 •

edited

Loading

github-actions bot commented Oct 15, 2025 •

edited

Loading

cehalliwell commented Oct 22, 2025 •

edited

Loading

cehalliwell commented Oct 23, 2025 •

edited

Loading