EPPT_3209_total_precipitable_water_assurance to master

.mailmap

@@ -16,6 +16,7 @@ Caroline Sandford <caroline.sandford@metoffice.gov.uk> <35029690+cgsandford@user
 Carwyn Pelley <carwyn.pelley@metoffice.gov.uk> <carwyn.pelley@metoffice.gov.uk>
 Chris Sampson <christopher.sampson@metoffice.gov.uk> <44228125+BelligerG@users.noreply.github.com>
 Daniel Mentiplay <dmentipl@users.noreply.github.com> <dmentipl@users.noreply.github.com>
+David Johnston <david.johnston@metoffice.gov.uk> <david.johnston@metoffice.gov.uk>
 Eleanor Smith <40183561+ellesmith88@users.noreply.github.com> <40183561+ellesmith88@users.noreply.github.com>
 Fiona Rust <fiona.rust@metoffice.gov.uk> <fiona.rust@metoffice.gov.uk>
 Gavin Evans <gavin.evans@metoffice.gov.uk> <gavin.evans@metoffice.gov.uk>

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)

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
@@ -461,6 +462,24 @@ def qsat_differential(qs, t, q, p):
     return temperature, humidity
 
 
+def get_pressure_points(cube: Cube) -> np.ndarray:
+    """
+    Get the pressure points from a <diagnostic>_on_pressure_levels cube.
+    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" == coord.name():
+            return coord.points
+    return np.array([])
+
+
 class HumidityMixingRatio(BasePlugin):
     """Returns the humidity mass mixing ratio from temperature, pressure and relative humidity"""
 
@@ -535,6 +554,28 @@ 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)
+            """
+            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:
+                pressure_cube = iris.util.reverse(pressure_cube, "pressure")
+
         except iris.exceptions.ConcatenateError as error:
             raise RuntimeError(
                 "Unable to concatenate pressure cubelist with input ",

improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py

@@ -4,16 +4,24 @@
 # See LICENSE in the root of the repository for full licensing details.
 """Tests for the HumidityMixingRatio plugin"""
 
+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,
 )
 from improver.synthetic_data.set_up_test_cubes import set_up_variable_cube
 
@@ -166,6 +174,238 @@ 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).
+
+    Args:
+        temp_cube: input temperature on pressure cube
+
+    Returns: 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.
+
+    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",
+        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.
+
+    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
+
+    """
+
+    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 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.
+
+    :return: None
+    """
+    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)
+
+    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 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([]))
+
+    # 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.
+
+    Args:
+        cube:  Cube to add attributes to
+
+    Returns: 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, 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)
+
+    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 test_height_levels():
     """Check that the plugin works with height level data"""