Height of max (#2020)

* Plugin for calculating the height of the maximum vertical velocity once the maximum vertical velocity has been calculated.

* Minor code updates.

* Deleting or adding required spaces.

* Updating some comments.

* Adds masked_add to cube combiner

* Remove unnecessary import

* formatting

* Updating and improving unit tests, and adding cli.

* Isort changes.

* Formatting changes.

* Sorting out trailing whitespaces.

* Sorting out formatting.

* flake8 change.

* Further changes to get acceptance tests to work.

* Changes to get acceptance test data to run.

* add docstrings and simplifications

* formatting

* update checksums and correct acceptance test bug

* Update default name and change low_or_high to boolean

* updates acceptance tests

* formatting

---------

Co-authored-by: Kathryn.Howard <kathryn.howard@metoffice.gov.uk>
Co-authored-by: Marcus Spelman <marcus.spelman@metoffice.gov.uk>

improver/cli/height_of_max_vertical_velocity.py

@@ -0,0 +1,46 @@
+#!/usr/bin/env python
+# (C) Crown copyright, Met Office. All rights reserved.
+#
+# This file is part of IMPROVER and is released under a BSD 3-Clause license.
+# See LICENSE in the root of the repository for full licensing details.
+"""Script to calculate the height at which the maximum vertical velocity
+occurs. """
+
+from improver import cli
+
+
+@cli.clizefy
+@cli.with_output
+def process(
+    cube: cli.inputcube,
+    max_cube: cli.inputcube,
+    *,
+    find_lowest: bool = True,
+    new_name: str = "height_of_maximum_upward_air_velocity",
+):
+    """Calculates the height level at which the maximum vertical velocity occurs for each
+    grid point. It requires an input cube of vertical velocity and a cube with the maximum
+    vertical velocity values at each grid point. For this case we are looking for the
+    lowest height at which the maximum occurs.
+
+    Args:
+        cube (iris.cube.Cube):
+            A cube containing vertical velocity.
+        max_cube (iris.cube.Cube):
+            A cube containing the maximum values of vertical velocity over all heights.
+        find_lowest (bool):
+            If true then the lowest maximum height will be found (for cases where
+            there are two heights with the maximum vertical velocity.) Otherwise the highest
+            height will be found.
+        new_name (str):
+            The new name to be assigned to the output cube. In this case it will become
+            height_of_maximum_upward_air_velocity. If unspecified the name of the original
+            cube is used.
+
+    Returns:
+        A cube of heights at which the maximum value occurs.
+    """
+
+    from improver.utilities.cube_manipulation import height_of_maximum
+
+    return height_of_maximum(cube, max_cube, find_lowest, new_name)

improver/utilities/cube_manipulation.py

@@ -774,3 +774,53 @@ def maximum_in_height(
         max_cube.rename(new_name)
 
     return max_cube
+
+
+def height_of_maximum(
+    cube: Cube, max_cube: Cube, find_lowest: bool = True, new_name: str = None,
+) -> Cube:
+    """Calculates the height level at which the maximum value has been calculated. This
+    takes in a cube with values at different heights, and also a cube with the maximum
+    of these heights. It compares these (default is to start at the lowest height and
+    work down through the height levels), and then outputs the height it reaches the
+    maximum value.
+
+    Args:
+        cube:
+            A cube with a height coordinate.
+        max_cube:
+            A cube of the maximum value over the height coordinate.
+        find_lowest:
+            If true then the lowest maximum height will be found (for cases where
+            there are two heights with the maximum vertical velocity.) Otherwise the highest
+            height will be found.
+        new_name:
+            The new name to be assigned to the output cube. If unspecified the name of the
+            original cube is used.
+    Returns:
+        A cube of heights at which the maximum values occur.
+
+    Raises:
+        ValueError:
+            If the cube has only 1 height level or if an input other than high or low is
+            tried for the high_or_low value.
+    """
+    height_of_max = max_cube.copy()
+    height_range = range(len(cube.coord("height").points))
+    if len(cube.coord("height").points) == 1:
+        raise ValueError("More than 1 height level is required.")
+    if find_lowest:
+        height_points = height_range
+    else:
+        height_points = reversed(height_range)
+
+    for height in height_points:
+        height_of_max.data = np.where(
+            cube[height].data == max_cube.data,
+            cube[height].coord("height").points[0],
+            height_of_max.data,
+        )
+    if new_name:
+        height_of_max.rename(new_name)
+    height_of_max.units = cube.coord("height").units
+    return height_of_max

improver_tests/acceptance/SHA256SUMS

@@ -539,6 +539,9 @@ caa759d708afc535223bcf35924e98962675f6e3f690fe1f82c5de5ba08302b4  ./hail-size/te
 c4451e5a94a946215e5a3e09787b0d25c215aa4511110f46b6015c1c5aab13c9  ./hail-size/wet_bulb_freezing_altitude.nc
 0af18a52713334627696679bb369bb00332e5cb8f8a1a82ca9a2a7612071e6d3  ./hail-size/with_id_attr/kgo.nc
 76d84b674d8c0c9bed8bd27fad6697be7559c9ebe1c13296363717cbcd888add  ./hail-size/without_id_attr/kgo.nc
+e4ad2774923662fad733e3d95730416993abd94486aa9e032256e9d60d4b1bc0  ./height-of-max-vertical-velocity/kgo.nc
+90ac17c415ba2b0249de3f304bf2f511a9a294710e0577dac9231b6ab822660d  ./height-of-max-vertical-velocity/max_vertical_velocity.nc
+929f98fa947ca8b635d76f6d3509b368fe7780019af76172dddbae4fef21822d  ./height-of-max-vertical-velocity/vertical_velocity_on_height_levels.nc
 cebc2ccbe6c8e33b3e39d65f07d189172133a3fc6c22e3567c61572e14750cef  ./integrate-time-bounds/basic/kgo.nc
 b8ae3b9d3db05fc0c8479bb707465aeaf140202ab4477d025f5c793e2d287dc8  ./integrate-time-bounds/basic/kgo_renamed.nc
 5aaa03199faf9df5fda699936b33df862b071b3790b04791fef31d8cc0fd074a  ./integrate-time-bounds/basic/lightning_frequency.nc

improver_tests/acceptance/test_height_of_max_vertical_velocity.py

@@ -0,0 +1,31 @@
+# (C) Crown copyright, Met Office. All rights reserved.
+#
+# This file is part of IMPROVER and is released under a BSD 3-Clause license.
+# See LICENSE in the root of the repository for full licensing details.
+"""Tests the height_of_max_vertical_velocity cli"""
+
+import pytest
+
+from . import acceptance as acc
+
+pytestmark = [pytest.mark.acc, acc.skip_if_kgo_missing]
+CLI = acc.cli_name_with_dashes(__file__)
+run_cli = acc.run_cli(CLI)
+
+
+def test_basic(tmp_path):
+    """Test height_of_max_vertical_velocity computation"""
+    kgo_dir = acc.kgo_root() / "height-of-max-vertical-velocity"
+    input_file1 = kgo_dir / "vertical_velocity_on_height_levels.nc"
+    input_file2 = kgo_dir / "max_vertical_velocity.nc"
+    output_path = tmp_path / "output.nc"
+    args = [
+        input_file1,
+        input_file2,
+        "--output",
+        f"{output_path}",
+    ]
+
+    kgo_path = kgo_dir / "kgo.nc"
+    run_cli(args)
+    acc.compare(output_path, kgo_path)

improver_tests/utilities/cube_manipulation/test_height_of_maximum.py

@@ -0,0 +1,80 @@
+# (C) Crown copyright, Met Office. All rights reserved.
+#
+# This file is part of IMPROVER and is released under a BSD 3-Clause license.
+# See LICENSE in the root of the repository for full licensing details.
+"""
+Unit tests for the function "cube_manipulation.height_of_maximum".
+"""
+
+import numpy as np
+import pytest
+from iris.cube import Cube
+from numpy.testing import assert_allclose
+
+from improver.synthetic_data.set_up_test_cubes import set_up_variable_cube
+from improver.utilities.cube_manipulation import height_of_maximum
+
+
+@pytest.fixture(name="input_cube")
+def input_cube() -> Cube:
+    """Test cube of vertical velocity  on height levels"""
+    data = np.array(
+        [[[2, 4, 9], [3, 4, 8]], [[5, 3, 3], [4, 2, 7]], [[9, 5, 1], [2, 5, 8]]]
+    )
+    cube = set_up_variable_cube(
+        data=data, name="vertical_velocity", height_levels=[5, 75, 300]
+    )
+    return cube
+
+
+@pytest.fixture(name="max_cube")
+def max_cube() -> Cube:
+    """Test cube of maximum vertical velocities over the height levels"""
+    data = np.array([[9, 5, 9], [4, 5, 8]])
+    cube = set_up_variable_cube(data=data, name="vertical_velocity", height_levels=[1])
+    return cube
+
+
+@pytest.fixture(name="high_cube")
+def high_cube() -> Cube:
+    """Test cube when we want the highest maximum"""
+    data_high = np.array([[300, 300, 5], [75, 300, 300]])
+    cube = set_up_variable_cube(
+        data=data_high, name="vertical_velocity", height_levels=[1]
+    )
+    return cube
+
+
+@pytest.fixture(name="low_cube")
+def low_cube() -> Cube:
+    """Test cube when we want the lowest maximum"""
+    data_low = np.array([[300, 300, 5], [75, 300, 5]])
+    cube = set_up_variable_cube(
+        data=data_low, name="vertical_velocity", height_levels=[1]
+    )
+    return cube
+
+
+@pytest.mark.parametrize("new_name", [None, "height_of_maximum"])
+@pytest.mark.parametrize("find_lowest", ["True", "False"])
+def test_basic(input_cube, max_cube, new_name, find_lowest, high_cube, low_cube):
+    """Tests that the name of the cube will be correctly updated. Test that
+    if find_lowest is true the lowest maximum height will be found"""
+
+    expected_name = new_name if new_name else input_cube.name()
+    expected_cube = high_cube if find_lowest else low_cube
+
+    output_cube = height_of_maximum(
+        input_cube, max_cube, new_name=new_name, find_lowest=find_lowest
+    )
+
+    assert expected_name == output_cube.name()
+    assert output_cube.units == "m"
+    assert_allclose(output_cube.data, expected_cube.data)
+
+
+def test_one_height(input_cube):
+    one_height = input_cube[0]
+    msg = "More than 1 height level is required."
+    with pytest.raises(ValueError, match=msg):
+        height_of_maximum(one_height, one_height)