Adds Contour Frequency by Altitude Diagrams (CFAD) calculation code

src/CSET/operators/plot.py

@@ -2003,6 +2003,73 @@ def _validate_cubes_coords(
         )
 
 
+def _calculate_CFAD(
+    cube: iris.cube.Cube, vertical_coordinate: str, bin_edges: list[float]
+) -> iris.cube.Cube:
+    """Calculate a Contour Frequency by Altitude Diagram (CFAD).
+
+    Parameters
+    ----------
+    cube: iris.cube.Cube
+        A cube of the data to be turned into a CFAD. It should be a minimum
+        of two dimensions with one being a user specified vertical coordinate.
+    vertical_coordinate: str
+        The vertical coordinate of the cube for the CFAD to be calculated over.
+    bin_edges: list[float]
+        The bin edges for the histogram. The bins need to be specified to
+        ensure consistency across the CFAD, otherwise it cannot be interpreted.
+
+    Notes
+    -----
+    Contour Frequency by Altitude Diagrams (CFADs) were first designed by
+    Yuter and Houze (1995)[YuterandHouze95]. They are calculated by binning the
+    data by altitude and then by variable bins (e.g. temperature). The variable
+    bins are then normalised by each altitude. This essenitally creates a
+    normalised frequency distribution for each altitude. These are then stacked
+    and combined in a single plot.
+
+    References
+    ----------
+    .. [YuterandHouze95] Yuter S.E., and Houze, R.A. (1995) "Three-Dimensional
+       Kinematic and Microphysical Evolution of Florida Cumulonimbus. Part II:
+       Frequency Distributions of Vertical Velocity, Reflectivity, and
+       Differential Reflectivity" Monthly Weather Review,  vol. 123, 1941-1963,
+       doi: 10.1175/1520-0493(1995)123<1941:TDKAME>2.0.CO;2
+    """
+    # Setup empty array for containing the CFAD data.
+    CFAD_values = np.zeros(
+        (len(cube.coord(vertical_coordinate).points), len(bin_edges) - 1)
+    )
+
+    # Calculate the CFAD as a histogram summing to one for each level.
+    for i, level_cube in enumerate(cube.slices_over(vertical_coordinate)):
+        # Note setting density to True does not produce the correct
+        # normalization for a CFAD, where each row must sum to one.
+        CFAD_values[i, :] = (
+            np.histogram(level_cube.data.reshape(level_cube.data.size), bins=bin_edges)[
+                0
+            ]
+            / level_cube.data.size
+        )
+    # Calculate central points for bins.
+    bins = (np.array(bin_edges[:-1]) + np.array(bin_edges[1:])) / 2.0
+    bin_bounds = np.array((bin_edges[:-1], bin_edges[1:])).T
+    # Now construct the coordinates for the cube.
+    vert_coord = cube.coord(vertical_coordinate)
+    bin_coord = iris.coords.DimCoord(
+        bins, bounds=bin_bounds, standard_name=cube.standard_name, units=cube.units
+    )
+    # Now construct the cube that is to be output.
+    CFAD = iris.cube.Cube(
+        CFAD_values,
+        dim_coords_and_dims=[(vert_coord, 0), (bin_coord, 1)],
+        long_name=f"{cube.name()}_cfad",
+        units="1",
+    )
+    CFAD.rename(f"{cube.name()}_cfad")
+    return CFAD
+
+
 ####################
 # Public functions #
 ####################

tests/conftest.py

@@ -323,3 +323,15 @@ def orography_4D_cube_read_only():
 def orography_4D_cube(orography_4D_cube_read_only):
     """Get 4D orography cube to run tests on. It is safe to modify."""
     return orography_4D_cube_read_only.copy()
+
+
+@pytest.fixture()
+def xwind_read_only():
+    """Get regridded xwind to run tests on. It is NOT safe to modify."""
+    return read.read_cube("tests/test_data/ageofair/aoa_in_rgd.nc", "x_wind")
+
+
+@pytest.fixture()
+def xwind(xwind_read_only):
+    """Get regridded xwind to run tests on. It is safe to modify."""
+    return xwind_read_only.copy()

tests/operators/test_ageofair.py

@@ -22,12 +22,6 @@
 from CSET.operators import ageofair, read
 
 
-@pytest.fixture()
-def xwind() -> iris.cube.Cube:
-    """Get regridded xwind to run tests on."""
-    return read.read_cube("tests/test_data/ageofair/aoa_in_rgd.nc", "x_wind")
-
-
 @pytest.fixture()
 def ywind() -> iris.cube.Cube:
     """Get regridded ywind to run tests on."""

tests/operators/test_plot.py

@@ -18,6 +18,7 @@
 import logging
 from pathlib import Path
 
+import cf_units
 import iris.coords
 import iris.cube
 import matplotlib as mpl
@@ -1253,3 +1254,44 @@ def test_qq_plot_grid_staggering_regrid(cube, tmp_working_dir):
         model_names=["a", "b"],
     )
     assert Path("untitled.png").is_file()
+
+
+def test_calculate_CFAD(xwind):
+    """Test calculating a CFAD."""
+    bins = np.array([-50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50])
+    calculated_CFAD = np.zeros((len(xwind.coord("pressure").points), len(bins) - 1))
+    for j, level_cube in enumerate(xwind.slices_over("pressure")):
+        calculated_CFAD[j, :] = (
+            np.histogram(level_cube.data.reshape(level_cube.data.size), bins=bins)[0]
+            / level_cube.data.size
+        )
+    assert np.allclose(
+        plot._calculate_CFAD(
+            xwind, "pressure", [-50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50]
+        ).data,
+        calculated_CFAD,
+        rtol=1e-06,
+        atol=1e-02,
+    )
+
+
+def test_name_CFAD(xwind):
+    """Test naming of CFAD cube."""
+    expected_name = f"{xwind.name()}_cfad"
+    assert (
+        plot._calculate_CFAD(
+            xwind, "pressure", [-50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50]
+        ).name()
+        == expected_name
+    )
+
+
+def test_units_CFAD(xwind):
+    """Test units of CFAD cube."""
+    expected_units = cf_units.Unit("1")
+    assert (
+        plot._calculate_CFAD(
+            xwind, "pressure", [-50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50]
+        ).units
+        == expected_units
+    )