diff --git a/src/CSET/operators/misc.py b/src/CSET/operators/misc.py index 9767b16f1..d3a9ae58f 100644 --- a/src/CSET/operators/misc.py +++ b/src/CSET/operators/misc.py @@ -19,6 +19,7 @@ from collections.abc import Iterable import iris +import iris.analysis.calculus import numpy as np from iris.cube import Cube, CubeList @@ -426,3 +427,43 @@ def convert_units(cubes: iris.cube.Cube | iris.cube.CubeList, units: str): return new_cubelist[0] else: return new_cubelist + + +def differentiate( + cubes: iris.cube.Cube | iris.cube.CubeList, coordinate: str, **kwargs +) -> iris.cube.Cube | iris.cube.CubeList: + """Differentiate a cube on a specified coordinate. + + Arguments + --------- + cubes: iris.cube.Cube | iris.cube.CubeList + A Cube or CubeList of a field that is to be differentiated. + + coordinate: str + The coordinate that is to be differentiated over. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + The differential of the cube along the specified coordinate. + + Notes + ----- + The differential is calculated based on a carteisan grid. This calculation + is then suitable for vertical and temporal derivatives. It is not sensible + for horizontal derivatives if they are based on spherical coordinates (e.g. + latitude and longitude). In essensce this operator is a CSET wrapper around + iris.analysis.calculus.differentiate. + + Examples + -------- + >>> dT_dz = misc.differentiate(temperature, "altitude") + """ + new_cubelist = iris.cube.CubeList([]) + for cube in iter_maybe(cubes): + dcube = iris.analysis.calculus.differentiate(cube, coordinate) + new_cubelist.append(dcube) + if len(new_cubelist) == 1: + return new_cubelist[0] + else: + return new_cubelist diff --git a/tests/operators/test_misc.py b/tests/operators/test_misc.py index 23dd1b6ee..3ff36e02f 100644 --- a/tests/operators/test_misc.py +++ b/tests/operators/test_misc.py @@ -17,6 +17,7 @@ import datetime import iris +import iris.analysis.calculus import iris.coords import iris.cube import iris.exceptions @@ -355,3 +356,24 @@ def test_convert_units_cubelist(cube): for actual, expected in zip(new_cubelist, expected_cubelist, strict=True): assert actual.units == expected.units assert np.allclose(actual.data, expected.data, rtol=1e-6, atol=1e-2) + + +def test_differentitate(vertical_profile_cube): + """Test a differentitation of a vertical profile cube.""" + expected_cube = iris.analysis.calculus.differentiate( + vertical_profile_cube, "pressure" + ) + actual_cube = misc.differentiate(vertical_profile_cube, coordinate="pressure") + assert np.allclose(actual_cube.data, expected_cube.data, rtol=1e-6, atol=1e-2) + + +def test_differentitate_cubelist(long_forecast): + """Test a differentiation of a CubeList.""" + # Create input CubeList. + input_cubes = iris.cube.CubeList([long_forecast, long_forecast]) + # Create expected cube and then convert to a CubeList + expectedcube = iris.analysis.calculus.differentiate(long_forecast, "time") + expected_cubelist = iris.cube.CubeList([expectedcube, expectedcube]) + new_cubelist = misc.differentiate(input_cubes, "time") + for actual, expected in zip(new_cubelist, expected_cubelist, strict=True): + assert np.allclose(actual.data, expected.data, rtol=1e-6, atol=1e-2)