diff --git a/.mailmap b/.mailmap index 8d9e1a0bdd..fc211aaab4 100644 --- a/.mailmap +++ b/.mailmap @@ -15,6 +15,7 @@ Caroline Sandford <35029690+cgsandford@users.noreply.github.com> Carwyn Pelley Chris Sampson <44228125+BelligerG@users.noreply.github.com> +David Johnston Daniel Mentiplay Eleanor Smith <40183561+ellesmith88@users.noreply.github.com> <40183561+ellesmith88@users.noreply.github.com> Fiona Rust diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index f8f84c05b5..57969579a3 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -58,6 +58,7 @@ below: - Timothy Hume (Bureau of Meteorology, Australia) - Katharine Hurst (Met Office, UK) - Simon Jackson (Met Office, UK) + - David Johnston (Met Office, UK) - Caroline Jones (Met Office, UK) - Peter Jordan (Met Office, UK) - Anzer Khan (Met Office, UK) diff --git a/improver/temperature/layer_mean_temperature.py b/improver/temperature/layer_mean_temperature.py new file mode 100644 index 0000000000..fbba872e62 --- /dev/null +++ b/improver/temperature/layer_mean_temperature.py @@ -0,0 +1,135 @@ +# (C) Crown Copyright, Met Office. All rights reserved. +# +# This file is part of 'IMPROVER' and is released under the BSD 3-Clause license. +# See LICENSE in the root of the repository for full licensing details. +import iris +import numpy as np +from cf_units import Unit +from iris.cube import Cube + +from improver import BasePlugin + + +class LayerTemperatureInterpolation(BasePlugin): + """ + Plugin to interpolate temperature values at specified layer boundaries. + + This plugin extracts all temperature levels within a specified vertical layer + (between `bottom` and `top` heights, in feet), and interpolates temperature + at the exact base and top of the layer. The output is a cube containing + temperature at all interior levels plus the interpolated base and top. + """ + + def process( + self, temp_cube: Cube, bottom: float, top: float, verbosity: int = 0 + ) -> Cube: + """ + Interpolate temperature values at layer boundaries. + + Args: + temp_cube: Input temperature cube with a height coordinate in metres. + bottom: Lower boundary of the layer in feet. + top: Upper boundary of the layer in feet. + verbosity: Verbosity level for printing debug information. + + Returns: + Cube containing temperature at all layer heights + (base, interior, and top). + """ + if bottom >= top: + raise ValueError(f"Bottom ({bottom} ft) must be less than top ({top} ft).") + + # Convert layer bounds from feet to metres using cf_units + bottom_m = Unit("ft").convert(bottom, "m") + top_m = Unit("ft").convert(top, "m") + + if verbosity: + print(f"Interpolating temperature at base: {bottom} ft") + + # Extract cube of temperature levels within layer + between_layer_temp_cube = temp_cube.extract( + iris.Constraint(height=lambda point: bottom_m < point < top_m) + ) + + # Interpolate temperature at base of layer + base_temp = temp_cube.interpolate( + [("height", np.array([bottom_m], dtype=np.float32))], + iris.analysis.Linear(), + collapse_scalar=False, + ) + + if verbosity: + print(f"Interpolating temperature at top: {top} ft") + + # Interpolate temperature at top of layer + top_temp = temp_cube.interpolate( + [("height", np.array([top_m], dtype=np.float32))], + iris.analysis.Linear(), + collapse_scalar=False, + ) + + # Merge cubes of temperature at top, bottom and within layer + cubes_to_merge = [base_temp, between_layer_temp_cube, top_temp] + cubes_to_merge = [cube for cube in cubes_to_merge if cube is not None] + layer_levels_temp_cube = iris.cube.CubeList(cubes_to_merge).concatenate_cube() + + return layer_levels_temp_cube + + +class CalculateLayerMeanTemperature(BasePlugin): + """Calculate the vertically weighted mean temperature for a layer.""" + + def process(self, layer_cube: Cube, verbosity: int = 0) -> Cube: + """ + Calculate the altitude-weighted mean temperature across the layer. + + Args: + layer_cube: Cube containing temperature at all heights within + the specified layer (including interpolated base and top). + verbosity: Set level of output to print. + + Returns: + 2D cube of layer mean temperature. + """ + # Set up array for holding sum of products of temperature and vertical distance + layer_temp_product = np.zeros(layer_cube.data.shape[1:]) + + # Estimate mean temperature of layers and + # Weight by vertical extent of layer + altitude_array = layer_cube.coord("height").points + for alt_index in range(1, len(altitude_array) - 1): + layer_thickness = ( + altitude_array[alt_index + 1] - altitude_array[alt_index - 1] + ) / 2 + layer_temp_product += layer_cube.data[alt_index, :, :] * layer_thickness + + # Add contributions from base and top + layer_temp_product += ( + layer_cube.data[0, :, :] * (altitude_array[1] - altitude_array[0]) / 2 + ) + layer_temp_product += ( + layer_cube.data[-1, :, :] * (altitude_array[-1] - altitude_array[-2]) / 2 + ) + + # Divide by total thickness to get mean + lmt_array = layer_temp_product / (altitude_array[-1] - altitude_array[0]) + + if verbosity: + print("Layer mean temperature array:", lmt_array) + + # Wrap result in a cube and add required metadata + lmt_cube = iris.cube.Cube( + lmt_array, + var_name="air_temperature", + units="K", + dim_coords_and_dims=( + (layer_cube.coord("projection_y_coordinate"), 0), + (layer_cube.coord("projection_x_coordinate"), 1), + ), + aux_coords_and_dims=( + (layer_cube.coord("forecast_period"), ()), + (layer_cube.coord("forecast_reference_time"), ()), + (layer_cube.coord("time"), ()), + ), + ) + return lmt_cube diff --git a/improver_tests/temperature/layer_mean_temperature/__init__.py b/improver_tests/temperature/layer_mean_temperature/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/improver_tests/temperature/layer_mean_temperature/test_calculate_layer_mean_temperature.py b/improver_tests/temperature/layer_mean_temperature/test_calculate_layer_mean_temperature.py new file mode 100644 index 0000000000..eb880ad347 --- /dev/null +++ b/improver_tests/temperature/layer_mean_temperature/test_calculate_layer_mean_temperature.py @@ -0,0 +1,171 @@ +# (C) Crown Copyright, Met Office. All rights reserved. +# +# This file is part of 'IMPROVER' and is released under the BSD 3-Clause license. +# See LICENSE in the root of the repository for full licensing details. +import iris +import numpy as np + +from improver.temperature.layer_mean_temperature import CalculateLayerMeanTemperature + + +def make_layer_cube(data): + """Create a simple 3D temperature cube with shape (height, y, x) for testing. + + Args: + data (np.ndarray): Temperature data with shape (height, y, x). + + Heights are in metres: + - 100 m + - 200 m + - 300 m + """ + # produces height arrays of the form [ 100.0, 200.0, 300.0, 400.0, ...] + height_points = np.array( + np.arange(start=100.0, stop=(data.shape[0] + 1) * 100, step=100), + dtype=np.float32, + ) + # produces y-coordinate arrays of the form [ 0.0, 1.0, 2.0, 3.0, ...] + y_points = np.array(range(data.shape[1]), dtype=np.float32) + # produces x-coordinate arrays of the form [ 0.0, 1.0, 2.0, 3.0, ...] + x_points = np.array(range(data.shape[2]), dtype=np.float32) + cube = iris.cube.Cube( + data, + dim_coords_and_dims=[ + (iris.coords.DimCoord(height_points, standard_name="height", units="m"), 0), + (iris.coords.DimCoord(y_points, "projection_y_coordinate", units="m"), 1), + (iris.coords.DimCoord(x_points, "projection_x_coordinate", units="m"), 2), + ], + var_name="air_temperature", + units="K", + ) + # Add scalar coordinates expected by CalculateLayerMeanTemperature + cube.add_aux_coord( + iris.coords.AuxCoord(0, standard_name="forecast_period", units="seconds") + ) + cube.add_aux_coord( + iris.coords.AuxCoord( + 0, standard_name="forecast_reference_time", units="seconds since epoch" + ) + ) + cube.add_aux_coord( + iris.coords.AuxCoord(0, standard_name="time", units="seconds since epoch") + ) + return cube + + +def test_layer_mean_temperature_uniform(): + """Test that CalculateLayerMeanTemperature returns the correct mean + for a uniform temperature field. + + If all temperature values are constant (280 K) at every height level, + the weighted mean must also be 280 K regardless of height spacing or weighting. + """ + # All data values set to 280 K across all heights and grid points + data = np.full((3, 2, 2), 280, dtype=np.float32) + cube = make_layer_cube(data) + plugin = CalculateLayerMeanTemperature() + result = plugin.process(cube, verbosity=0) + + # Mean of a uniform field must equal the uniform value + np.testing.assert_allclose(result.data, 280, rtol=1e-5) + + +def test_layer_mean_temperature_output_shape_and_coordinates_exists(): + """Test that CalculateLayerMeanTemperature returns a 2D cube with correct + dimension and auxiliary coordinates. + + Checks: + - Output cube is 2D (height dimension collapsed). + - projection_y_coordinate and projection_x_coordinate are present. + - Scalar coordinates (forecast_period, time, forecast_reference_time) are preserved. + """ + data = np.full((3, 2, 2), 280, dtype=np.float32) + cube = make_layer_cube(data) + plugin = CalculateLayerMeanTemperature() + result = plugin.process(cube, verbosity=0) + + # Check output is 2D - height dimension has been collapsed + assert result.shape == (2, 2) + + # Check dimension coordinates are present + assert result.coords("projection_y_coordinate") + assert result.coords("projection_x_coordinate") + + # Check height coordinate is no longer a dimension coordinate + assert not result.coords("height", dim_coords=True) + + # Check scalar coordinates are preserved + assert result.coords("forecast_period") + assert result.coords("forecast_reference_time") + assert result.coords("time") + + # Check mean values are as expected + assert np.allclose(result.data, 280) + + +def test_layer_mean_temperature_scalar_coordinate_values_preserved(): + """Test that CalculateLayerMeanTemperature preserves scalar coordinate + values from the input cube in the output cube. + + Checks that forecast_period, forecast_reference_time and time coordinate + values are unchanged between input and output. + """ + data = np.full((3, 2, 2), 280, dtype=np.float32) + cube = make_layer_cube(data) + plugin = CalculateLayerMeanTemperature() + result = plugin.process(cube, verbosity=0) + + # Check scalar coordinate values are preserved from input to output + assert ( + result.coord("forecast_period").points == cube.coord("forecast_period").points + ) + assert ( + result.coord("forecast_reference_time").points + == cube.coord("forecast_reference_time").points + ) + assert result.coord("time").points == cube.coord("time").points + + +def test_layer_mean_temperature_values(): + """Test that CalculateLayerMeanTemperature returns + the correct mean temperature values where the temperature + field is varying + + Reference data was collected from a prior run. + + """ + + # create a varying temperature field, so the resulting + # layer mean temperature will have varying values too + # the input is deterministically calculated (i.e. no random numbers) + data = np.fromfunction( + lambda layer, y, x: 280.0 + 10.0 * np.sin(x) - 5.0 * np.cos(y) + np.sqrt(layer), + (3, 2, 2), + dtype=np.float32, + ) + + cube = make_layer_cube(data) + plugin = CalculateLayerMeanTemperature() + + actual_result = plugin.process(cube, verbosity=0) + + expected_result = [ + [275.8535546875, 284.26826171875], + [278.15205078125, 286.5667529296875], + ] + + assert np.allclose(actual_result.data, expected_result) + + +def test_verbosity_layer_mean_temperature(capsys): + """Test that CalculateLayerMeanTemperature prints expected output + when verbosity is set to 1. + + Checks that the layer mean temperature array is printed to stdout. + """ + data = np.full((3, 2, 2), 280, dtype=np.float32) + cube = make_layer_cube(data) + plugin = CalculateLayerMeanTemperature() + plugin.process(cube, verbosity=1) + captured = capsys.readouterr() + assert "Layer mean temperature array" in captured.out diff --git a/improver_tests/temperature/layer_mean_temperature/test_layer_temperature_interpolation.py b/improver_tests/temperature/layer_mean_temperature/test_layer_temperature_interpolation.py new file mode 100644 index 0000000000..57cfd0d28e --- /dev/null +++ b/improver_tests/temperature/layer_mean_temperature/test_layer_temperature_interpolation.py @@ -0,0 +1,209 @@ +# (C) Crown Copyright, Met Office. All rights reserved. +# +# This file is part of 'IMPROVER' and is released under the BSD 3-Clause license. +# See LICENSE in the root of the repository for full licensing details. +import iris +import numpy as np +import pytest + +from improver.temperature.layer_mean_temperature import LayerTemperatureInterpolation + +METRES_TO_FT = 3.28084 + + +def make_test_cube(): + """Create a simple 3D temperature cube with shape (height, y, x) for testing. + + Heights are in metres: + - 200 m (~656 ft) + - 220 m (~722 ft) + - 240 m (~787 ft) + """ + data = np.array( + [ + [[280, 281], [282, 283]], # height = 200 m + [[285, 286], [287, 288]], # height = 220 m + [[290, 291], [292, 293]], # height = 240 m + ], + dtype=np.float32, + ) + height_points = np.array([200, 220, 240], dtype=np.float32) + y_points = np.array([0, 1], dtype=np.float32) + x_points = np.array([0, 1], dtype=np.float32) + cube = iris.cube.Cube( + data, + dim_coords_and_dims=[ + (iris.coords.DimCoord(height_points, standard_name="height", units="m"), 0), + (iris.coords.DimCoord(y_points, "projection_y_coordinate", units="m"), 1), + (iris.coords.DimCoord(x_points, "projection_x_coordinate", units="m"), 2), + ], + var_name="air_temperature", + units="K", + ) + return cube + + +def test_extract_and_interpolate_layer(): + """Test that LayerTemperatureInterpolation correctly extracts and interpolates + temperature levels within a specified layer. + + Layer bounds: 600 ft (≈182.88 m) to 800 ft (≈243.84 m). + + Expected output: + - Interpolated base at 600 ft (≈182.88 m) + - Interior levels at 200 m, 220 m, 240 m + - Interpolated top at 800 ft (≈243.84 m) + - Total: 5 height levels, shape (5, 2, 2) + """ + cube = make_test_cube() + plugin = LayerTemperatureInterpolation() + # bottom=600 ft converts to ~182.88 m, top=800 ft converts to ~243.84 m + # All three data levels (200, 220, 240 m) fall within the layer bounds, + # so the result should contain 5 heights: interpolated base + 3 interior + interpolated top + result = plugin.process(cube, bottom=600, top=800, verbosity=0) + # Should include interpolated base (600ft), all interior levels, and interpolated top (800ft) + # In this synthetic example, all heights are present, so result should have 5 heights + assert result.shape == (5, 2, 2) + # Check height coordinate values + expected_heights = np.array([182.88, 200, 220, 240, 243.84]) + np.testing.assert_allclose( + result.coord("height").points, expected_heights, rtol=1e-2 + ) + + +def test_extract_and_interpolate_layer_coordinates(): + """Test that LayerTemperatureInterpolation preserves and returns + the correct coordinates in the output cube. + + Checks: + - Height coordinate exists with correct units. + - projection_y_coordinate exists and matches input. + - projection_x_coordinate exists and matches input. + """ + cube = make_test_cube() + plugin = LayerTemperatureInterpolation() + result = plugin.process(cube, bottom=600, top=800, verbosity=0) + + # Check height coordinate exists and has correct units + height_coord = result.coord("height") + assert height_coord.units == "m" + + # Check y coordinate exists and matches input + np.testing.assert_array_equal( + result.coord("projection_y_coordinate").points, + cube.coord("projection_y_coordinate").points, + ) + + # Check x coordinate exists and matches input + np.testing.assert_array_equal( + result.coord("projection_x_coordinate").points, + cube.coord("projection_x_coordinate").points, + ) + + +def test_interpolated_base_and_top_values(): + """Test that LayerTemperatureInterpolation correctly interpolates + temperature values at the base and top of the layer. + + Layer bounds: 600 ft (≈182.88 m) to 800 ft (≈243.84 m). + + Expected values are derived directly from iris.analysis.Linear() interpolation + to ensure we are testing our plugin's behaviour, not reimplementing the science. + """ + cube = make_test_cube() + plugin = LayerTemperatureInterpolation() + result = plugin.process(cube, bottom=600, top=800, verbosity=0) + + # Use Iris directly to get the expected interpolated values at base and top + expected_base = cube.interpolate( + [("height", np.array([600 / METRES_TO_FT], dtype=np.float32))], + iris.analysis.Linear(), + ) + expected_top = cube.interpolate( + [("height", np.array([800 / METRES_TO_FT], dtype=np.float32))], + iris.analysis.Linear(), + ) + + # Check base temperature values match Iris interpolation at all grid points + np.testing.assert_allclose( + result.data[0, :, :], expected_base.data[0, :, :], rtol=1e-5 + ) + + # Check top temperature values match Iris interpolation at all grid points + np.testing.assert_allclose( + result.data[-1, :, :], expected_top.data[0, :, :], rtol=1e-5 + ) + + +def test_verbosity_layer_extraction(capsys): + """Test that LayerTemperatureInterpolation prints expected output + when verbosity is set to 1. + + Checks that the bottom and top layer bounds are printed to stdout. + """ + cube = make_test_cube() + plugin = LayerTemperatureInterpolation() + plugin.process(cube, bottom=600, top=800, verbosity=1) + captured = capsys.readouterr() + assert "600" in captured.out + assert "800" in captured.out + + +def test_layer_bounds_match_data_level(): + """Test that LayerTemperatureInterpolation handles the case where + the layer bounds are both very close to an existing data level. + + Layer bounds: 656 ft (≈200 m) to 787 ft (≈240 m). + Both 200 m and 240 m are exact data levels in the test cube. + + Expected output: + - No duplicate height points in the output cube. + """ + cube = make_test_cube() + plugin = LayerTemperatureInterpolation() + + # Pass bounds directly in feet - the plugin handles the conversion internally + result = plugin.process(cube, bottom=656, top=787, verbosity=0) + + # Check no duplicate height points exist in the output + height_points = result.coord("height").points + assert len(height_points) == len( + np.unique(height_points) + ), f"Duplicate height points found: {height_points}" + + +def test_no_interior_levels(): + """Test that LayerTemperatureInterpolation returns at least the interpolated + base and top when no interior levels exist within the layer bounds. + + Layer bounds: 610 ft (≈185 m) to 640 ft (≈195 m). + No data levels fall between 185 m and 195 m in the test cube + (nearest levels are 200 m, 220 m, 240 m). + + Expected output: + - Interpolated base at 610 ft (≈185 m) + - Interpolated top at 640 ft (≈195 m) + - Total: 2 height levels, shape (2, 2, 2) + """ + cube = make_test_cube() + plugin = LayerTemperatureInterpolation() + + # These bounds have no data levels between them + result = plugin.process(cube, bottom=610, top=640, verbosity=0) + + # Should have exactly 2 levels: interpolated base and top only + assert result.shape == (2, 2, 2) + + +def test_bottom_greater_than_top_raises_error(): + """Test that LayerTemperatureInterpolation raises a ValueError when + the bottom bound is greater than the top bound. + + This is physically nonsensical and should be caught early with a + clear error message rather than producing garbage output silently. + """ + cube = make_test_cube() + plugin = LayerTemperatureInterpolation() + + with pytest.raises(ValueError, match="Bottom .* must be less than top"): + plugin.process(cube, bottom=800, top=600, verbosity=0)