From 15a1e0ebe23f2b710cbf6ac9befd58e26f1dc409 Mon Sep 17 00:00:00 2001 From: mo-DavidJohnJohnston Date: Wed, 1 Apr 2026 22:43:39 +0100 Subject: [PATCH 1/9] revert to 6 out of 7 passes --- .../psychrometric_calculations.py | 48 +++ .../test_HumidityMixingRatio.py | 305 +++++++++++++++++- 2 files changed, 352 insertions(+), 1 deletion(-) diff --git a/improver/psychrometric_calculations/psychrometric_calculations.py b/improver/psychrometric_calculations/psychrometric_calculations.py index 70d52e9dbe..6c312f4e0f 100644 --- a/improver/psychrometric_calculations/psychrometric_calculations.py +++ b/improver/psychrometric_calculations/psychrometric_calculations.py @@ -7,6 +7,7 @@ import functools from typing import List, Optional, Tuple, Union +import iris import iris._constraints import numpy as np from iris.cube import Cube, CubeList @@ -460,6 +461,42 @@ def qsat_differential(qs, t, q, p): humidity[sub_saturated] = humidity_in[sub_saturated] return temperature, humidity +def get_pressure_points(cube: Cube) -> np.ndarray: + """ + Get the pressure points from a _on_pressure_levels cube. + :param cube: input cube + :return: pressure points + """ + for coord in cube.dim_coords: + if "pressure" in coord.name().lower(): + return coord.points + return np.array([]) + +def flip_cube_in_place(cube: Cube, axis:int=0) -> None: + """ + flip an Iris cube in-place along the specified axis + flips bounds if present + + :param cube: input cube to modify in-place + :param axis: axis along which to flip the cube + :return: None + + """ + # ---- Flip data in-place ---- + cube.data[:] = np.flip(cube.data, axis=axis) + + # ---- Find the dimension coord for the axis ---- + coord = cube.dim_coords[axis] + + # ---- Flip coord points in-place ---- + coord.points = coord.points[::-1] # 1D reversal + + # ---- Flip coord bounds if present ---- + if coord.bounds is not None: + # bounds flip is more complex than expected - flip on both axes + # as iris bounds are [low, high] or [high, low] + # depending on direction of points + coord.bounds = np.flip(np.flip(coord.bounds,axis=0),axis=1) class HumidityMixingRatio(BasePlugin): """Returns the humidity mass mixing ratio from temperature, pressure and relative humidity""" @@ -535,6 +572,17 @@ def generate_pressure_cube(self, temperature_cube) -> Cube: try: pressure_cube = expanded_pressure_list.concatenate_cube() + """ + the Iris concatenate_cube function can reverse the list order when forming the cube + so the pressure cube is flipped vertically compared to the temperature cube + check if this is the case and then re-flip + """ + pressure_points_for_pressure = get_pressure_points(pressure_cube) + pressure_points_for_temperature = get_pressure_points(temperature_cube) + flip_required = np.array_equal(np.flip(pressure_points_for_pressure), pressure_points_for_temperature) + if flip_required: + flip_cube_in_place(pressure_cube, axis=0) + except iris.exceptions.ConcatenateError as error: raise RuntimeError( "Unable to concatenate pressure cubelist with input ", diff --git a/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py b/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py index a68ac932c8..458c238f9f 100644 --- a/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py +++ b/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py @@ -6,15 +6,28 @@ from unittest.mock import patch, sentinel -import numpy as np import pytest from iris.coords import AncillaryVariable from iris.cube import Cube +import iris +import numpy as np +import iris.cube as icube + +from typing import List, Tuple +from copy import deepcopy from improver.metadata.constants.attributes import MANDATORY_ATTRIBUTES from improver.psychrometric_calculations.psychrometric_calculations import ( HumidityMixingRatio, + get_pressure_points, + flip_cube_in_place + +) +from improver.psychrometric_calculations.precipitable_water import ( + PrecipitableWater, ) + +from improver.constants import EARTH_SURFACE_GRAVITY_ACCELERATION, WATER_DENSITY from improver.synthetic_data.set_up_test_cubes import set_up_variable_cube LOCAL_MANDATORY_ATTRIBUTES = { @@ -166,6 +179,296 @@ def test_zero_humidity( assert np.isclose(result.data, expected, atol=1e-7).all() +def make_pressure_cube(temp_cube: Cube) -> Cube: + """ + Create a 3D pressure cube from a temperature_on_pressure_levels cube. + The resulting cube has shape (levels, y, x). + + :param temp_cube: input temperature on pressure cube + :return: a 3D pressure cube + + """ + + # ---------------------------------------- + # 1. Extract the 1D pressure coordinate + # ---------------------------------------- + p_coord = temp_cube.coord("pressure") # DimCoord + p_vals = p_coord.points # 1D array (nz,) + + # ---------------------------------------- + # 2. Broadcast pressure to match the cube grid + # ---------------------------------------- + + p_3d = p_vals[:, None, None] # (nz, 1, 1) + p_3d = np.broadcast_to(p_3d, temp_cube.shape) # (nz, ny, nx) + + # ---------------------------------------- + # 3. Build a new pressure cube + # ---------------------------------------- + pressure_cube = icube.Cube( + p_3d, + standard_name="air_pressure", + long_name="air_pressure unit test", + units=p_coord.units, + dim_coords_and_dims=[ + (p_coord, 0), # vertical dimension + (temp_cube.coord(axis='y'), 1), # y coordinate + (temp_cube.coord(axis='x'), 2), # x coordinate + ], + attributes = temp_cube.attributes, + ) + + return pressure_cube + +def set_up_temperature_cube(shape: Tuple[int], temperature_value: float, vertical_levels: List[float]) -> Cube: + """ + Create a temperature on pressure cube. + :param: shape: Shape of the temperature cube + :param: temperature_value: temperature value + :param: vertical_levels: List of vertical levels + :return: a temperature cube + + """ + temperature = set_up_variable_cube( + np.full(shape, temperature_value, dtype=np.float32), + "latlon", + name="air_temperature", + x_grid_spacing = 1.0, + y_grid_spacing = 1.0, + vertical_levels = vertical_levels, + pressure = True + ) + add_attribute_dictionary(temperature) + return temperature + +def set_up_rel_humidity_cube(shape: Tuple[int], rel_humidity_value: float, vertical_levels: List[float]) -> Cube: + """ + Create a relative humidity on pressure cube. + :param: shape: Shape of the relative humidity cube + :param: rel_humidity_value: relative humidity value + :param: vertical_levels: List of vertical levels + :return: a relative humidity cube + """ + rel_humidity = set_up_variable_cube( + np.full(shape, rel_humidity_value, dtype=np.float32), + "latlon", + name="relative_humidity", + units="1", + x_grid_spacing=1.0, + y_grid_spacing=1.0, + vertical_levels=vertical_levels, + pressure=True + ) + add_attribute_dictionary(rel_humidity) + return rel_humidity + +def test_get_pressure_points() -> None: + """ + tests function "get_pressure_points" which is a support function + written to check if a pressure cube has been inadvertantly flipped + within the Improver implementation of PrecipitableWater + + :return: None + """ + temperature_value, pressure_value, rel_humidity_value, expected = (293, 100000, 0.1, 1.459832e-3) + + # set up cubes + vertical_levels = [100000.0, 50000.0, 100.0] + shape = (len(vertical_levels),3,3) + + temperature = set_up_temperature_cube(shape, temperature_value, vertical_levels) + pressure = make_pressure_cube(temperature) + rel_humidity = set_up_rel_humidity_cube(shape, rel_humidity_value, vertical_levels) + + assert np.allclose( get_pressure_points(temperature), np.array(vertical_levels) ) + assert np.allclose( get_pressure_points(pressure), np.array(vertical_levels) ) + assert np.allclose( get_pressure_points(rel_humidity), np.array(vertical_levels) ) + + # check captialisation has no affect + rel_humidity.coord('pressure').rename("Pressure") + assert np.allclose(get_pressure_points(rel_humidity), np.array(vertical_levels)) + + # check null result when no "pressure" dimension + rel_humidity.coord('Pressure').rename("Pr3ssure") + assert np.allclose(get_pressure_points(rel_humidity), np.array([])) + +def add_attribute_dictionary(cube: Cube) -> None: + """ + Adds attributes dictionary to cube attributes + to allow pre-existing checking function "metadata_ok" + to be used. + + :param cube: Cube to add attributes to + :return: None + """ + # set up meta-data required by testing + attributes_dictionary = { + "title":"unit test data", + "source":"unit test", + "institution":"somewhere", + "least_significant_digit": 4 + } + for k, v in attributes_dictionary.items(): + cube.attributes[k] = v + +def test_mixing_ratio_without_pressure_parameter() -> None: + """ + the HumidityMixingRatio calculation will generate its own pressure cube + if one is not supplied. This unit tests verifies that the results are the + same with/without an explicit pressure parameter. + + This ticket reports values for the total precipitable water (TPW) being far too high. + + https://metoffice.atlassian.net/browse/EPPT-3209 + + The reason was that HumidityMixingRadio generated a pressure cube that was wrongly + flipped veritically. This unit test re-creates the failing scenario to test and + exercise the bug fix + + The unit test then does a very simple total precipitable water calculation + ensuring the output from HumidityMixingRadio is suitable. + + The improver calculation is then compared against a DIY calculation as a sanity check. + + """ + iris.FUTURE.save_split_attrs = True # to stop Iris warning + temperature_value, pressure_value, rel_humidity_value, expected = (293, 100000, 0.1, 1.459832e-3) + + # set up input cubes + vertical_levels = [100000.0, 50000.0, 100.0] + shape = (len(vertical_levels),3,3) + + temperature = set_up_temperature_cube(shape, temperature_value, vertical_levels) + pressure = make_pressure_cube(temperature) + rel_humidity = set_up_rel_humidity_cube(shape, rel_humidity_value, vertical_levels) + + # mixing ratio calculation with 3 parameters + w3 = HumidityMixingRatio()([temperature, pressure, rel_humidity]) + metadata_ok(w3, temperature) # asserts in function call + # check results on single layer are as expected where pressure is 100000 Pa + assert np.isclose(w3.data[0], expected, atol=1e-7).all() + + # mixing ratio calculation with 2 parameters + w2 = HumidityMixingRatio()([temperature, rel_humidity]) + metadata_ok(w2, temperature) # asserts in function call + + # check 2 parameter calculation gives same results as 3 parameter calculation + assert np.isclose(w3.data, w2.data).all() + + # use w3 to calculate precipitable water + pw = PrecipitableWater().process(w3) + + # perform integration step (summing water in vertical atmosphere column) for Improver + improver_tpw = np.sum(pw.data, axis=0) + + # perform DIY TPW calculation + delta = pressure.data[1:,:,:] - pressure.data[:-1,:,:] + mid_w = ( w3.data[1:,:,:] + w3.data[:-1,:,:] ) / 2.0 + integral_terms = (delta * mid_w) + unit_test_tpw_1 = -np.sum(integral_terms, axis=0) / (EARTH_SURFACE_GRAVITY_ACCELERATION * WATER_DENSITY) + # numpy's trapezium rule integration + unit_test_tpw_2 = -np.trapezoid(w3.data, x=pressure.data, axis=0) / (EARTH_SURFACE_GRAVITY_ACCELERATION * WATER_DENSITY) + + # verify DIY TPW integrations produce same results + # N.B. Improver uses np.trapezoid + assert np.isclose(unit_test_tpw_1.data, unit_test_tpw_2.data).all() + + # note for such a small cube the DIY and Improver calculations + # for TPW are somewhat different. + # one presumes they will converge for cubes with more cells + # give a 200% latitude for testing + + assert np.isclose(improver_tpw, unit_test_tpw_1, rtol=2).all() + + +def generate_cube_bounds_for_axis(cube : Cube, axis: int) -> np.ndarray: + """ + automatically generates a suitable set of cube bounds for an axis + support function to test "flip_cube_in_place" when the cube has bounds + + :param cube: Cube to generate bounds for + :param axis: axis to generate bounds for + :return: numpy array of bounds of shape (Nz,2) + """ + points = cube.dim_coords[axis].points + delta = points[1:] - points[:-1] + delta = np.concatenate((delta[:1], delta, delta[-1:] )) + # N.B. delta can be negative so "low" and "high" are nominal + low_bounds = points - delta[:-1] / 2 + high_bounds = points + delta[1: ] / 2 + np_bounds = np.transpose(np.vstack((low_bounds, high_bounds))) + # N.B. bounds are not necessarily [low,high] but will be rectified by Iris to follow ordering of points so can be [high, low] + # i.e. the following line to give a genuine [low, high] ordering is not required and may be confusingly over-ridden in any case + # np_bounds = np.sort(np_bounds, axis=-1) + return np_bounds + +@pytest.mark.parametrize( + "bounds", + ( + (False), + (True), + ), +) +def test_flip_cube_in_place(bounds: bool) -> None: + """ + this tests support function "flip_cube_in_place" for implementation of + Improver's PrecipitableWater. It is used when a pressure cube + is found to be inadvertantly flipped in the Z-dimension. + + In this test, a cube is flipped once to check the data, and other components + have individually been flipped. + + Then the cube is flipped again to check we come back to the starting point + i.e. the operation is its own inverse. + + :param bounds: if True then the cube tested will have bounds + if False then the cube tested will not have bounds + + :return: None + """ + # to pacify an Iris warning + iris.FUTURE.date_microseconds = True + + axis = 0 + temperature_value, pressure_value, rel_humidity_value, expected = (293, 100000, 0.1, 1.459832e-3) + vertical_levels = [100000.0, 50000.0, 100.0] + shape = (len(vertical_levels),3,3) + temperature = set_up_temperature_cube(shape, temperature_value, vertical_levels) + temperature.data[:] = np.random.uniform(273, 280, size=temperature.data.shape) + if bounds: + # set up pressure bounds for temperature on pressure cube + b = generate_cube_bounds_for_axis(temperature, axis) + temperature.dim_coords[axis].bounds = b + + # first flip + temperature_before_first_flip = deepcopy(temperature) + flip_cube_in_place(temperature) + + # check data flipped + np.allclose(np.flip(temperature_before_first_flip.data, axis=axis), temperature.data) + + # check coordinate points flipped + assert list(reversed(temperature.dim_coords[axis].points)) == list(temperature_before_first_flip.dim_coords[axis].points) + + # check bounds flipped + if bounds: + # performing an a, b comparison + a = temperature_before_first_flip.dim_coords[axis].bounds + b = temperature.dim_coords[axis].bounds + # note that the bounds flip is unexpectedly complex + # flipping on both axes is required. + a_f = np.flip( np.flip(a, axis=0), axis= 1) + + assert np.allclose( a_f, b ) + + # second flip + flip_cube_in_place(temperature) + + # cube equality is generally not a good idea - but the flip transform is simple + exactly_equal = temperature == temperature_before_first_flip + assert exactly_equal + + def test_height_levels(): """Check that the plugin works with height level data""" From 5d0206200fe35996d722ccf03a414208e0874280 Mon Sep 17 00:00:00 2001 From: mo-DavidJohnJohnston Date: Wed, 1 Apr 2026 23:01:05 +0100 Subject: [PATCH 2/9] ruff formatting --- .../psychrometric_calculations.py | 15 +- .../test_HumidityMixingRatio.py | 138 ++++++++++-------- 2 files changed, 90 insertions(+), 63 deletions(-) diff --git a/improver/psychrometric_calculations/psychrometric_calculations.py b/improver/psychrometric_calculations/psychrometric_calculations.py index 6c312f4e0f..10fd039332 100644 --- a/improver/psychrometric_calculations/psychrometric_calculations.py +++ b/improver/psychrometric_calculations/psychrometric_calculations.py @@ -461,6 +461,7 @@ def qsat_differential(qs, t, q, p): humidity[sub_saturated] = humidity_in[sub_saturated] return temperature, humidity + def get_pressure_points(cube: Cube) -> np.ndarray: """ Get the pressure points from a _on_pressure_levels cube. @@ -472,7 +473,8 @@ def get_pressure_points(cube: Cube) -> np.ndarray: return coord.points return np.array([]) -def flip_cube_in_place(cube: Cube, axis:int=0) -> None: + +def flip_cube_in_place(cube: Cube, axis: int = 0) -> None: """ flip an Iris cube in-place along the specified axis flips bounds if present @@ -496,7 +498,8 @@ def flip_cube_in_place(cube: Cube, axis:int=0) -> None: # bounds flip is more complex than expected - flip on both axes # as iris bounds are [low, high] or [high, low] # depending on direction of points - coord.bounds = np.flip(np.flip(coord.bounds,axis=0),axis=1) + coord.bounds = np.flip(np.flip(coord.bounds, axis=0), axis=1) + class HumidityMixingRatio(BasePlugin): """Returns the humidity mass mixing ratio from temperature, pressure and relative humidity""" @@ -575,11 +578,13 @@ def generate_pressure_cube(self, temperature_cube) -> Cube: """ the Iris concatenate_cube function can reverse the list order when forming the cube so the pressure cube is flipped vertically compared to the temperature cube - check if this is the case and then re-flip + check if this is the case and then re-flip """ - pressure_points_for_pressure = get_pressure_points(pressure_cube) + pressure_points_for_pressure = get_pressure_points(pressure_cube) pressure_points_for_temperature = get_pressure_points(temperature_cube) - flip_required = np.array_equal(np.flip(pressure_points_for_pressure), pressure_points_for_temperature) + flip_required = np.array_equal( + np.flip(pressure_points_for_pressure), pressure_points_for_temperature + ) if flip_required: flip_cube_in_place(pressure_cube, axis=0) diff --git a/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py b/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py index 458c238f9f..4345601970 100644 --- a/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py +++ b/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py @@ -4,30 +4,27 @@ # See LICENSE in the root of the repository for full licensing details. """Tests for the HumidityMixingRatio plugin""" +from copy import deepcopy +from typing import List, Tuple from unittest.mock import patch, sentinel +import iris +import iris.cube as icube +import numpy as np import pytest from iris.coords import AncillaryVariable from iris.cube import Cube -import iris -import numpy as np -import iris.cube as icube - -from typing import List, Tuple -from copy import deepcopy +from improver.constants import EARTH_SURFACE_GRAVITY_ACCELERATION, WATER_DENSITY from improver.metadata.constants.attributes import MANDATORY_ATTRIBUTES +from improver.psychrometric_calculations.precipitable_water import ( + PrecipitableWater, +) from improver.psychrometric_calculations.psychrometric_calculations import ( HumidityMixingRatio, + flip_cube_in_place, get_pressure_points, - flip_cube_in_place - ) -from improver.psychrometric_calculations.precipitable_water import ( - PrecipitableWater, -) - -from improver.constants import EARTH_SURFACE_GRAVITY_ACCELERATION, WATER_DENSITY from improver.synthetic_data.set_up_test_cubes import set_up_variable_cube LOCAL_MANDATORY_ATTRIBUTES = { @@ -192,14 +189,14 @@ def make_pressure_cube(temp_cube: Cube) -> Cube: # ---------------------------------------- # 1. Extract the 1D pressure coordinate # ---------------------------------------- - p_coord = temp_cube.coord("pressure") # DimCoord - p_vals = p_coord.points # 1D array (nz,) + p_coord = temp_cube.coord("pressure") # DimCoord + p_vals = p_coord.points # 1D array (nz,) # ---------------------------------------- # 2. Broadcast pressure to match the cube grid # ---------------------------------------- - p_3d = p_vals[:, None, None] # (nz, 1, 1) + p_3d = p_vals[:, None, None] # (nz, 1, 1) p_3d = np.broadcast_to(p_3d, temp_cube.shape) # (nz, ny, nx) # ---------------------------------------- @@ -211,16 +208,19 @@ def make_pressure_cube(temp_cube: Cube) -> Cube: long_name="air_pressure unit test", units=p_coord.units, dim_coords_and_dims=[ - (p_coord, 0), # vertical dimension - (temp_cube.coord(axis='y'), 1), # y coordinate - (temp_cube.coord(axis='x'), 2), # x coordinate + (p_coord, 0), # vertical dimension + (temp_cube.coord(axis="y"), 1), # y coordinate + (temp_cube.coord(axis="x"), 2), # x coordinate ], - attributes = temp_cube.attributes, + attributes=temp_cube.attributes, ) return pressure_cube -def set_up_temperature_cube(shape: Tuple[int], temperature_value: float, vertical_levels: List[float]) -> Cube: + +def set_up_temperature_cube( + shape: Tuple[int], temperature_value: float, vertical_levels: List[float] +) -> Cube: """ Create a temperature on pressure cube. :param: shape: Shape of the temperature cube @@ -233,15 +233,18 @@ def set_up_temperature_cube(shape: Tuple[int], temperature_value: float, vertica np.full(shape, temperature_value, dtype=np.float32), "latlon", name="air_temperature", - x_grid_spacing = 1.0, - y_grid_spacing = 1.0, - vertical_levels = vertical_levels, - pressure = True + x_grid_spacing=1.0, + y_grid_spacing=1.0, + vertical_levels=vertical_levels, + pressure=True, ) add_attribute_dictionary(temperature) return temperature -def set_up_rel_humidity_cube(shape: Tuple[int], rel_humidity_value: float, vertical_levels: List[float]) -> Cube: + +def set_up_rel_humidity_cube( + shape: Tuple[int], rel_humidity_value: float, vertical_levels: List[float] +) -> Cube: """ Create a relative humidity on pressure cube. :param: shape: Shape of the relative humidity cube @@ -257,11 +260,12 @@ def set_up_rel_humidity_cube(shape: Tuple[int], rel_humidity_value: float, verti x_grid_spacing=1.0, y_grid_spacing=1.0, vertical_levels=vertical_levels, - pressure=True + pressure=True, ) add_attribute_dictionary(rel_humidity) return rel_humidity + def test_get_pressure_points() -> None: """ tests function "get_pressure_points" which is a support function @@ -270,28 +274,32 @@ def test_get_pressure_points() -> None: :return: None """ - temperature_value, pressure_value, rel_humidity_value, expected = (293, 100000, 0.1, 1.459832e-3) + temperature_value, rel_humidity_value = ( + 293, + 0.1, + ) # set up cubes vertical_levels = [100000.0, 50000.0, 100.0] - shape = (len(vertical_levels),3,3) + shape = (len(vertical_levels), 3, 3) temperature = set_up_temperature_cube(shape, temperature_value, vertical_levels) pressure = make_pressure_cube(temperature) rel_humidity = set_up_rel_humidity_cube(shape, rel_humidity_value, vertical_levels) - assert np.allclose( get_pressure_points(temperature), np.array(vertical_levels) ) - assert np.allclose( get_pressure_points(pressure), np.array(vertical_levels) ) - assert np.allclose( get_pressure_points(rel_humidity), np.array(vertical_levels) ) + assert np.allclose(get_pressure_points(temperature), np.array(vertical_levels)) + assert np.allclose(get_pressure_points(pressure), np.array(vertical_levels)) + assert np.allclose(get_pressure_points(rel_humidity), np.array(vertical_levels)) # check captialisation has no affect - rel_humidity.coord('pressure').rename("Pressure") + rel_humidity.coord("pressure").rename("Pressure") assert np.allclose(get_pressure_points(rel_humidity), np.array(vertical_levels)) # check null result when no "pressure" dimension - rel_humidity.coord('Pressure').rename("Pr3ssure") + rel_humidity.coord("Pressure").rename("Pr3ssure") assert np.allclose(get_pressure_points(rel_humidity), np.array([])) + def add_attribute_dictionary(cube: Cube) -> None: """ Adds attributes dictionary to cube attributes @@ -303,14 +311,15 @@ def add_attribute_dictionary(cube: Cube) -> None: """ # set up meta-data required by testing attributes_dictionary = { - "title":"unit test data", - "source":"unit test", - "institution":"somewhere", - "least_significant_digit": 4 + "title": "unit test data", + "source": "unit test", + "institution": "somewhere", + "least_significant_digit": 4, } for k, v in attributes_dictionary.items(): cube.attributes[k] = v + def test_mixing_ratio_without_pressure_parameter() -> None: """ the HumidityMixingRatio calculation will generate its own pressure cube @@ -331,12 +340,16 @@ def test_mixing_ratio_without_pressure_parameter() -> None: The improver calculation is then compared against a DIY calculation as a sanity check. """ - iris.FUTURE.save_split_attrs = True # to stop Iris warning - temperature_value, pressure_value, rel_humidity_value, expected = (293, 100000, 0.1, 1.459832e-3) + iris.FUTURE.save_split_attrs = True # to stop Iris warning + temperature_value, rel_humidity_value, expected = ( + 293, + 0.1, + 1.459832e-3, + ) # set up input cubes vertical_levels = [100000.0, 50000.0, 100.0] - shape = (len(vertical_levels),3,3) + shape = (len(vertical_levels), 3, 3) temperature = set_up_temperature_cube(shape, temperature_value, vertical_levels) pressure = make_pressure_cube(temperature) @@ -344,13 +357,13 @@ def test_mixing_ratio_without_pressure_parameter() -> None: # mixing ratio calculation with 3 parameters w3 = HumidityMixingRatio()([temperature, pressure, rel_humidity]) - metadata_ok(w3, temperature) # asserts in function call + metadata_ok(w3, temperature) # asserts in function call # check results on single layer are as expected where pressure is 100000 Pa assert np.isclose(w3.data[0], expected, atol=1e-7).all() # mixing ratio calculation with 2 parameters w2 = HumidityMixingRatio()([temperature, rel_humidity]) - metadata_ok(w2, temperature) # asserts in function call + metadata_ok(w2, temperature) # asserts in function call # check 2 parameter calculation gives same results as 3 parameter calculation assert np.isclose(w3.data, w2.data).all() @@ -362,12 +375,16 @@ def test_mixing_ratio_without_pressure_parameter() -> None: improver_tpw = np.sum(pw.data, axis=0) # perform DIY TPW calculation - delta = pressure.data[1:,:,:] - pressure.data[:-1,:,:] - mid_w = ( w3.data[1:,:,:] + w3.data[:-1,:,:] ) / 2.0 - integral_terms = (delta * mid_w) - unit_test_tpw_1 = -np.sum(integral_terms, axis=0) / (EARTH_SURFACE_GRAVITY_ACCELERATION * WATER_DENSITY) + delta = pressure.data[1:, :, :] - pressure.data[:-1, :, :] + mid_w = (w3.data[1:, :, :] + w3.data[:-1, :, :]) / 2.0 + integral_terms = delta * mid_w + unit_test_tpw_1 = -np.sum(integral_terms, axis=0) / ( + EARTH_SURFACE_GRAVITY_ACCELERATION * WATER_DENSITY + ) # numpy's trapezium rule integration - unit_test_tpw_2 = -np.trapezoid(w3.data, x=pressure.data, axis=0) / (EARTH_SURFACE_GRAVITY_ACCELERATION * WATER_DENSITY) + unit_test_tpw_2 = -np.trapezoid(w3.data, x=pressure.data, axis=0) / ( + EARTH_SURFACE_GRAVITY_ACCELERATION * WATER_DENSITY + ) # verify DIY TPW integrations produce same results # N.B. Improver uses np.trapezoid @@ -381,7 +398,7 @@ def test_mixing_ratio_without_pressure_parameter() -> None: assert np.isclose(improver_tpw, unit_test_tpw_1, rtol=2).all() -def generate_cube_bounds_for_axis(cube : Cube, axis: int) -> np.ndarray: +def generate_cube_bounds_for_axis(cube: Cube, axis: int) -> np.ndarray: """ automatically generates a suitable set of cube bounds for an axis support function to test "flip_cube_in_place" when the cube has bounds @@ -392,16 +409,17 @@ def generate_cube_bounds_for_axis(cube : Cube, axis: int) -> np.ndarray: """ points = cube.dim_coords[axis].points delta = points[1:] - points[:-1] - delta = np.concatenate((delta[:1], delta, delta[-1:] )) + delta = np.concatenate((delta[:1], delta, delta[-1:])) # N.B. delta can be negative so "low" and "high" are nominal - low_bounds = points - delta[:-1] / 2 - high_bounds = points + delta[1: ] / 2 + low_bounds = points - delta[:-1] / 2 + high_bounds = points + delta[1:] / 2 np_bounds = np.transpose(np.vstack((low_bounds, high_bounds))) # N.B. bounds are not necessarily [low,high] but will be rectified by Iris to follow ordering of points so can be [high, low] # i.e. the following line to give a genuine [low, high] ordering is not required and may be confusingly over-ridden in any case # np_bounds = np.sort(np_bounds, axis=-1) return np_bounds + @pytest.mark.parametrize( "bounds", ( @@ -430,9 +448,9 @@ def test_flip_cube_in_place(bounds: bool) -> None: iris.FUTURE.date_microseconds = True axis = 0 - temperature_value, pressure_value, rel_humidity_value, expected = (293, 100000, 0.1, 1.459832e-3) + temperature_value = (293,) vertical_levels = [100000.0, 50000.0, 100.0] - shape = (len(vertical_levels),3,3) + shape = (len(vertical_levels), 3, 3) temperature = set_up_temperature_cube(shape, temperature_value, vertical_levels) temperature.data[:] = np.random.uniform(273, 280, size=temperature.data.shape) if bounds: @@ -445,10 +463,14 @@ def test_flip_cube_in_place(bounds: bool) -> None: flip_cube_in_place(temperature) # check data flipped - np.allclose(np.flip(temperature_before_first_flip.data, axis=axis), temperature.data) + np.allclose( + np.flip(temperature_before_first_flip.data, axis=axis), temperature.data + ) # check coordinate points flipped - assert list(reversed(temperature.dim_coords[axis].points)) == list(temperature_before_first_flip.dim_coords[axis].points) + assert list(reversed(temperature.dim_coords[axis].points)) == list( + temperature_before_first_flip.dim_coords[axis].points + ) # check bounds flipped if bounds: @@ -457,9 +479,9 @@ def test_flip_cube_in_place(bounds: bool) -> None: b = temperature.dim_coords[axis].bounds # note that the bounds flip is unexpectedly complex # flipping on both axes is required. - a_f = np.flip( np.flip(a, axis=0), axis= 1) + a_f = np.flip(np.flip(a, axis=0), axis=1) - assert np.allclose( a_f, b ) + assert np.allclose(a_f, b) # second flip flip_cube_in_place(temperature) From 84659c80079c46a951e0a964fb573dc6b9175f38 Mon Sep 17 00:00:00 2001 From: mo-DavidJohnJohnston Date: Wed, 1 Apr 2026 23:36:13 +0100 Subject: [PATCH 3/9] add my name as an Improver developer so I can commit changes --- .mailmap | 1 + CONTRIBUTING.md | 1 + 2 files changed, 2 insertions(+) diff --git a/.mailmap b/.mailmap index 462db7e116..5b1d0522b2 100644 --- a/.mailmap +++ b/.mailmap @@ -16,6 +16,7 @@ Caroline Sandford <35029690+cgsandford@user Carwyn Pelley Chris Sampson <44228125+BelligerG@users.noreply.github.com> Daniel Mentiplay +David Johnston Eleanor Smith <40183561+ellesmith88@users.noreply.github.com> <40183561+ellesmith88@users.noreply.github.com> Fiona Rust Gavin Evans diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 51261b641f..2d60a4d42a 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -60,6 +60,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) From eb4067a893c45096db1d09964e3472112b1a5e85 Mon Sep 17 00:00:00 2001 From: mo-DavidJohnJohnston Date: Thu, 2 Apr 2026 09:22:02 +0100 Subject: [PATCH 4/9] white space fix! --- CONTRIBUTING.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 2d60a4d42a..5655537b90 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -60,7 +60,7 @@ below: - Timothy Hume (Bureau of Meteorology, Australia) - Katharine Hurst (Met Office, UK) - Simon Jackson (Met Office, UK) - - David Johnston (Met Office, UK) + - David Johnston (Met Office, UK) - Caroline Jones (Met Office, UK) - Peter Jordan (Met Office, UK) - Anzer Khan (Met Office, UK) From 194620ff829a17779ca94749a14eb98018934d1b Mon Sep 17 00:00:00 2001 From: mo-DavidJohnJohnston Date: Thu, 2 Apr 2026 16:08:04 +0100 Subject: [PATCH 5/9] use Iris "reverse" utility function instead of writing my own "flip" one. thanks to Stephen Moseley for telling me about this. --- .../psychrometric_calculations.py | 29 +----- .../test_HumidityMixingRatio.py | 95 ------------------- 2 files changed, 1 insertion(+), 123 deletions(-) diff --git a/improver/psychrometric_calculations/psychrometric_calculations.py b/improver/psychrometric_calculations/psychrometric_calculations.py index 10fd039332..71205324f3 100644 --- a/improver/psychrometric_calculations/psychrometric_calculations.py +++ b/improver/psychrometric_calculations/psychrometric_calculations.py @@ -474,33 +474,6 @@ def get_pressure_points(cube: Cube) -> np.ndarray: return np.array([]) -def flip_cube_in_place(cube: Cube, axis: int = 0) -> None: - """ - flip an Iris cube in-place along the specified axis - flips bounds if present - - :param cube: input cube to modify in-place - :param axis: axis along which to flip the cube - :return: None - - """ - # ---- Flip data in-place ---- - cube.data[:] = np.flip(cube.data, axis=axis) - - # ---- Find the dimension coord for the axis ---- - coord = cube.dim_coords[axis] - - # ---- Flip coord points in-place ---- - coord.points = coord.points[::-1] # 1D reversal - - # ---- Flip coord bounds if present ---- - if coord.bounds is not None: - # bounds flip is more complex than expected - flip on both axes - # as iris bounds are [low, high] or [high, low] - # depending on direction of points - coord.bounds = np.flip(np.flip(coord.bounds, axis=0), axis=1) - - class HumidityMixingRatio(BasePlugin): """Returns the humidity mass mixing ratio from temperature, pressure and relative humidity""" @@ -586,7 +559,7 @@ def generate_pressure_cube(self, temperature_cube) -> Cube: np.flip(pressure_points_for_pressure), pressure_points_for_temperature ) if flip_required: - flip_cube_in_place(pressure_cube, axis=0) + pressure_cube = iris.util.reverse(pressure_cube, "pressure") except iris.exceptions.ConcatenateError as error: raise RuntimeError( diff --git a/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py b/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py index 4345601970..1d6f260b3c 100644 --- a/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py +++ b/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py @@ -4,7 +4,6 @@ # See LICENSE in the root of the repository for full licensing details. """Tests for the HumidityMixingRatio plugin""" -from copy import deepcopy from typing import List, Tuple from unittest.mock import patch, sentinel @@ -22,7 +21,6 @@ ) from improver.psychrometric_calculations.psychrometric_calculations import ( HumidityMixingRatio, - flip_cube_in_place, get_pressure_points, ) from improver.synthetic_data.set_up_test_cubes import set_up_variable_cube @@ -398,99 +396,6 @@ def test_mixing_ratio_without_pressure_parameter() -> None: assert np.isclose(improver_tpw, unit_test_tpw_1, rtol=2).all() -def generate_cube_bounds_for_axis(cube: Cube, axis: int) -> np.ndarray: - """ - automatically generates a suitable set of cube bounds for an axis - support function to test "flip_cube_in_place" when the cube has bounds - - :param cube: Cube to generate bounds for - :param axis: axis to generate bounds for - :return: numpy array of bounds of shape (Nz,2) - """ - points = cube.dim_coords[axis].points - delta = points[1:] - points[:-1] - delta = np.concatenate((delta[:1], delta, delta[-1:])) - # N.B. delta can be negative so "low" and "high" are nominal - low_bounds = points - delta[:-1] / 2 - high_bounds = points + delta[1:] / 2 - np_bounds = np.transpose(np.vstack((low_bounds, high_bounds))) - # N.B. bounds are not necessarily [low,high] but will be rectified by Iris to follow ordering of points so can be [high, low] - # i.e. the following line to give a genuine [low, high] ordering is not required and may be confusingly over-ridden in any case - # np_bounds = np.sort(np_bounds, axis=-1) - return np_bounds - - -@pytest.mark.parametrize( - "bounds", - ( - (False), - (True), - ), -) -def test_flip_cube_in_place(bounds: bool) -> None: - """ - this tests support function "flip_cube_in_place" for implementation of - Improver's PrecipitableWater. It is used when a pressure cube - is found to be inadvertantly flipped in the Z-dimension. - - In this test, a cube is flipped once to check the data, and other components - have individually been flipped. - - Then the cube is flipped again to check we come back to the starting point - i.e. the operation is its own inverse. - - :param bounds: if True then the cube tested will have bounds - if False then the cube tested will not have bounds - - :return: None - """ - # to pacify an Iris warning - iris.FUTURE.date_microseconds = True - - axis = 0 - temperature_value = (293,) - vertical_levels = [100000.0, 50000.0, 100.0] - shape = (len(vertical_levels), 3, 3) - temperature = set_up_temperature_cube(shape, temperature_value, vertical_levels) - temperature.data[:] = np.random.uniform(273, 280, size=temperature.data.shape) - if bounds: - # set up pressure bounds for temperature on pressure cube - b = generate_cube_bounds_for_axis(temperature, axis) - temperature.dim_coords[axis].bounds = b - - # first flip - temperature_before_first_flip = deepcopy(temperature) - flip_cube_in_place(temperature) - - # check data flipped - np.allclose( - np.flip(temperature_before_first_flip.data, axis=axis), temperature.data - ) - - # check coordinate points flipped - assert list(reversed(temperature.dim_coords[axis].points)) == list( - temperature_before_first_flip.dim_coords[axis].points - ) - - # check bounds flipped - if bounds: - # performing an a, b comparison - a = temperature_before_first_flip.dim_coords[axis].bounds - b = temperature.dim_coords[axis].bounds - # note that the bounds flip is unexpectedly complex - # flipping on both axes is required. - a_f = np.flip(np.flip(a, axis=0), axis=1) - - assert np.allclose(a_f, b) - - # second flip - flip_cube_in_place(temperature) - - # cube equality is generally not a good idea - but the flip transform is simple - exactly_equal = temperature == temperature_before_first_flip - assert exactly_equal - - def test_height_levels(): """Check that the plugin works with height level data""" From 63ee61620c86f5216322f8b2c0bc105def307a97 Mon Sep 17 00:00:00 2001 From: mo-DavidJohnJohnston Date: Wed, 8 Apr 2026 09:08:49 +0100 Subject: [PATCH 6/9] use Iris "reverse" utility function instead of writing my own "flip" one. thanks to Stephen Moseley for telling me about this. --- .../psychrometric_calculations.py | 29 ++++++-- .../test_HumidityMixingRatio.py | 68 +++++++++++-------- 2 files changed, 61 insertions(+), 36 deletions(-) diff --git a/improver/psychrometric_calculations/psychrometric_calculations.py b/improver/psychrometric_calculations/psychrometric_calculations.py index 71205324f3..c11cfd1c22 100644 --- a/improver/psychrometric_calculations/psychrometric_calculations.py +++ b/improver/psychrometric_calculations/psychrometric_calculations.py @@ -465,11 +465,17 @@ def qsat_differential(qs, t, q, p): def get_pressure_points(cube: Cube) -> np.ndarray: """ Get the pressure points from a _on_pressure_levels cube. - :param cube: input cube - :return: pressure points + If no pressure coordinate is found, an empty array is returned. + + Args: + cube: input cube + + Returns: + array of presssure points """ + for coord in cube.dim_coords: - if "pressure" in coord.name().lower(): + if "pressure" == coord.name(): return coord.points return np.array([]) @@ -549,13 +555,22 @@ def generate_pressure_cube(self, temperature_cube) -> Cube: try: pressure_cube = expanded_pressure_list.concatenate_cube() """ - the Iris concatenate_cube function can reverse the list order when forming the cube - so the pressure cube is flipped vertically compared to the temperature cube - check if this is the case and then re-flip + The Iris concatenate_cube function can reverse the list order when forming the cube + so the pressure cube is flipped vertically compared to the temperature cube. + Check if this is the case and then re-flip. """ pressure_points_for_pressure = get_pressure_points(pressure_cube) pressure_points_for_temperature = get_pressure_points(temperature_cube) - flip_required = np.array_equal( + """ + To stop flipping occurring if there are no pressure coordinates, + check that one array is not empty via the size test, also testing + if a flip is worthwhile (i.e. at least 2 elements). + + However, the input (temperature_on_pressure cube) argument is + preconditioned to have pressure points, so this test should not be + required so it also serves as defensive programming. + """ + flip_required = (pressure_points_for_pressure.size > 1) and np.allclose( np.flip(pressure_points_for_pressure), pressure_points_for_temperature ) if flip_required: diff --git a/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py b/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py index 1d6f260b3c..0eaebc3d51 100644 --- a/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py +++ b/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py @@ -175,12 +175,13 @@ def test_zero_humidity( def make_pressure_cube(temp_cube: Cube) -> Cube: - """ - Create a 3D pressure cube from a temperature_on_pressure_levels cube. + """Create a 3D pressure cube from a temperature_on_pressure_levels cube. The resulting cube has shape (levels, y, x). - :param temp_cube: input temperature on pressure cube - :return: a 3D pressure cube + Args: + temp_cube: input temperature on pressure cube + + Returns: a 3D pressure cube """ @@ -219,14 +220,17 @@ def make_pressure_cube(temp_cube: Cube) -> Cube: def set_up_temperature_cube( shape: Tuple[int], temperature_value: float, vertical_levels: List[float] ) -> Cube: - """ - Create a temperature on pressure cube. - :param: shape: Shape of the temperature cube - :param: temperature_value: temperature value - :param: vertical_levels: List of vertical levels - :return: a temperature cube + """Create a temperature on pressure cube. + + Args: + shape: Shape of the temperature cube + temperature_value: temperature value + vertical_levels: List of vertical levels + + Returns: a temperature cube """ + temperature = set_up_variable_cube( np.full(shape, temperature_value, dtype=np.float32), "latlon", @@ -243,13 +247,17 @@ def set_up_temperature_cube( def set_up_rel_humidity_cube( shape: Tuple[int], rel_humidity_value: float, vertical_levels: List[float] ) -> Cube: + """Create a relative humidity on pressure cube. + + Args: + shape: shape of the relative humidity cube + rel_humidity_value: relative humidity value + vertical_levels: list of vertical levels + + Returns: a relative humidity cube + """ - Create a relative humidity on pressure cube. - :param: shape: Shape of the relative humidity cube - :param: rel_humidity_value: relative humidity value - :param: vertical_levels: List of vertical levels - :return: a relative humidity cube - """ + rel_humidity = set_up_variable_cube( np.full(shape, rel_humidity_value, dtype=np.float32), "latlon", @@ -265,10 +273,9 @@ def set_up_rel_humidity_cube( def test_get_pressure_points() -> None: - """ - tests function "get_pressure_points" which is a support function + """Tests for function "get_pressure_points" which is a support function written to check if a pressure cube has been inadvertantly flipped - within the Improver implementation of PrecipitableWater + within the Improver implementation of PrecipitableWater. :return: None """ @@ -289,24 +296,28 @@ def test_get_pressure_points() -> None: assert np.allclose(get_pressure_points(pressure), np.array(vertical_levels)) assert np.allclose(get_pressure_points(rel_humidity), np.array(vertical_levels)) - # check captialisation has no affect + # check captialisation has an effect (i.e. the function gives a null result) + # the meta-data should be CF compliant and not be capitalised in any way rel_humidity.coord("pressure").rename("Pressure") - assert np.allclose(get_pressure_points(rel_humidity), np.array(vertical_levels)) + assert np.allclose(get_pressure_points(rel_humidity), np.array([])) # check null result when no "pressure" dimension - rel_humidity.coord("Pressure").rename("Pr3ssure") + rel_humidity.coord("Pressure").rename("pr3ssure") assert np.allclose(get_pressure_points(rel_humidity), np.array([])) def add_attribute_dictionary(cube: Cube) -> None: - """ - Adds attributes dictionary to cube attributes + """Adds attributes dictionary to cube attributes to allow pre-existing checking function "metadata_ok" to be used. - :param cube: Cube to add attributes to - :return: None + Args: + cube: Cube to add attributes to + + Returns: None + """ + # set up meta-data required by testing attributes_dictionary = { "title": "unit test data", @@ -319,8 +330,7 @@ def add_attribute_dictionary(cube: Cube) -> None: def test_mixing_ratio_without_pressure_parameter() -> None: - """ - the HumidityMixingRatio calculation will generate its own pressure cube + """The HumidityMixingRatio calculation will generate its own pressure cube if one is not supplied. This unit tests verifies that the results are the same with/without an explicit pressure parameter. @@ -330,7 +340,7 @@ def test_mixing_ratio_without_pressure_parameter() -> None: The reason was that HumidityMixingRadio generated a pressure cube that was wrongly flipped veritically. This unit test re-creates the failing scenario to test and - exercise the bug fix + exercise the bug fix. The unit test then does a very simple total precipitable water calculation ensuring the output from HumidityMixingRadio is suitable. From 21ad0bb2e83b5f14a792b05806092447870141cf Mon Sep 17 00:00:00 2001 From: mo-DavidJohnJohnston Date: Mon, 13 Apr 2026 23:14:27 +0100 Subject: [PATCH 7/9] changes made in response to second review --- .../psychrometric_calculations.py | 36 +++++--- .../test_HumidityMixingRatio.py | 85 +++---------------- 2 files changed, 36 insertions(+), 85 deletions(-) diff --git a/improver/psychrometric_calculations/psychrometric_calculations.py b/improver/psychrometric_calculations/psychrometric_calculations.py index c11cfd1c22..018230561b 100644 --- a/improver/psychrometric_calculations/psychrometric_calculations.py +++ b/improver/psychrometric_calculations/psychrometric_calculations.py @@ -537,6 +537,24 @@ def generate_pressure_cube(self, temperature_cube) -> Cube: See https://scitools-iris.readthedocs.io/en/stable/further_topics/controlling_merge.html for more information + + Note that the underlying Iris function concatenate_cube, which constructs the pressure + cube from a list of component layers, can have the unexpected behaviour of effectively + reversing the supplied list order when forming the cube. The result is that the pressure + cube is flipped vertically compared to the temperature cube. The function iris.util.reverse + is used to reflip the cube back to its expected configuration. + + As security, the presence of the unexpected flip is checked for before the reflip is executed. + + To stop flipping occurring if there are no pressure coordinates, + check that one coordinate array is not empty via a size test, which + also tests if a flip is worthwhile (i.e. at least 2 elements). + + However, the input (temperature_on_pressure cube) argument is + preconditioned to have pressure points, so the size test should not be + required as defensive programming but it does serve a useful purpose + in the single layer case to avoid unnecessary processing. + """ coord_list = [coord.name() for coord in temperature_cube.coords()] @@ -553,27 +571,17 @@ def generate_pressure_cube(self, temperature_cube) -> Cube: ) try: + # the cube constructed here may be unexpectedly flipped vertically pressure_cube = expanded_pressure_list.concatenate_cube() - """ - The Iris concatenate_cube function can reverse the list order when forming the cube - so the pressure cube is flipped vertically compared to the temperature cube. - Check if this is the case and then re-flip. - """ + pressure_points_for_pressure = get_pressure_points(pressure_cube) pressure_points_for_temperature = get_pressure_points(temperature_cube) - """ - To stop flipping occurring if there are no pressure coordinates, - check that one array is not empty via the size test, also testing - if a flip is worthwhile (i.e. at least 2 elements). - - However, the input (temperature_on_pressure cube) argument is - preconditioned to have pressure points, so this test should not be - required so it also serves as defensive programming. - """ + flip_required = (pressure_points_for_pressure.size > 1) and np.allclose( np.flip(pressure_points_for_pressure), pressure_points_for_temperature ) if flip_required: + # reflip the cube to its correct configuration pressure_cube = iris.util.reverse(pressure_cube, "pressure") except iris.exceptions.ConcatenateError as error: diff --git a/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py b/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py index 0eaebc3d51..35b85e3184 100644 --- a/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py +++ b/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py @@ -7,18 +7,12 @@ from typing import List, Tuple from unittest.mock import patch, sentinel -import iris -import iris.cube as icube import numpy as np import pytest from iris.coords import AncillaryVariable from iris.cube import Cube -from improver.constants import EARTH_SURFACE_GRAVITY_ACCELERATION, WATER_DENSITY from improver.metadata.constants.attributes import MANDATORY_ATTRIBUTES -from improver.psychrometric_calculations.precipitable_water import ( - PrecipitableWater, -) from improver.psychrometric_calculations.psychrometric_calculations import ( HumidityMixingRatio, get_pressure_points, @@ -185,34 +179,13 @@ def make_pressure_cube(temp_cube: Cube) -> Cube: """ - # ---------------------------------------- - # 1. Extract the 1D pressure coordinate - # ---------------------------------------- p_coord = temp_cube.coord("pressure") # DimCoord p_vals = p_coord.points # 1D array (nz,) + p_3d = np.broadcast_to(p_vals[:, None, None], temp_cube.shape) # (nz, ny, nx) - # ---------------------------------------- - # 2. Broadcast pressure to match the cube grid - # ---------------------------------------- - - p_3d = p_vals[:, None, None] # (nz, 1, 1) - p_3d = np.broadcast_to(p_3d, temp_cube.shape) # (nz, ny, nx) - - # ---------------------------------------- - # 3. Build a new pressure cube - # ---------------------------------------- - pressure_cube = icube.Cube( - p_3d, - standard_name="air_pressure", - long_name="air_pressure unit test", - units=p_coord.units, - dim_coords_and_dims=[ - (p_coord, 0), # vertical dimension - (temp_cube.coord(axis="y"), 1), # y coordinate - (temp_cube.coord(axis="x"), 2), # x coordinate - ], - attributes=temp_cube.attributes, - ) + pressure_cube = temp_cube.copy(p_3d) + pressure_cube.rename("air_pressure") + pressure_cube.units = p_coord.units return pressure_cube @@ -235,8 +208,6 @@ def set_up_temperature_cube( np.full(shape, temperature_value, dtype=np.float32), "latlon", name="air_temperature", - x_grid_spacing=1.0, - y_grid_spacing=1.0, vertical_levels=vertical_levels, pressure=True, ) @@ -263,8 +234,6 @@ def set_up_rel_humidity_cube( "latlon", name="relative_humidity", units="1", - x_grid_spacing=1.0, - y_grid_spacing=1.0, vertical_levels=vertical_levels, pressure=True, ) @@ -348,17 +317,20 @@ def test_mixing_ratio_without_pressure_parameter() -> None: The improver calculation is then compared against a DIY calculation as a sanity check. """ - iris.FUTURE.save_split_attrs = True # to stop Iris warning - temperature_value, rel_humidity_value, expected = ( + + temperature_value, rel_humidity_value = ( 293, 0.1, - 1.459832e-3, ) - # set up input cubes vertical_levels = [100000.0, 50000.0, 100.0] - shape = (len(vertical_levels), 3, 3) + shape = (len(vertical_levels), 4, 5) + # create expected result as a numpy array + w_data_col_expected = np.array([0.0014598317, 0.0029387323, 0.1]) + w_data_expected = np.broadcast_to(w_data_col_expected[:, None, None], shape) + + # set up input cubes temperature = set_up_temperature_cube(shape, temperature_value, vertical_levels) pressure = make_pressure_cube(temperature) rel_humidity = set_up_rel_humidity_cube(shape, rel_humidity_value, vertical_levels) @@ -367,43 +339,14 @@ def test_mixing_ratio_without_pressure_parameter() -> None: w3 = HumidityMixingRatio()([temperature, pressure, rel_humidity]) metadata_ok(w3, temperature) # asserts in function call # check results on single layer are as expected where pressure is 100000 Pa - assert np.isclose(w3.data[0], expected, atol=1e-7).all() + assert np.isclose(w3.data, w_data_expected, atol=1e-7).all() # mixing ratio calculation with 2 parameters w2 = HumidityMixingRatio()([temperature, rel_humidity]) metadata_ok(w2, temperature) # asserts in function call # check 2 parameter calculation gives same results as 3 parameter calculation - assert np.isclose(w3.data, w2.data).all() - - # use w3 to calculate precipitable water - pw = PrecipitableWater().process(w3) - - # perform integration step (summing water in vertical atmosphere column) for Improver - improver_tpw = np.sum(pw.data, axis=0) - - # perform DIY TPW calculation - delta = pressure.data[1:, :, :] - pressure.data[:-1, :, :] - mid_w = (w3.data[1:, :, :] + w3.data[:-1, :, :]) / 2.0 - integral_terms = delta * mid_w - unit_test_tpw_1 = -np.sum(integral_terms, axis=0) / ( - EARTH_SURFACE_GRAVITY_ACCELERATION * WATER_DENSITY - ) - # numpy's trapezium rule integration - unit_test_tpw_2 = -np.trapezoid(w3.data, x=pressure.data, axis=0) / ( - EARTH_SURFACE_GRAVITY_ACCELERATION * WATER_DENSITY - ) - - # verify DIY TPW integrations produce same results - # N.B. Improver uses np.trapezoid - assert np.isclose(unit_test_tpw_1.data, unit_test_tpw_2.data).all() - - # note for such a small cube the DIY and Improver calculations - # for TPW are somewhat different. - # one presumes they will converge for cubes with more cells - # give a 200% latitude for testing - - assert np.isclose(improver_tpw, unit_test_tpw_1, rtol=2).all() + assert np.isclose(w2.data, w_data_expected).all() def test_height_levels(): From fb0448a6f5ebd22b4604d20a86526e0059e25d05 Mon Sep 17 00:00:00 2001 From: mo-DavidJohnJohnston Date: Tue, 14 Apr 2026 11:21:06 +0100 Subject: [PATCH 8/9] changes made in response to second second review. --- .../psychrometric_calculations.py | 2 +- .../test_HumidityMixingRatio.py | 22 ++++--------------- 2 files changed, 5 insertions(+), 19 deletions(-) diff --git a/improver/psychrometric_calculations/psychrometric_calculations.py b/improver/psychrometric_calculations/psychrometric_calculations.py index 018230561b..e20ccb7a46 100644 --- a/improver/psychrometric_calculations/psychrometric_calculations.py +++ b/improver/psychrometric_calculations/psychrometric_calculations.py @@ -527,7 +527,7 @@ def generate_pressure_cube(self, temperature_cube) -> Cube: If there is a pressure coordinate in the temperature and relative humidity cubes, and no pressure cube has been provided (as is the case for calculating virtual temperature - on pressure levels, for example) the pressure cube is generated from the from the + on pressure levels, for example) the pressure cube is generated from the pressure coordinate on the temperature cube. The temperature cube has a status flag that indicates where the data were derived by StaGE for data points that fell below the model orography, the flag meaning is above_surface_pressure below_surface_pressure. diff --git a/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py b/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py index 35b85e3184..5d776fc70e 100644 --- a/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py +++ b/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py @@ -179,9 +179,9 @@ def make_pressure_cube(temp_cube: Cube) -> Cube: """ - p_coord = temp_cube.coord("pressure") # DimCoord - p_vals = p_coord.points # 1D array (nz,) - p_3d = np.broadcast_to(p_vals[:, None, None], temp_cube.shape) # (nz, ny, nx) + p_coord = temp_cube.coord("pressure") + p_vals = p_coord.points + p_3d = np.broadcast_to(p_vals[:, None, None], temp_cube.shape) pressure_cube = temp_cube.copy(p_3d) pressure_cube.rename("air_pressure") @@ -300,22 +300,8 @@ def add_attribute_dictionary(cube: Cube) -> None: def test_mixing_ratio_without_pressure_parameter() -> None: """The HumidityMixingRatio calculation will generate its own pressure cube - if one is not supplied. This unit tests verifies that the results are the + if one is not supplied. This unit test verifies that the results are the same with/without an explicit pressure parameter. - - This ticket reports values for the total precipitable water (TPW) being far too high. - - https://metoffice.atlassian.net/browse/EPPT-3209 - - The reason was that HumidityMixingRadio generated a pressure cube that was wrongly - flipped veritically. This unit test re-creates the failing scenario to test and - exercise the bug fix. - - The unit test then does a very simple total precipitable water calculation - ensuring the output from HumidityMixingRadio is suitable. - - The improver calculation is then compared against a DIY calculation as a sanity check. - """ temperature_value, rel_humidity_value = ( From c772c09a36b310a7ca415fc8e02ffac2c829a703 Mon Sep 17 00:00:00 2001 From: robertplatt-mo Date: Tue, 14 Apr 2026 11:38:04 +0100 Subject: [PATCH 9/9] Apply suggestion from @robertplatt-mo Bring return docstring in line with Improver style --- .../psychrometric_calculations/test_HumidityMixingRatio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py b/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py index 5d776fc70e..05cfe0091c 100644 --- a/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py +++ b/improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py @@ -246,7 +246,7 @@ def test_get_pressure_points() -> None: written to check if a pressure cube has been inadvertantly flipped within the Improver implementation of PrecipitableWater. - :return: None + Returns: None """ temperature_value, rel_humidity_value = ( 293,