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"""
 
@@ -508,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.
@@ -518,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()]
@@ -534,7 +571,19 @@ 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()
+
+            pressure_points_for_pressure = get_pressure_points(pressure_cube)
+            pressure_points_for_temperature = get_pressure_points(temperature_cube)
+
+            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:
             raise RuntimeError(
                 "Unable to concatenate pressure cubelist with input ",

improver_tests/psychrometric_calculations/test_HumidityMixingRatio.py

@@ -4,6 +4,7 @@
 # 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 numpy as np
@@ -14,6 +15,7 @@
 from improver.metadata.constants.attributes import MANDATORY_ATTRIBUTES
 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 +168,173 @@ 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
+
+    """
+
+    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")
+    pressure_cube.units = p_coord.units
+
+    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",
+        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",
+        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.
+
+    Returns: 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 test verifies that the results are the
+    same with/without an explicit pressure parameter.
+    """
+
+    temperature_value, rel_humidity_value = (
+        293,
+        0.1,
+    )
+
+    vertical_levels = [100000.0, 50000.0, 100.0]
+    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)
+
+    # 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, 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(w2.data, w_data_expected).all()
+
+
 def test_height_levels():
     """Check that the plugin works with height level data"""