From e85b47e524b19542a5ce2179d1918b423412f1c0 Mon Sep 17 00:00:00 2001 From: daflack Date: Thu, 19 Dec 2024 16:11:13 +0000 Subject: [PATCH 1/3] Adds basic CFAD calculation and tests Fixes #1010 --- src/CSET/operators/plot.py | 54 ++++++++++++++++++++++++++++++++++++ tests/operators/test_plot.py | 27 ++++++++++++++++++ 2 files changed, 81 insertions(+) diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index 7aaff96a8..f35ca7aca 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -2013,6 +2013,60 @@ 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. + """ + # Setup empty array for containing the CFAD data. + CFAD_values = np.zeros( + (len(cube.coord(vertical_coordinate).points), len(bin_edges) - 1) + ) + + # Set iterator for CFAD values. + i = 0 + + # Calculate the CFAD as a histogram summing to one for each level. + for level_cube in 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 + ) + i += 1 + # 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)], + standard_name=cube.standard_name, + units="1", + ) + CFAD.attributes["type"] = "Contour Frequency by Altitude Diagram (CFAD)" + return CFAD + + #################### # Public functions # #################### diff --git a/tests/operators/test_plot.py b/tests/operators/test_plot.py index 1616edb5a..1d36389ff 100644 --- a/tests/operators/test_plot.py +++ b/tests/operators/test_plot.py @@ -28,6 +28,12 @@ from CSET.operators import collapse, misc, plot, read +@pytest.fixture() +def xwind() -> iris.cube.Cube: + """Get regridded xwind to run tests on.""" + return iris.load_cube("tests/test_data/ageofair/aoa_in_rgd.nc", "x_wind") + + def test_check_single_cube(): """Conversion to a single cube, and rejection where not possible.""" cube = iris.cube.Cube([0.0]) @@ -1253,3 +1259,24 @@ 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)) + j = 0 + for level_cube in xwind.slices_over("pressure"): + calculated_CFAD[j, :] = ( + np.histogram(level_cube.data.reshape(level_cube.data.size), bins=bins)[0] + / level_cube.data.size + ) + j += 1 + 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, + ) From de036a8759c2ff2cc604ff2043b87375d05fb8bc Mon Sep 17 00:00:00 2001 From: daflack Date: Mon, 11 Aug 2025 13:39:38 +0100 Subject: [PATCH 2/3] Technical review corrections --- src/CSET/operators/plot.py | 28 ++++++++++++++++++++------- tests/conftest.py | 12 ++++++++++++ tests/operators/test_ageofair.py | 6 ------ tests/operators/test_plot.py | 33 +++++++++++++++++++++++--------- 4 files changed, 57 insertions(+), 22 deletions(-) diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index f35ca7aca..68458a114 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -2028,17 +2028,31 @@ def _calculate_CFAD( 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) ) - # Set iterator for CFAD values. - i = 0 - # Calculate the CFAD as a histogram summing to one for each level. - for level_cube in cube.slices_over(vertical_coordinate): + 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, :] = ( @@ -2047,8 +2061,7 @@ def _calculate_CFAD( ] / level_cube.data.size ) - i += 1 - # calculate central points for bins + # 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. @@ -2060,10 +2073,11 @@ def _calculate_CFAD( CFAD = iris.cube.Cube( CFAD_values, dim_coords_and_dims=[(vert_coord, 0), (bin_coord, 1)], + long_name=f"{cube.name()}_cfad", standard_name=cube.standard_name, units="1", ) - CFAD.attributes["type"] = "Contour Frequency by Altitude Diagram (CFAD)" + CFAD.rename(f"{cube.name()}_cfad") return CFAD 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 1d36389ff..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 @@ -28,12 +29,6 @@ from CSET.operators import collapse, misc, plot, read -@pytest.fixture() -def xwind() -> iris.cube.Cube: - """Get regridded xwind to run tests on.""" - return iris.load_cube("tests/test_data/ageofair/aoa_in_rgd.nc", "x_wind") - - def test_check_single_cube(): """Conversion to a single cube, and rejection where not possible.""" cube = iris.cube.Cube([0.0]) @@ -1265,13 +1260,11 @@ 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)) - j = 0 - for level_cube in xwind.slices_over("pressure"): + 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 ) - j += 1 assert np.allclose( plot._calculate_CFAD( xwind, "pressure", [-50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50] @@ -1280,3 +1273,25 @@ def test_calculate_CFAD(xwind): 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 + ) From a7737d0658537d242e2c5f8417d23b3d1b16845b Mon Sep 17 00:00:00 2001 From: daflack Date: Tue, 12 Aug 2025 09:40:17 +0100 Subject: [PATCH 3/3] Remove naming of standard name as renaming long name --- src/CSET/operators/plot.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index 68458a114..d0805a13d 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -2074,7 +2074,6 @@ def _calculate_CFAD( CFAD_values, dim_coords_and_dims=[(vert_coord, 0), (bin_coord, 1)], long_name=f"{cube.name()}_cfad", - standard_name=cube.standard_name, units="1", ) CFAD.rename(f"{cube.name()}_cfad")