Adds differentiation operator

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
 
@@ -480,3 +481,43 @@ def rename_cube(cubes: iris.cube.Cube | iris.cube.CubeList, name: 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 essence this operator is a CSET wrapper around
+    `iris.analysis.calculus.differentiate <https://scitools-iris.readthedocs.io/en/stable/generated/api/iris.analysis.calculus.html#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

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
@@ -407,3 +408,24 @@ def test_rename_cube_for_cubelist(cube):
     new_cubelist = misc.rename_cube(cube_list, "air_temperature_at_screen_level")
     for new in new_cubelist:
         assert new.name() == "air_temperature_at_screen_level"
+
+
+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)