diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index ab0415558..3387e4d99 100644 --- a/src/CSET/operators/plot.py +++ b/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 # #################### diff --git a/tests/conftest.py b/tests/conftest.py index 7bc1c40bb..8ba9d0c79 100644 --- a/tests/conftest.py +++ b/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() diff --git a/tests/operators/test_ageofair.py b/tests/operators/test_ageofair.py index a2fae03e1..b494b1807 100644 --- a/tests/operators/test_ageofair.py +++ b/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.""" diff --git a/tests/operators/test_plot.py b/tests/operators/test_plot.py index 1616edb5a..1e395ee2b 100644 --- a/tests/operators/test_plot.py +++ b/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 + )