From 009831c661bff6cb981bb2e4278b4434c7bd5d8a Mon Sep 17 00:00:00 2001 From: mspelman07 <99179165+mspelman07@users.noreply.github.com> Date: Tue, 16 Jul 2024 15:17:54 +0100 Subject: [PATCH 1/7] Make difference module handle circular cubes (#2016) * Resolved merge conflict * Added another failing test for difference module. Differences calculated correctly but need to add extra dimension for the wrap-around distance. * Resolved merge conflicts * Auto-formatting * Import Bugfix * Ran isort against spatial.py * Added test for appropriate error for unhandled cubes and made code match it. * iSort on spatial.py again * Another attempt at autoformatting * Tidied up * Removed confusing variable name * Autoformatting * Tried manually sorting import order as isort still failing in github actions despite running locally and pushing * And again * And again * Linting suggestion * Review Actions * Ran isort to fix gha * Added test for cubes with flipped axes (lonitude increasing along columns not rows and vice versa for latitude). Fixed code to handle this case. * black, isort, flake8 * Ran isort in isolation * Fixed test failures resulting from breaking changes made by upstream PR * Black * updates as per review comments and adds test for 0-360 degree case * update comment * formatting --------- Co-authored-by: Peter Jordan Co-authored-by: MO-PeterJordan <52462411+MO-PeterJordan@users.noreply.github.com> Co-authored-by: Marcus Spelman --- improver/utilities/spatial.py | 90 ++++++++++++--- ...st_DifferenceBetweenAdjacentGridSquares.py | 105 +++++++++++++++--- 2 files changed, 168 insertions(+), 27 deletions(-) diff --git a/improver/utilities/spatial.py b/improver/utilities/spatial.py index 72ec144b3d..3c98dd660b 100644 --- a/improver/utilities/spatial.py +++ b/improver/utilities/spatial.py @@ -5,6 +5,7 @@ """ Provides support utilities.""" import copy +import warnings from typing import List, Optional, Tuple, Union import cartopy.crs as ccrs @@ -13,11 +14,13 @@ import numpy as np from cartopy.crs import CRS from cf_units import Unit -from iris.coords import AuxCoord, CellMethod +from iris.coord_systems import GeogCS +from iris.coords import AuxCoord, CellMethod, Coord from iris.cube import Cube, CubeList from numpy import ndarray from numpy.ma import MaskedArray from scipy.ndimage.filters import maximum_filter +from scipy.stats import circmean from improver import BasePlugin, PostProcessingPlugin from improver.metadata.amend import update_diagnostic_name @@ -164,6 +167,45 @@ class DifferenceBetweenAdjacentGridSquares(BasePlugin): individually. """ + @staticmethod + def _axis_wraps_around_meridian(axis: Coord, cube: Cube) -> bool: + """Returns true if the cube is 'circular' with the given axis wrapping around, i.e. if there + is a smooth transition between 180 degrees and -180 degrees on the axis. + + Args: + axis: + Axis to check for circularity. + cube: + The cube to which the axis belongs. + + Returns: + True if the axis wraps around the meridian; false otherwise. + """ + return axis.circular and axis == cube.coord(axis="x") + + @staticmethod + def _get_wrap_around_mean_point(points: ndarray) -> float: + """ + Calculates the midpoint between the two x coordinate points nearest the meridian. + + args: + points: + The x coordinate points of the cube. + + returns: + The x value of the midpoint between the two x coordinate points nearest the meridian. + """ + # The values of max and min azimuth doesn't matter as long as there is 360 degrees + # between them. + min_azimuth = -180 + max_azimuth = 180 + extra_mean_point = circmean([points[-1], points[0]], max_azimuth, min_azimuth) + extra_mean_point = np.round(extra_mean_point, 4) + if extra_mean_point < points[-1]: + # Ensures that the longitudinal coordinate is monotonically increasing + extra_mean_point += 360 + return extra_mean_point + @staticmethod def _update_metadata(diff_cube: Cube, coord_name: str, cube_name: str) -> None: """Rename cube, add attribute and cell method to describe difference. @@ -183,13 +225,11 @@ def _update_metadata(diff_cube: Cube, coord_name: str, cube_name: str) -> None: diff_cube.attributes["form_of_difference"] = "forward_difference" diff_cube.rename("difference_of_" + cube_name) - @staticmethod def create_difference_cube( - cube: Cube, coord_name: str, diff_along_axis: ndarray + self, cube: Cube, coord_name: str, diff_along_axis: ndarray ) -> Cube: """ - Put the difference array into a cube with the appropriate - metadata. + Put the difference array into a cube with the appropriate metadata. Args: cube: @@ -204,8 +244,21 @@ def create_difference_cube( Cube containing the differences calculated along the specified axis. """ - points = cube.coord(coord_name).points + axis = cube.coord(coord_name) + points = axis.points mean_points = (points[1:] + points[:-1]) / 2 + if self._axis_wraps_around_meridian(axis, cube): + if type(axis.coord_system) != GeogCS: + warnings.warn( + "DifferenceBetweenAdjacentGridSquares does not fully support cubes with " + "circular x-axis that do not use a geographic (i.e. latlon) coordinate system. " + "Such cubes will be handled as if they were not circular, meaning that the " + "differences cube returned will have one fewer points along the specified axis" + "than the input cube." + ) + else: + extra_mean_point = self._get_wrap_around_mean_point(points) + mean_points = np.hstack([mean_points, extra_mean_point]) # Copy cube metadata and coordinates into a new cube. # Create a new coordinate for the coordinate along which the @@ -226,8 +279,7 @@ def create_difference_cube( diff_cube.add_aux_coord(coord.copy(), dims) return diff_cube - @staticmethod - def calculate_difference(cube: Cube, coord_name: str) -> ndarray: + def calculate_difference(self, cube: Cube, coord_name: str) -> ndarray: """ Calculate the difference along the axis specified by the coordinate. @@ -242,8 +294,21 @@ def calculate_difference(cube: Cube, coord_name: str) -> ndarray: Array after the differences have been calculated along the specified axis. """ - diff_axis = cube.coord_dims(coord_name)[0] - diff_along_axis = np.diff(cube.data, axis=diff_axis) + diff_axis = cube.coord(name_or_coord=coord_name) + diff_axis_number = cube.coord_dims(coord_name)[0] + diff_along_axis = np.diff(cube.data, axis=diff_axis_number) + if self._axis_wraps_around_meridian(diff_axis, cube): + # Get wrap-around difference: + first_column = np.take(cube.data, indices=0, axis=diff_axis_number) + last_column = np.take(cube.data, indices=-1, axis=diff_axis_number) + wrap_around_diff = first_column - last_column + # Apply wrap-around difference vector to diff array: + if diff_axis_number == 0: + diff_along_axis = np.vstack([diff_along_axis, wrap_around_diff]) + elif diff_axis_number == 1: + diff_along_axis = np.hstack( + [diff_along_axis, wrap_around_diff.reshape([-1, 1])] + ) return diff_along_axis def process(self, cube: Cube) -> Tuple[Cube, Cube]: @@ -265,9 +330,8 @@ def process(self, cube: Cube) -> Tuple[Cube, Cube]: diffs = [] for axis in ["x", "y"]: coord_name = cube.coord(axis=axis).name() - diff_cube = self.create_difference_cube( - cube, coord_name, self.calculate_difference(cube, coord_name) - ) + difference = self.calculate_difference(cube, coord_name) + diff_cube = self.create_difference_cube(cube, coord_name, difference) self._update_metadata(diff_cube, coord_name, cube.name()) diffs.append(diff_cube) return tuple(diffs) diff --git a/improver_tests/utilities/test_DifferenceBetweenAdjacentGridSquares.py b/improver_tests/utilities/test_DifferenceBetweenAdjacentGridSquares.py index 1441dab528..4ff663261e 100644 --- a/improver_tests/utilities/test_DifferenceBetweenAdjacentGridSquares.py +++ b/improver_tests/utilities/test_DifferenceBetweenAdjacentGridSquares.py @@ -33,26 +33,74 @@ def setUp(self): def test_y_dimension(self): """Test differences calculated along the y dimension.""" points = self.cube.coord(axis="y").points - expected_y = (points[1:] + points[:-1]) / 2 + expected_y_coords = (points[1:] + points[:-1]) / 2 result = self.plugin.create_difference_cube( self.cube, "projection_y_coordinate", self.diff_in_y_array ) self.assertIsInstance(result, Cube) - self.assertArrayAlmostEqual(result.coord(axis="y").points, expected_y) + self.assertArrayAlmostEqual(result.coord(axis="y").points, expected_y_coords) self.assertArrayEqual(result.data, self.diff_in_y_array) def test_x_dimension(self): """Test differences calculated along the x dimension.""" diff_array = np.array([[1, 1], [2, 2], [5, 5]]) points = self.cube.coord(axis="x").points - expected_x = (points[1:] + points[:-1]) / 2 + expected_x_coords = (points[1:] + points[:-1]) / 2 result = self.plugin.create_difference_cube( self.cube, "projection_x_coordinate", diff_array ) self.assertIsInstance(result, Cube) - self.assertArrayAlmostEqual(result.coord(axis="x").points, expected_x) + self.assertArrayAlmostEqual(result.coord(axis="x").points, expected_x_coords) self.assertArrayEqual(result.data, diff_array) + def test_x_dimension_for_circular_latlon_cube(self): + """Test differences calculated along the x dimension for a cube which is circular in x.""" + test_cube_data = np.array([[1, 2, 3], [2, 4, 6], [5, 10, 15]]) + test_cube_x_grid_spacing = 120 + test_cube = set_up_variable_cube( + test_cube_data, + "latlon", + x_grid_spacing=test_cube_x_grid_spacing, + name="wind_speed", + units="m s-1", + ) + test_cube.coord(axis="x").circular = True + expected_diff_array = np.array([[1, 1, -2], [2, 2, -4], [5, 5, -10]]) + expected_x_coords = np.array( + [-60, 60, 180] + ) # Original data are at [-120, 0, 120], therefore differences are at [-60, 60, 180]. + result = self.plugin.create_difference_cube( + test_cube, "longitude", expected_diff_array + ) + self.assertIsInstance(result, Cube) + self.assertArrayAlmostEqual(result.coord(axis="x").points, expected_x_coords) + self.assertArrayEqual(result.data, expected_diff_array) + + def test_x_dimension_for_circular_latlon_cube_360_degree_coord(self): + """Test differences calculated along the x dimension for a cube which is circular in x.""" + test_cube_data = np.array([[1, 2, 3], [2, 4, 6], [5, 10, 15]]) + test_cube_x_grid_spacing = 120 + test_cube = set_up_variable_cube( + test_cube_data, + "latlon", + x_grid_spacing=test_cube_x_grid_spacing, + name="wind_speed", + units="m s-1", + ) + test_cube.coord(axis="x").bounds = [[0, 120], [120, 240], [240, 360]] + test_cube.coord(axis="x").points = [60, 120, 300] + test_cube.coord(axis="x").circular = True + expected_diff_array = np.array([[1, 1, -2], [2, 2, -4], [5, 5, -10]]) + expected_x_coords = np.array( + [90, 210, 360] + ) # Original data are at [60, 120, 300], therefore differences are at [90, 210, 360]. + result = self.plugin.create_difference_cube( + test_cube, "longitude", expected_diff_array + ) + self.assertIsInstance(result, Cube) + self.assertArrayAlmostEqual(result.coord(axis="x").points, expected_x_coords) + self.assertArrayEqual(result.data, expected_diff_array) + def test_othercoords(self): """Test that other coords are transferred properly""" time_coord = self.cube.coord("time") @@ -65,20 +113,40 @@ def test_othercoords(self): class Test_calculate_difference(IrisTest): - """Test the calculate_difference method.""" def setUp(self): """Set up cube.""" - data = np.array([[1, 2, 3], [2, 4, 6], [5, 10, 15]]) + data = np.array([[1, 2, 3, 4], [2, 4, 6, 8], [5, 10, 15, 20]]) self.cube = set_up_variable_cube( - data, name="wind_speed", units="m s-1", spatial_grid="equalarea", + data, "equalarea", name="wind_speed", units="m s-1", ) self.plugin = DifferenceBetweenAdjacentGridSquares() def test_x_dimension(self): """Test differences calculated along the x dimension.""" - expected = np.array([[1, 1], [2, 2], [5, 5]]) + expected = np.array([[1, 1, 1], [2, 2, 2], [5, 5, 5]]) + result = self.plugin.calculate_difference( + self.cube, self.cube.coord(axis="x").name() + ) + self.assertIsInstance(result, np.ndarray) + self.assertArrayEqual(result, expected) + + def test_x_dimension_wraps_around_meridian(self): + """Test differences calculated along the x dimension for a cube which is circular in x.""" + self.cube.coord(axis="x").circular = True + expected = np.array([[1, 1, 1, -3], [2, 2, 2, -6], [5, 5, 5, -15]]) + result = self.plugin.calculate_difference( + self.cube, self.cube.coord(axis="x").name() + ) + self.assertIsInstance(result, np.ndarray) + self.assertArrayEqual(result, expected) + + def test_x_dimension_wraps_around_meridian_cube_axes_flipped(self): + """Test differences calculated along the x dimension for a cube which is circular in x.""" + self.cube.coord(axis="x").circular = True + self.cube.transpose() + expected = np.array([[1, 1, 1, -3], [2, 2, 2, -6], [5, 5, 5, -15]]).transpose() result = self.plugin.calculate_difference( self.cube, self.cube.coord(axis="x").name() ) @@ -87,7 +155,7 @@ def test_x_dimension(self): def test_y_dimension(self): """Test differences calculated along the y dimension.""" - expected = np.array([[1, 2, 3], [3, 6, 9]]) + expected = np.array([[1, 2, 3, 4], [3, 6, 9, 12]]) result = self.plugin.calculate_difference( self.cube, self.cube.coord(axis="y").name() ) @@ -96,9 +164,11 @@ def test_y_dimension(self): def test_missing_data(self): """Test that the result is as expected when data is missing.""" - data = np.array([[1, 2, 3], [np.nan, 4, 6], [5, 10, 15]], dtype=np.float32) + data = np.array( + [[1, 2, 3, 4], [np.nan, 4, 6, 8], [5, 10, 15, 20]], dtype=np.float32 + ) self.cube.data = data - expected = np.array([[np.nan, 2, 3], [np.nan, 6, 9]]) + expected = np.array([[np.nan, 2, 3, 4], [np.nan, 6, 9, 12]]) result = self.plugin.calculate_difference( self.cube, self.cube.coord(axis="y").name() ) @@ -108,10 +178,13 @@ def test_missing_data(self): def test_masked_data(self): """Test that the result is as expected when data is masked.""" data = ma.array( - [[1, 2, 3], [2, 4, 6], [5, 10, 15]], mask=[[0, 0, 0], [1, 0, 0], [0, 0, 0]] + [[1, 2, 3, 4], [2, 4, 6, 8], [5, 10, 15, 20]], + mask=[[0, 0, 0, 0], [1, 0, 0, 0], [0, 0, 0, 0]], ) self.cube.data = data - expected = ma.array([[1, 2, 3], [3, 6, 9]], mask=[[1, 0, 0], [1, 0, 0]]) + expected = ma.array( + [[1, 2, 3, 4], [3, 6, 9, 12]], mask=[[1, 0, 0, 0], [1, 0, 0, 0]] + ) result = self.plugin.calculate_difference( self.cube, self.cube.coord(axis="y").name() ) @@ -121,7 +194,6 @@ def test_masked_data(self): class Test_process(IrisTest): - """Test the process method.""" def setUp(self): @@ -185,6 +257,11 @@ def test_3d_cube(self): self.assertIsInstance(result[1], iris.cube.Cube) self.assertArrayEqual(result[1].data, expected_y) + def test_circular_non_geographic_cube_raises_approprate_exception(self): + self.cube.coord(axis="x").circular = True + with self.assertRaises(ValueError): + self.plugin.process(self.cube) + if __name__ == "__main__": unittest.main() From 707436b7b6ebff19535a65843d0de11caa2f9842 Mon Sep 17 00:00:00 2001 From: mspelman07 <99179165+mspelman07@users.noreply.github.com> Date: Tue, 23 Jul 2024 12:54:16 +0100 Subject: [PATCH 2/7] Apply mask to cube (#2014) * Add apply_mask cli and function * remove print statement * Updates to make sure plugin doesn't require cubelist input * Update Docstring --------- Co-authored-by: Marcus Spelman --- improver/cli/apply_mask.py | 34 +++++++++ improver/utilities/mask.py | 72 ++++++++++++++++++ improver_tests/acceptance/SHA256SUMS | 4 + improver_tests/acceptance/test_apply_mask.py | 38 ++++++++++ improver_tests/utilities/test_mask.py | 78 ++++++++++++++++++++ 5 files changed, 226 insertions(+) create mode 100644 improver/cli/apply_mask.py create mode 100644 improver/utilities/mask.py create mode 100644 improver_tests/acceptance/test_apply_mask.py create mode 100644 improver_tests/utilities/test_mask.py diff --git a/improver/cli/apply_mask.py b/improver/cli/apply_mask.py new file mode 100644 index 0000000000..84f1c9959f --- /dev/null +++ b/improver/cli/apply_mask.py @@ -0,0 +1,34 @@ +#!/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 apply provided mask to cube data.""" + +from improver import cli + + +@cli.clizefy +@cli.with_output +def process(*cubes: cli.inputcube, mask_name: str, invert_mask: bool = "False"): + """ + Applies provided mask to cube data. The mask_name is used to extract the mask cube + from the input cubelist. The other cube in the cubelist is then masked using the + mask data. If invert_mask is True, the mask will be inverted before it is applied. + + Args: + cubes (iris.cube.CubeList): + A list of iris cubes that should contain exactly two cubes: a mask to be applied + and a cube to apply the mask to. The cubes should have the same dimensions. + mask_name (str): + The name of the cube containing the mask data. This should match with exactly one + of the cubes in the input cubelist. + invert_mask (bool): + Use to select whether the mask should be inverted before being applied to the data. + + Returns: + A cube with the mask applied to the data. The metadata will exactly match the input cube. + """ + from improver.utilities.mask import apply_mask + + return apply_mask(*cubes, mask_name=mask_name, invert_mask=invert_mask) diff --git a/improver/utilities/mask.py b/improver/utilities/mask.py new file mode 100644 index 0000000000..db3b4af94a --- /dev/null +++ b/improver/utilities/mask.py @@ -0,0 +1,72 @@ +# (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. +"""Module for applying mask to a cube.""" + +from typing import Union + +import iris +import numpy as np + +from improver.utilities.common_input_handle import as_cubelist +from improver.utilities.cube_checker import find_dimension_coordinate_mismatch +from improver.utilities.cube_manipulation import ( + enforce_coordinate_ordering, + get_coord_names, +) + + +def apply_mask( + *cubes: Union[iris.cube.CubeList, iris.cube.Cube], + mask_name: str, + invert_mask: bool = False, +) -> iris.cube.Cube: + """ + Apply a provided mask to a cube. If invert_mask is True, the mask will be inverted. + + Args: + cubes: + A list of iris cubes that should contain exactly two cubes: a mask and a cube + to apply the mask to. The cubes should have the same dimensions. + mask_name: + The name of the mask cube. It should match with exactly one of the cubes in + the input cubelist. + invert_mask: + If True, the mask will be inverted before it is applied. + Returns: + A cube with a mask applied to the data. + + Raises: + ValueError: If the number of cubes provided is not equal to 2. + ValueError: If the input cube and mask cube have different dimensions. + + """ + cubes = as_cubelist(*cubes) + cube_names = [cube.name() for cube in cubes] + if len(cubes) != 2: + raise ValueError( + f"""Two cubes are required for masking, a mask and the cube to be masked. + Provided cubes are {cube_names}""" + ) + + mask = cubes.extract_cube(mask_name) + cubes.remove(mask) + cube = cubes[0] + + # Ensure mask is in a boolean form and invert if requested + mask.data = mask.data.astype(np.bool) + if invert_mask: + mask.data = ~mask.data + + coord_list = get_coord_names(cube) + enforce_coordinate_ordering(mask, coord_list) + + # This error is required to stop the mask from being broadcasted to a new shape by numpy. When + # the mask and cube have different shapes numpy will try to broadcast the mask to be the same + # shape as the cube data. This might succeed but masks unexpected data points. + if find_dimension_coordinate_mismatch(cube, mask): + raise ValueError("Input cube and mask cube must have the same dimensions") + + cube.data = np.ma.array(cube.data, mask=mask.data) + return cube diff --git a/improver_tests/acceptance/SHA256SUMS b/improver_tests/acceptance/SHA256SUMS index a766e5d979..a7dacfd47e 100644 --- a/improver_tests/acceptance/SHA256SUMS +++ b/improver_tests/acceptance/SHA256SUMS @@ -70,6 +70,10 @@ c7eb9bab2ad43ac19ecc071730479a9f27a58992fcb22050743d34cdf2ad9639 ./apply-lapse- 6eb89848a2a36007d9a2210e2ff4cb459921502329eb4bb6ef0bde74701d0c27 ./apply-lapse-rate/realizations/enukx_orography.nc 25eeaa6b86b95354efd9406b1f9f1b5b43f0a7dd46f75fc5f3d144f52a30059e ./apply-lapse-rate/realizations/enukx_temperature.nc a747cf91ea8a2b720e1a22a9569cc1a2c22bb7ad28282e661127413470ba0392 ./apply-lapse-rate/realizations/kgo.nc +e455d391d15c30fe3cf3680fa0c355a3d4f0a9a84450f685f618655eb962ff52 ./apply-mask/kgo.nc +a3aca2f841ff96a7780cb98a5ccf874001905f18db0c6166dffcf116a3f1b319 ./apply-mask/kgo_inverted.nc +f33c0f3d10846ee69edd24fb318bebd4a2645363c89a434f56aeff24313fc886 ./apply-mask/mask.nc +9fba2a41268f30c530e9362bf15ef34c781d24b5e580589cec7715b1a6d25b27 ./apply-mask/wind_speed.nc 3a3d1b902e9e97f6ec7c74937e43dfce460c6599ef1890ac90f711959415f6ed ./apply-night-mask/global_basic/input.nc 65ba596b39ae589e1033b8e35fcb71f341b2fe64337fc2e4d88edff4331925e9 ./apply-night-mask/global_basic/kgo.nc c5860ce0801d097e1705130d6e3f49b19c43fa399041bd7726b7a1257335f2fb ./apply-night-mask/uk_basic/input.nc diff --git a/improver_tests/acceptance/test_apply_mask.py b/improver_tests/acceptance/test_apply_mask.py new file mode 100644 index 0000000000..050712bdbb --- /dev/null +++ b/improver_tests/acceptance/test_apply_mask.py @@ -0,0 +1,38 @@ +# (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 for the apply-mask 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) + + +@pytest.mark.parametrize("invert", [True, False]) +def test_apply_mask(tmp_path, invert): + """Test apply-mask CLI.""" + kgo_dir = acc.kgo_root() / "apply-mask/" + kgo_path = kgo_dir / "kgo.nc" + if invert: + kgo_path = kgo_dir / "kgo_inverted.nc" + + output_path = tmp_path / "output.nc" + args = [ + kgo_dir / "wind_speed.nc", + kgo_dir / "mask.nc", + "--mask-name", + "land_binary_mask", + "--invert-mask", + f"{invert}", + "--output", + output_path, + ] + run_cli(args) + acc.compare(output_path, kgo_path) diff --git a/improver_tests/utilities/test_mask.py b/improver_tests/utilities/test_mask.py new file mode 100644 index 0000000000..eddbfbe6fb --- /dev/null +++ b/improver_tests/utilities/test_mask.py @@ -0,0 +1,78 @@ +# (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 apply_mask. +""" + +import iris +import numpy as np +import pytest + +from improver.synthetic_data.set_up_test_cubes import set_up_variable_cube +from improver.utilities.cube_manipulation import enforce_coordinate_ordering +from improver.utilities.mask import apply_mask + + +@pytest.fixture +def wind_gust_cube(): + data = np.full((2, 3), 10) + return set_up_variable_cube( + data=data, attributes={"wind_gust_type": "10m_ratio"}, name="wind_gust" + ) + + +@pytest.fixture +def mask(): + data = np.array([[0, 0, 1], [1, 1, 0]]) + return set_up_variable_cube(data=data, name="land_sea_mask") + + +@pytest.mark.parametrize("invert_mask", [True, False]) +@pytest.mark.parametrize("switch_coord_order", [True, False]) +def test_basic(wind_gust_cube, mask, switch_coord_order, invert_mask): + """ + Test the basic functionality of the apply_mask plugin. Checks that the + mask is correctly applied and inverted if requested. Also checks plugin + can cope with different orderings of coordinates on the input cubes.""" + + expected_data = np.full((2, 3), 10) + expected_mask = np.array([[False, False, True], [True, True, False]]) + if switch_coord_order: + enforce_coordinate_ordering(wind_gust_cube, ["longitude", "latitude"]) + expected_data = expected_data.transpose() + expected_mask = expected_mask.transpose() + if invert_mask: + expected_mask = np.invert(expected_mask) + + input_list = [wind_gust_cube, mask] + + result = apply_mask( + iris.cube.CubeList(input_list), + mask_name="land_sea_mask", + invert_mask=invert_mask, + ) + + assert np.allclose(result.data, expected_data) + assert np.allclose(result.data.mask, expected_mask) + + +def test_different_dimensions(wind_gust_cube, mask): + """ Test that the function will raise an error if the mask cube has different + dimensions to other cube.""" + mask = mask[0] + input_list = [wind_gust_cube, mask] + with pytest.raises( + ValueError, match="Input cube and mask cube must have the same dimensions" + ): + apply_mask(iris.cube.CubeList(input_list), mask_name="land_sea_mask") + + +def test_too_many_cubes(wind_gust_cube, mask): + """ + Test that the function will raise an error if more than two cubes are provided. + """ + input_list = [wind_gust_cube, wind_gust_cube, wind_gust_cube] + with pytest.raises(ValueError, match="Two cubes are required for masking"): + apply_mask(iris.cube.CubeList(input_list), mask_name="land_sea_mask") From d6e643edd7063bdeadeccf5573ce6ce40dccb19b Mon Sep 17 00:00:00 2001 From: mspelman07 <99179165+mspelman07@users.noreply.github.com> Date: Wed, 24 Jul 2024 16:08:42 +0100 Subject: [PATCH 3/7] Adds masked_add to cube combiner (#2015) * Adds masked_add to cube combiner * Remove unnecessary import * formatting * add docstrings and simplifications * formatting --------- Co-authored-by: Marcus Spelman --- improver/cli/combine.py | 2 +- improver/cube_combiner.py | 31 +++++++++- .../cube_combiner/test_CubeCombiner.py | 57 ++++++++++++++++++- 3 files changed, 86 insertions(+), 4 deletions(-) diff --git a/improver/cli/combine.py b/improver/cli/combine.py index e297995473..0580b25103 100755 --- a/improver/cli/combine.py +++ b/improver/cli/combine.py @@ -32,7 +32,7 @@ def process( An iris CubeList to be combined. operation (str): An operation to use in combining input cubes. One of: - +, -, \*, add, subtract, multiply, min, max, mean + +, -, \*, add, subtract, multiply, min, max, mean, masked_add new_name (str): New name for the resulting dataset. broadcast (str): diff --git a/improver/cube_combiner.py b/improver/cube_combiner.py index f43326c9d5..6968ed8146 100644 --- a/improver/cube_combiner.py +++ b/improver/cube_combiner.py @@ -128,6 +128,32 @@ def process(self, *cubes: Union[Cube, CubeList]) -> Cube: return self.plugin(CubeList(filtered_cubes), self.new_name) +def masked_add( + masked_array: np.ma.MaskedArray, masked_array_2: np.ma.MaskedArray +) -> np.ma.MaskedArray: + """ + Operation to add two masked arrays treating masked points as 0. + + Args: + masked_array (numpy.ma.MaskedArray): + An array that may be masked. + masked_array_2 (numpy.ma.MaskedArray): + An array that may be masked. + + Returns: + numpy.ma.MaskedArray: + The sum of the two masked arrays with masked points treated as 0. + """ + new_array_1 = np.ma.filled(masked_array, 0) + new_array_2 = np.ma.filled(masked_array_2, 0) + + new_mask = np.ma.getmask(masked_array) * np.ma.getmask(masked_array_2) + + summed_cube = np.ma.MaskedArray(np.add(new_array_1, new_array_2), mask=new_mask) + + return summed_cube + + class CubeCombiner(BasePlugin): """Plugin for combining cubes using linear operators""" @@ -140,8 +166,9 @@ class CubeCombiner(BasePlugin): "multiply": np.multiply, "max": np.maximum, "min": np.minimum, - "mean": np.add, - } # mean is calculated in two steps: sum and normalise + "mean": np.add, # mean is calculated in two steps: sum and normalise + "masked_add": masked_add, # masked_add sums arrays but treats masked points as 0 + } def __init__( self, diff --git a/improver_tests/cube_combiner/test_CubeCombiner.py b/improver_tests/cube_combiner/test_CubeCombiner.py index da865bbfa6..baee8076f9 100644 --- a/improver_tests/cube_combiner/test_CubeCombiner.py +++ b/improver_tests/cube_combiner/test_CubeCombiner.py @@ -9,12 +9,13 @@ import iris import numpy as np +import pytest from iris.coords import CellMethod from iris.cube import Cube from iris.exceptions import CoordinateNotFoundError from iris.tests import IrisTest -from improver.cube_combiner import Combine, CubeCombiner +from improver.cube_combiner import Combine, CubeCombiner, masked_add from improver.synthetic_data.set_up_test_cubes import ( add_coordinate, set_up_probability_cube, @@ -497,3 +498,57 @@ def test_with_Combine(self): if __name__ == "__main__": unittest.main() + + +@pytest.fixture +def cube1(): + """Set up a probability cube with data for testing""" + data = np.full((1, 2, 2), 0.5, dtype=np.float32) + cube1 = set_up_probability_cube( + data, + np.array([0.001], dtype=np.float32), + variable_name="lwe_thickness_of_precipitation_amount", + time=datetime(2015, 11, 19, 0), + time_bounds=(datetime(2015, 11, 18, 23), datetime(2015, 11, 19, 0)), + frt=datetime(2015, 11, 18, 22), + ) + return cube1 + + +@pytest.fixture +def cube2(): + """Set up a second probability cube with data for testing""" + data = np.full((1, 2, 2), 0.6, dtype=np.float32) + cube2 = set_up_probability_cube( + data, + np.array([0.001], dtype=np.float32), + variable_name="lwe_thickness_of_precipitation_amount", + time=datetime(2015, 11, 19, 0), + time_bounds=(datetime(2015, 11, 18, 23), datetime(2015, 11, 19, 0)), + frt=datetime(2015, 11, 18, 22), + ) + return cube2 + + +@pytest.mark.parametrize("cube1_mask", [False, True]) +@pytest.mark.parametrize("cube2_mask", [False, True]) +def test_masked_add(cube1, cube2, cube1_mask, cube2_mask): + """Tests the plugin works with the masked_add option""" + mask = [[False, True], [False, False]] + expected_output = np.array(np.full((2, 2), 1.1, dtype=np.float32)) + expected_mask = [[False, False], [False, False]] + + if cube1_mask: + cube1.data = np.ma.MaskedArray(cube1.data, mask=mask) + expected_output[0][1] = 0.6 + if cube2_mask: + cube2.data = np.ma.MaskedArray(cube2.data, mask=mask) + if cube1_mask: + expected_mask = [[False, True], [False, False]] + expected_output[0][1] = 0.0 + else: + expected_output[0][1] = 0.5 + result = masked_add(cube1.data, cube2.data) + assert np.allclose(result.data, expected_output) + assert np.allclose(result.mask, expected_mask) + assert result.dtype == np.float32 From ba3f4e02e19df491460ec369c2f01b0cb277536a Mon Sep 17 00:00:00 2001 From: Carwyn Pelley Date: Mon, 29 Jul 2024 16:25:04 +0100 Subject: [PATCH 4/7] Plugin discovery (#2009) --- improver/api/__init__.py | 132 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 132 insertions(+) create mode 100644 improver/api/__init__.py diff --git a/improver/api/__init__.py b/improver/api/__init__.py new file mode 100644 index 0000000000..add901b267 --- /dev/null +++ b/improver/api/__init__.py @@ -0,0 +1,132 @@ +# (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. +# flake8: noqa +""" +This module contains the plugins for the IMPROVER project. This aids in discoverability +by making them available to a single flat namespace. This also protects end-users from +changes in structure to IMPROVER impacting their use of the plugins. +""" +from importlib import import_module + +# alphabetically sorted IMPROVER plugin lookup +PROCESSING_MODULES = { + "Accumulation": "improver.nowcasting.accumulation", + "AdjustLandSeaPoints": "improver.regrid.landsea", + "AdvectField": "improver.nowcasting.forecasting", + "AggregateReliabilityCalibrationTables": "improver.calibration.reliability_calibration", + "apply_mask": "improver.utilities.mask", + "ApplyBiasCorrection": "improver.calibration.simple_bias_correction", + "ApplyDecisionTree": "improver.categorical.decision_tree", + "ApplyDzRescaling": "improver.calibration.dz_rescaling", + "ApplyEMOS": "improver.calibration.ensemble_calibration", + "ApplyGriddedLapseRate": "improver.lapse_rate", + "ApplyNeighbourhoodProcessingWithAMask": "improver.nbhood.use_nbhood", + "ApplyOrographicEnhancement": "improver.nowcasting.utilities", + "ApplyRainForestsCalibration": "improver.calibration.rainforest_calibration", + "ApplyReliabilityCalibration": "improver.calibration.reliability_calibration", + "BaseNeighbourhoodProcessing": "improver.nbhood.nbhood", + "CalculateForecastBias": "improver.calibration.simple_bias_correction", + "CalibratedForecastDistributionParameters": "improver.calibration.ensemble_calibration", + "ChooseDefaultWeightsLinear": "improver.blending.weights", + "ChooseDefaultWeightsNonLinear": "improver.blending.weights", + "ChooseDefaultWeightsTriangular": "improver.blending.weights", + "ChooseWeightsLinear": "improver.blending.weights", + "CloudCondensationLevel": "improver.psychrometric_calculations.cloud_condensation_level", + "CloudTopTemperature": "improver.psychrometric_calculations.cloud_top_temperature", + "Combine": "improver.cube_combiner", + "ConstructReliabilityCalibrationTables": "improver.calibration.reliability_calibration", + "ContinuousRankedProbabilityScoreMinimisers": "improver.calibration.ensemble_calibration", + "ConvectionRatioFromComponents": "improver.precipitation_type.convection", + "ConvertProbabilitiesToPercentiles": "improver.ensemble_copula_coupling.ensemble_copula_coupling", + "CopyAttributes": "improver.utilities.copy_attributes", + "CorrectLandSeaMask": "improver.generate_ancillaries.generate_ancillary", + "CreateExtrapolationForecast": "improver.nowcasting.forecasting", + "CubeCombiner": "improver.cube_combiner", + "DayNightMask": "improver.utilities.solar", + "DifferenceBetweenAdjacentGridSquares": "improver.utilities.spatial", + "EnforceConsistentForecasts": "improver.utilities.forecast_reference_enforcement", + "EnsembleReordering": "improver.ensemble_copula_coupling.ensemble_copula_coupling", + "EstimateCoefficientsForEnsembleCalibration": "improver.calibration.ensemble_calibration", + "EstimateDzRescaling": "improver.calibration.dz_rescaling", + "ExpectedValue": "improver.expected_value", + "ExtendRadarMask": "improver.nowcasting.utilities", + "ExtractLevel": "improver.utilities.cube_extraction", + "ExtractSubCube": "improver.utilities.cube_extraction", + "FieldTexture": "improver.utilities.textural", + "FillRadarHoles": "improver.nowcasting.utilities", + "FreezingRain": "improver.precipitation_type.freezing_rain", + "FrictionVelocity": "improver.wind_calculations.wind_downscaling", + "GenerateClearskySolarRadiation": "improver.generate_ancillaries.generate_derived_solar_fields", + "GenerateOrographyBandAncils": "improver.generate_ancillaries.generate_ancillary", + "GenerateSolarTime": "improver.generate_ancillaries.generate_derived_solar_fields", + "GenerateTimeLaggedEnsemble": "improver.utilities.time_lagging", + "GenerateTopographicZoneWeights": "improver.generate_ancillaries.generate_topographic_zone_weights", + "GradientBetweenAdjacentGridSquares": "improver.utilities.spatial", + "HailFraction": "improver.precipitation_type.hail_fraction", + "HailSize": "improver.psychrometric_calculations.hail_size", + "HumidityMixingRatio": "improver.psychrometric_calculations.psychrometric_calculations", + "Integration": "improver.utilities.mathematical_operations", + "InterpolateUsingDifference": "improver.utilities.interpolation", + "LapseRate": "improver.lapse_rate", + "LightningFromCapePrecip": "improver.lightning", + "LightningMultivariateProbability_USAF2024": "improver.lightning", + "ManipulateReliabilityTable": "improver.calibration.reliability_calibration", + "MaxInTimeWindow": "improver.cube_combiner", + "MergeCubes": "improver.utilities.cube_manipulation", + "MergeCubesForWeightedBlending": "improver.blending.weighted_blend", + "MetaCloudCondensationLevel": "improver.psychrometric_calculations.cloud_condensation_level", + "MetaNeighbourhood": "improver.nbhood.nbhood", + "MetaWetBulbFreezingLevel": "improver.psychrometric_calculations.wet_bulb_temperature", + "ModalCategory": "improver.categorical.modal_code", + "NeighbourSelection": "improver.spotdata.neighbour_finding", + "NowcastLightning": "improver.nowcasting.lightning", + "OccurrenceBetweenThresholds": "improver.between_thresholds", + "OccurrenceWithinVicinity": "improver.utilities.spatial", + "OpticalFlow": "improver.nowcasting.optical_flow", + "OrographicEnhancement": "improver.orographic_enhancement", + "OrographicSmoothingCoefficients": "improver.generate_ancillaries.generate_orographic_smoothing_coefficients", + "PercentileConverter": "improver.percentile", + "PhaseChangeLevel": "improver.psychrometric_calculations.psychrometric_calculations", + "PostProcessingPlugin": "improver.__init__", + "PrecipPhaseProbability": "improver.psychrometric_calculations.precip_phase_probability", + "PystepsExtrapolate": "improver.nowcasting.pysteps_advection", + "RebadgePercentilesAsRealizations": "improver.ensemble_copula_coupling.ensemble_copula_coupling", + "RebadgeRealizationsAsPercentiles": "improver.ensemble_copula_coupling.ensemble_copula_coupling", + "RecursiveFilter": "improver.nbhood.recursive_filter", + "RegridLandSea": "improver.regrid.landsea", + "RegridWithLandSeaMask": "improver.regrid.landsea2", + "ResamplePercentiles": "improver.ensemble_copula_coupling.ensemble_copula_coupling", + "ResolveWindComponents": "improver.wind_calculations.wind_components", + "RoughnessCorrection": "improver.wind_calculations.wind_downscaling", + "SaturatedVapourPressureTable": "improver.generate_ancillaries.generate_svp_table", + "ShowerConditionProbability": "improver.precipitation_type.shower_condition_probability", + "SignificantPhaseMask": "improver.psychrometric_calculations.significant_phase_mask", + "SnowFraction": "improver.precipitation_type.snow_fraction", + "SnowSplitter": "improver.precipitation_type.snow_splitter", + "SpatiallyVaryingWeightsFromMask": "improver.blending.spatial_weights", + "SpotExtraction": "improver.spotdata.spot_extraction", + "SpotHeightAdjustment": "improver.spotdata.height_adjustment", + "SpotLapseRateAdjust": "improver.spotdata.apply_lapse_rate", + "SpotManipulation": "improver.spotdata.spot_manipulation", + "StandardiseMetadata": "improver.standardise", + "TemporalInterpolation": "improver.utilities.temporal_interpolation", + "Threshold": "improver.threshold", + "TriangularWeightedBlendAcrossAdjacentPoints": "improver.blending.blend_across_adjacent_points", + "VerticalUpdraught": "improver.wind_calculations.vertical_updraught", + "VisibilityCombineCloudBase": "improver.visibility.visibility_combine_cloud_base", + "WeightAndBlend": "improver.blending.calculate_weights_and_blend", + "WeightedBlendAcrossWholeDimension": "improver.blending.weighted_blend", + "WetBulbTemperature": "improver.psychrometric_calculations.wet_bulb_temperature", + "WetBulbTemperatureIntegral": "improver.psychrometric_calculations.wet_bulb_temperature", + "WindDirection": "improver.wind_calculations.wind_direction", + "WindGustDiagnostic": "improver.wind_calculations.wind_gust_diagnostic", +} + + +def __getattr__(name): + if name.startswith("__") and name.endswith("__"): + raise AttributeError(f"{name} is not a valid attribute") + mod = import_module(PROCESSING_MODULES[name]) + return getattr(mod, name) From e0daf3881299272964a15146f203ce3605356d73 Mon Sep 17 00:00:00 2001 From: Carwyn Pelley Date: Tue, 30 Jul 2024 08:48:35 +0100 Subject: [PATCH 5/7] BUG: Fix Scheduled Tests Sphinx-Pytest-Coverage (#2021) Pulling over pinned packages introduced to the `conda-forge.yml` into the `latest.yml` file. https://github.com/metoppv/improver/pull/2011 --- envs/latest.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/envs/latest.yml b/envs/latest.yml index 72c61d8f7f..dbe90fc694 100644 --- a/envs/latest.yml +++ b/envs/latest.yml @@ -45,3 +45,6 @@ dependencies: - sphinx-autodoc-typehints - sphinx_rtd_theme - threadpoolctl + # Pinned dependencies of dependencies + - pillow<=10.0.1 # https://github.com/metoppv/improver/issues/2010 + - pandas<=2.0.0 # https://github.com/metoppv/improver/issues/2010 From bc43de415499720980fd49c0e29075abb70ba215 Mon Sep 17 00:00:00 2001 From: SamGriffithsMO <122271903+SamGriffithsMO@users.noreply.github.com> Date: Thu, 1 Aug 2024 11:09:47 +0100 Subject: [PATCH 6/7] Migrate CLI functionality: simple_bias_correction (#2018) * refactored simple_bias_correction cli * docstring and contributing * changes for actions * sphinx building locally without error * whitespace change * bugfix: test upper_bound * changes based on review * review changes * review - doc-strings --- .mailmap | 1 + CONTRIBUTING.md | 1 + .../calibration/simple_bias_correction.py | 115 +++++++++++++----- improver/cli/apply_bias_correction.py | 30 +---- .../acceptance/test_apply_bias_correction.py | 67 ---------- .../test_ApplyBiasCorrection.py | 72 ++++++++++- 6 files changed, 155 insertions(+), 131 deletions(-) diff --git a/.mailmap b/.mailmap index 3959ea3131..4a4ca333f4 100644 --- a/.mailmap +++ b/.mailmap @@ -29,6 +29,7 @@ Meabh NicGuidhir <43375279+neilCrosswaite@users.noreply.github.com> Paul Abernethy Peter Jordan <52462411+mo-peterjordan@users.noreply.github.com> +Sam Griffiths <122271903+SamGriffithsMO@users.noreply.github.com> Shafiat Dewan <87321907+ShafiatDewan@users.noreply.github.com> <87321907+ShafiatDewan@users.noreply.github.com> Shubhendra Singh Chauhan Simon Jackson diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 22cf5c57cc..8dd2566698 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -49,6 +49,7 @@ below: - Zhiliang Fan (Bureau of Meteorology, Australia) - Ben Fitzpatrick (Met Office, UK) - Tom Gale (Bureau of Meteorology, Australia) + - Sam Griffiths (Met Office, UK) - Ben Hooper (Met Office, UK) - Aaron Hopkinson (Met Office, UK) - Kathryn Howard (Met Office, UK) diff --git a/improver/calibration/simple_bias_correction.py b/improver/calibration/simple_bias_correction.py index 23b71ebe42..a350ec3a77 100644 --- a/improver/calibration/simple_bias_correction.py +++ b/improver/calibration/simple_bias_correction.py @@ -4,7 +4,8 @@ # See LICENSE in the root of the repository for full licensing details. """Simple bias correction plugins.""" -from typing import Dict, Optional +import warnings +from typing import Dict, Optional, Union import iris import numpy.ma as ma @@ -12,6 +13,7 @@ from numpy import ndarray from improver import BasePlugin +from improver.calibration import add_warning_comment, split_forecasts_and_bias_files from improver.calibration.utilities import ( check_forecast_consistency, create_unified_frt_coord, @@ -23,6 +25,7 @@ create_new_diagnostic_cube, generate_mandatory_attributes, ) +from improver.utilities.common_input_handle import as_cubelist from improver.utilities.cube_manipulation import ( clip_cube_data, collapsed, @@ -248,11 +251,59 @@ class ApplyBiasCorrection(BasePlugin): the specified bias values. """ - def __init__(self): + def __init__( + self, + lower_bound: Optional[float] = None, + upper_bound: Optional[float] = None, + fill_masked_bias_values: bool = False, + ): """ Initialise class for applying simple bias correction. + + Args: + lower_bound: + A lower bound below which all values will be remapped to + after the bias correction step. + upper_bound: + An upper bound above which all values will be remapped to + after the bias correction step. + fill_masked_bias_values: + Flag to specify whether masked areas in the bias data + should be filled to an appropriate fill value. + """ + self._correction_method = apply_additive_correction + self._lower_bound = lower_bound + self._upper_bound = upper_bound + self._fill_masked_bias_values = fill_masked_bias_values + + def _split_forecasts_and_bias(self, cubes: CubeList): + """ + Wrapper for the split_forecasts_and_bias_files function. + + Args: + cubes: + Cubelist containing the input forecast and bias cubes. + + Return: + - Cube containing the forecast data to be bias-corrected. + - Cubelist containing the bias data to use in bias-correction. + Or None if no bias data is provided. """ - self.correction_method = apply_additive_correction + forecast_cube, bias_cubes = split_forecasts_and_bias_files(cubes) + + # Check whether bias data supplied, if not then return unadjusted input cube. + # This behaviour is to allow spin-up of the bias-correction terms. + if not bias_cubes: + msg = ( + "There are no forecast_error (bias) cubes provided for calibration. " + "The uncalibrated forecast will be returned." + ) + warnings.warn(msg) + forecast_cube = add_warning_comment(forecast_cube) + return forecast_cube, None + else: + bias_cubes = as_cubelist(bias_cubes) + return forecast_cube, bias_cubes def _get_mean_bias(self, bias_values: CubeList) -> Cube: """ @@ -302,6 +353,11 @@ def _check_forecast_bias_consistent( """Check that forecast and bias values are defined over the same valid-hour and forecast-period. + Checks that between the bias_data Cubes there is a common hour value for the + forecast_reference_time and single coordinate value for forecast_period. Then check + forecast Cube contains the same hour value for the forecast_reference_time and the + same forecast_period coordinate value. + Args: forecast: Cube containing forecast data to be bias-corrected. @@ -339,16 +395,8 @@ def _check_forecast_bias_consistent( "Forecast period differ between forecast and bias datasets." ) - def process( - self, - forecast: Cube, - bias: CubeList, - lower_bound: Optional[float] = None, - upper_bound: Optional[float] = None, - fill_masked_bias_values: Optional[bool] = False, - ) -> Cube: - """ - Apply bias correction using the specified bias values. + def process(self, *cubes: Union[Cube, CubeList],) -> Cube: + """Split then apply bias correction using the specified bias values. Where the bias data is defined point-by-point, the bias-correction will also be applied in this way enabling a form of statistical downscaling where coherent @@ -362,36 +410,37 @@ def process( filled using an appropriate fill value to leave the forecast data unchanged in the masked areas. + If no bias correction is provided, then the forecast is returned, unaltered. + Args: - forecast: - The cube to which bias correction is to be applied. - bias: - The cubelist containing the bias values for which to use in - the bias correction. - lower_bound: - A lower bound below which all values will be remapped to - after the bias correction step. - upper_bound: - An upper bound above which all values will be remapped to - after the bias correction step. - fill_masked_bias_values: - Flag to specify whether masked areas in the bias data - should be filled to an appropriate fill value. + cubes: + A list of cubes containing: + - A Cube containing the forecast to be calibrated. The input format is expected + to be realizations. + - A cube or cubelist containing forecast bias data over a specified + set of forecast reference times. If a list of cubes is passed in, each cube + should represent the forecast error for a single forecast reference time; the + mean value will then be evaluated over the forecast_reference_time coordinate. Returns: Bias corrected forecast cube. """ - self._check_forecast_bias_consistent(forecast, bias) - bias = self._get_mean_bias(bias) + cubes = as_cubelist(*cubes) + forecast, bias_cubes = self._split_forecasts_and_bias(cubes) + if bias_cubes is None: + return forecast + + self._check_forecast_bias_consistent(forecast, bias_cubes) + bias = self._get_mean_bias(bias_cubes) corrected_forecast = forecast.copy() - corrected_forecast.data = self.correction_method( - forecast, bias, fill_masked_bias_values + corrected_forecast.data = self._correction_method( + forecast, bias, self._fill_masked_bias_values ) - if (lower_bound is not None) or (upper_bound is not None): + if (self._lower_bound is not None) or (self._upper_bound is not None): corrected_forecast = clip_cube_data( - corrected_forecast, lower_bound, upper_bound + corrected_forecast, self._lower_bound, self._upper_bound ) return corrected_forecast diff --git a/improver/cli/apply_bias_correction.py b/improver/cli/apply_bias_correction.py index 5998c50014..2b9ed5ae38 100644 --- a/improver/cli/apply_bias_correction.py +++ b/improver/cli/apply_bias_correction.py @@ -12,7 +12,7 @@ @cli.clizefy @cli.with_output def process( - *input_cubes: cli.inputcube, + *cubes: cli.inputcube, lower_bound: float = None, upper_bound: float = None, fill_masked_bias_data: bool = False, @@ -32,7 +32,7 @@ def process( sensible post-bias correction. Args: - input_cubes (iris.cube.Cube or list of iris.cube.Cube): + cubes (iris.cube.Cube or list of iris.cube.Cube): A list of cubes containing: - A Cube containing the forecast to be calibrated. The input format is expected to be realizations. @@ -52,28 +52,8 @@ def process( iris.cube.Cube: Forecast cube with bias correction applied on a per member basis. """ - import warnings - - import iris - - from improver.calibration import add_warning_comment, split_forecasts_and_bias_files from improver.calibration.simple_bias_correction import ApplyBiasCorrection - forecast_cube, bias_cubes = split_forecasts_and_bias_files(input_cubes) - - # Check whether bias data supplied, if not then return unadjusted input cube. - # This behaviour is to allow spin-up of the bias-correction terms. - if not bias_cubes: - msg = ( - "There are no forecast_error (bias) cubes provided for calibration. " - "The uncalibrated forecast will be returned." - ) - warnings.warn(msg) - forecast_cube = add_warning_comment(forecast_cube) - return forecast_cube - else: - bias_cubes = iris.cube.CubeList(bias_cubes) - plugin = ApplyBiasCorrection() - return plugin.process( - forecast_cube, bias_cubes, lower_bound, upper_bound, fill_masked_bias_data - ) + return ApplyBiasCorrection(lower_bound, upper_bound, fill_masked_bias_data).process( + *cubes + ) diff --git a/improver_tests/acceptance/test_apply_bias_correction.py b/improver_tests/acceptance/test_apply_bias_correction.py index a86d83c8c1..d579559b90 100644 --- a/improver_tests/acceptance/test_apply_bias_correction.py +++ b/improver_tests/acceptance/test_apply_bias_correction.py @@ -134,70 +134,3 @@ def test_fill_masked_bias_data(tmp_path): ] run_cli(args) acc.compare(output_path, kgo_path) - - -def test_no_bias_file(tmp_path): - """ - Test case where no bias values are passed in. Expected behaviour is to - return the forecast value. - """ - kgo_dir = acc.kgo_root() / "apply-bias-correction" - fcst_path = kgo_dir / "20220814T0300Z-PT0003H00M-wind_speed_at_10m.nc" - kgo_path = kgo_dir / "fcst_with_comment" / "kgo.nc" - output_path = tmp_path / "output.nc" - args = [ - fcst_path, - "--output", - output_path, - ] - with pytest.warns(UserWarning, match=".*no forecast_error.*"): - run_cli(args) - acc.compare(output_path, fcst_path, exclude_attributes="comment") - acc.compare(output_path, kgo_path) - - -def test_missing_fcst_file(tmp_path): - """ - Test case where no forecast value has been passed in. This should raise - a ValueError. - """ - kgo_dir = acc.kgo_root() / "apply-bias-correction" - bias_file_path = ( - kgo_dir - / "single_bias_file" - / "bias_data" - / "20220813T0300Z-PT0003H00M-wind_speed_at_10m.nc" - ) - output_path = tmp_path / "output.nc" - args = [ - bias_file_path, - "--output", - output_path, - ] - with pytest.raises(ValueError, match="No forecast"): - run_cli(args) - - -def test_multiple_fcst_files(tmp_path): - """ - Test case where multiple forecast values are passed in. This should raise a - ValueError. - """ - kgo_dir = acc.kgo_root() / "apply-bias-correction" - fcst_path = kgo_dir / "20220814T0300Z-PT0003H00M-wind_speed_at_10m.nc" - bias_file_path = ( - kgo_dir - / "single_bias_file" - / "bias_data" - / "20220813T0300Z-PT0003H00M-wind_speed_at_10m.nc" - ) - output_path = tmp_path / "output.nc" - args = [ - fcst_path, - fcst_path, - bias_file_path, - "--output", - output_path, - ] - with pytest.raises(ValueError, match="Multiple forecast"): - run_cli(args) diff --git a/improver_tests/calibration/simple_bias_correction/test_ApplyBiasCorrection.py b/improver_tests/calibration/simple_bias_correction/test_ApplyBiasCorrection.py index f618224bbc..b807dcfdc2 100644 --- a/improver_tests/calibration/simple_bias_correction/test_ApplyBiasCorrection.py +++ b/improver_tests/calibration/simple_bias_correction/test_ApplyBiasCorrection.py @@ -4,6 +4,7 @@ # See LICENSE in the root of the repository for full licensing details. from datetime import datetime, timedelta +from unittest.mock import patch, sentinel import iris import numpy as np @@ -155,7 +156,7 @@ def test_apply_additive_correction( def test__init__(): """Test that the class functions are set to the expected values.""" plugin = ApplyBiasCorrection() - assert plugin.correction_method == apply_additive_correction + assert plugin._correction_method == apply_additive_correction @pytest.mark.parametrize("single_input_frt", (False, True)) @@ -257,7 +258,7 @@ def test_inconsistent_bias_forecast_inputs(forecast_cube, num_bias_inputs): @pytest.mark.parametrize("num_bias_inputs", (1, 30)) @pytest.mark.parametrize("single_input_frt", (False, True)) @pytest.mark.parametrize("lower_bound", (None, 1)) -@pytest.mark.parametrize("upper_bound", (None, 5)) +@pytest.mark.parametrize("upper_bound", (None, 4)) @pytest.mark.parametrize("masked_input_data", (True, False)) @pytest.mark.parametrize("masked_bias_data", (True, False)) @pytest.mark.parametrize("fill_masked_bias_data", (True, False)) @@ -282,12 +283,12 @@ def test_process( forecast_cube.data, dtype=forecast_cube.data.dtype ) forecast_cube.data.mask = MASK - result = ApplyBiasCorrection().process( - forecast_cube, - input_bias_cubelist, + result = ApplyBiasCorrection( lower_bound=lower_bound, + upper_bound=upper_bound, fill_masked_bias_values=fill_masked_bias_data, - ) + ).process(forecast_cube, input_bias_cubelist) + expected = TEST_FCST_DATA - MEAN_BIAS_DATA if fill_masked_bias_data and masked_bias_data: expected = np.where(MASK, TEST_FCST_DATA, expected) @@ -314,3 +315,62 @@ def test_process( # Check coords and attributes are consistent assert result.coords() == forecast_cube.coords() assert result.attributes == forecast_cube.attributes + + +class HaltExecution(Exception): + pass + + +@patch("improver.calibration.simple_bias_correction.as_cubelist") +def test_as_cubelist_called(mock_as_cubelist): + mock_as_cubelist.side_effect = HaltExecution + try: + ApplyBiasCorrection()(sentinel.cube1, sentinel.cube2, sentinel.cube3) + except HaltExecution: + pass + mock_as_cubelist.assert_called_once_with( + sentinel.cube1, sentinel.cube2, sentinel.cube3 + ) + + +def test_no_bias_file(forecast_cube): + """ + Test case where no bias values are passed in. Expected behaviour is to + return the forecast value. + """ + with pytest.warns(UserWarning, match=".*no forecast_error.*"): + result = ApplyBiasCorrection()(forecast_cube) + assert ( + "Warning: Calibration of this forecast has been attempted" + in result.attributes["comment"] + ) + + +def test_missing_fcst_file(): + """ + Test case where no forecast value has been passed in. This should raise + a ValueError. + """ + bias_cubes = generate_bias_cubelist( + 3, + last_valid_time=VALID_TIME + timedelta(hours=3), + single_frt_with_bounds=False, + ) + + with pytest.raises(ValueError, match="No forecast"): + ApplyBiasCorrection()(bias_cubes) + + +def test_multiple_fcst_files(forecast_cube): + """ + Test case where multiple forecast values are passed in. This should raise a + ValueError. + """ + bias_cubes = generate_bias_cubelist( + 3, + last_valid_time=VALID_TIME + timedelta(hours=3), + single_frt_with_bounds=False, + ) + + with pytest.raises(ValueError, match="Multiple forecast"): + ApplyBiasCorrection()(forecast_cube, forecast_cube, bias_cubes) From 8daf3350df3525022d1657c9b604abf6f31f9450 Mon Sep 17 00:00:00 2001 From: mspelman07 <99179165+mspelman07@users.noreply.github.com> Date: Fri, 2 Aug 2024 09:10:10 +0100 Subject: [PATCH 7/7] 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 Co-authored-by: Marcus Spelman --- .../cli/height_of_max_vertical_velocity.py | 46 +++++++++++ improver/utilities/cube_manipulation.py | 50 ++++++++++++ improver_tests/acceptance/SHA256SUMS | 3 + .../test_height_of_max_vertical_velocity.py | 31 +++++++ .../test_height_of_maximum.py | 80 +++++++++++++++++++ 5 files changed, 210 insertions(+) create mode 100644 improver/cli/height_of_max_vertical_velocity.py create mode 100644 improver_tests/acceptance/test_height_of_max_vertical_velocity.py create mode 100644 improver_tests/utilities/cube_manipulation/test_height_of_maximum.py diff --git a/improver/cli/height_of_max_vertical_velocity.py b/improver/cli/height_of_max_vertical_velocity.py new file mode 100644 index 0000000000..67972877b4 --- /dev/null +++ b/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) diff --git a/improver/utilities/cube_manipulation.py b/improver/utilities/cube_manipulation.py index 1513887314..6a94ce1b45 100644 --- a/improver/utilities/cube_manipulation.py +++ b/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 diff --git a/improver_tests/acceptance/SHA256SUMS b/improver_tests/acceptance/SHA256SUMS index a7dacfd47e..7b21941570 100644 --- a/improver_tests/acceptance/SHA256SUMS +++ b/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 diff --git a/improver_tests/acceptance/test_height_of_max_vertical_velocity.py b/improver_tests/acceptance/test_height_of_max_vertical_velocity.py new file mode 100644 index 0000000000..8ea42ec1f7 --- /dev/null +++ b/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) diff --git a/improver_tests/utilities/cube_manipulation/test_height_of_maximum.py b/improver_tests/utilities/cube_manipulation/test_height_of_maximum.py new file mode 100644 index 0000000000..dc8c3b9850 --- /dev/null +++ b/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)