From ab7e892ad908e5236d88414bf398d2e80760c54a Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Thu, 8 Feb 2024 11:49:35 +0000 Subject: [PATCH 01/19] Modified to return the later of the two input cubes along with the new interpolated values. --- improver/utilities/temporal_interpolation.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index d93220a153..70237abefa 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -56,6 +56,15 @@ class TemporalInterpolation(BasePlugin): cubes. This can be used to fill in missing data (e.g. for radar fields) or to ensure data is available at the required intervals when model data is not available at these times. + + The plugin will return the interpolated times and the later of the two + input times. This allows us to modify the input diagnostics if they + represent periods. The IMPROVER convention is that period diagnostics + have their time coordinate point at the end of the period. The later of + the two inputs therefore covers the period that has been broken down into + shorter periods by the interpolation and must itself be modified. The + result of this approach is that in a long run of lead-times, e.g. T+0 to + T+120 all the lead-times will be available except T+0. """ def __init__( @@ -177,6 +186,9 @@ def construct_time_list( break time_list.append(time_entry) + time_list.append(final_time) + time_list = sorted(set(time_list)) + return [("time", time_list)] @staticmethod @@ -393,7 +405,7 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: (final_time,) = iris_time_to_datetime(cube_t1.coord("time")) except CoordinateNotFoundError: msg = ( - "Cube provided to TemporalInterpolation contains no time " "coordinate." + "Cube provided to TemporalInterpolation contains no time coordinate." ) raise CoordinateNotFoundError(msg) except ValueError: From 2df6e07bd9f0ff69f94ec43c3452272ae720661c Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Thu, 8 Feb 2024 12:09:54 +0000 Subject: [PATCH 02/19] Add identification of period diagnostics using time bounds on inputs. Add bounds to time and forecast period coordinates of interpolated outputs if the inputs are period diagnostics. --- improver/utilities/temporal_interpolation.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index 70237abefa..1f8d7b6e43 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -111,12 +111,14 @@ def __init__( self.interval_in_minutes = interval_in_minutes self.times = times known_interpolation_methods = ["linear", "solar", "daynight"] + self.period_interpolation_methods = ["linear"] if interpolation_method not in known_interpolation_methods: raise ValueError( "TemporalInterpolation: Unknown interpolation " "method {}. ".format(interpolation_method) ) self.interpolation_method = interpolation_method + self.period_inputs = False def __repr__(self) -> str: """Represent the configured plugin instance as a string.""" @@ -400,6 +402,15 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: ) raise TypeError(msg) + if cube_t0.coord("time").bounds is not None: + if not self.interpolation_method in self.period_interpolation_methods: + raise ValueError( + "Period diagnostics can only be temporally interpolated " + f"using these methods: {self.period_interpolation_methods}.\n" + f"Currently selected method is: {self.interpolation_method}." + ) + self.period_inputs = True + try: (initial_time,) = iris_time_to_datetime(cube_t0.coord("time")) (final_time,) = iris_time_to_datetime(cube_t1.coord("time")) @@ -442,6 +453,12 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: interpolated_cube.data ) + # Add bounds to the time coordinates of the interpolated outputs + # if the inputs were period diagnostics. + if self.period_inputs: + for crd in ["time", "forecast_period"]: + interpolated_cube.coord(crd).guess_bounds(bound_position=1.0) + self.enforce_time_coords_dtype(interpolated_cube) interpolated_cubes = iris.cube.CubeList() if self.interpolation_method == "solar": From 51d974c7fde51f5be9cd20c67bf44c0f6cc63edd Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Thu, 8 Feb 2024 14:27:46 +0000 Subject: [PATCH 03/19] Add renormalisation for period accumulations. These need to be adjusted to ensure the temporally interpolated values give the same total over the whole period as the original inputs. --- improver/utilities/temporal_interpolation.py | 26 ++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index 1f8d7b6e43..4d9177937f 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -72,6 +72,8 @@ def __init__( interval_in_minutes: Optional[int] = None, times: Optional[List[datetime]] = None, interpolation_method: str = "linear", + accumulation: bool = False, + ) -> None: """ Initialise class. @@ -91,6 +93,12 @@ def __init__( interpolation_method: Method of interpolation to use. Default is linear. Only methods in known_interpolation_methods can be used. + accumulation: + Set True if the diagnostic being temporally interpolated is a + period accumulation. The output will be renormalised to ensure + that the total across the period constructed from the shorter + intervals matches the total across the period from the coarser + intervals. Raises: ValueError: If neither interval_in_minutes nor times are set. @@ -119,6 +127,7 @@ def __init__( ) self.interpolation_method = interpolation_method self.period_inputs = False + self.accumulation = accumulation def __repr__(self) -> str: """Represent the configured plugin instance as a string.""" @@ -453,12 +462,25 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: interpolated_cube.data ) - # Add bounds to the time coordinates of the interpolated outputs - # if the inputs were period diagnostics. if self.period_inputs: + # Add bounds to the time coordinates of the interpolated outputs + # if the inputs were period diagnostics. for crd in ["time", "forecast_period"]: interpolated_cube.coord(crd).guess_bounds(bound_position=1.0) + # If the input is an accumulation the total must be renormalised + # to avoid double counting. The input cube the contains an + # accumulation spanning the whole period is adjusted down to + # represent only a suitable fraction of the total time. The + # use of interpolation, rather than just splitting up the period, + # allows for trends in the data to add some variation to the + # different periods that are created. + if self.accumulation: + time_coord, = interpolated_cube.coord_dims("time") + interpolated_total = np.sum(interpolated_cube.data, axis=time_coord) + renormalisation = cube_t1.data / interpolated_total + interpolated_cube.data *= renormalisation + self.enforce_time_coords_dtype(interpolated_cube) interpolated_cubes = iris.cube.CubeList() if self.interpolation_method == "solar": From 58437b1ef9c66eddcd8ed4947a831905a93768c5 Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Thu, 8 Feb 2024 14:28:13 +0000 Subject: [PATCH 04/19] Add a check to ensure the guessed bounds of the interpolated outputs cover the entire period of the input. --- improver/utilities/temporal_interpolation.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index 4d9177937f..f9f20db77c 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -468,6 +468,17 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: for crd in ["time", "forecast_period"]: interpolated_cube.coord(crd).guess_bounds(bound_position=1.0) + # Check the interpolated output bounds cover the expected period. + lowest_bound = interpolated_cube.coord(crd).bounds[0][0] + highest_bound = interpolated_cube.coord(crd).bounds[-1][-1] + lowest_match = lowest_bound == cube_t1.coord(crd).bounds[0][0] + highest_match = highest_bound == cube_t1.coord(crd).bounds[-1][-1] + if not lowest_match or not highest_match: + raise ValueError( + "The interpolated periods do not cover the entire " + "period of the input." + ) + # If the input is an accumulation the total must be renormalised # to avoid double counting. The input cube the contains an # accumulation spanning the whole period is adjusted down to From f72fd0b651ad75717c75c07dc1067db746c38322 Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Fri, 9 Feb 2024 08:40:09 +0000 Subject: [PATCH 05/19] Partial rewrite of unit tests to use pytest. --- improver/utilities/temporal_interpolation.py | 9 - .../utilities/test_TemporalInterpolation.py | 928 +++++++----------- 2 files changed, 351 insertions(+), 586 deletions(-) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index f9f20db77c..7e13eef2cd 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -129,15 +129,6 @@ def __init__( self.period_inputs = False self.accumulation = accumulation - def __repr__(self) -> str: - """Represent the configured plugin instance as a string.""" - result = ( - "" - ) - return result.format( - self.interval_in_minutes, self.times, self.interpolation_method - ) - def construct_time_list( self, initial_time: datetime, final_time: datetime ) -> List[Tuple[str, List[datetime]]]: diff --git a/improver_tests/utilities/test_TemporalInterpolation.py b/improver_tests/utilities/test_TemporalInterpolation.py index c272995483..8e9ed258a5 100644 --- a/improver_tests/utilities/test_TemporalInterpolation.py +++ b/improver_tests/utilities/test_TemporalInterpolation.py @@ -31,10 +31,16 @@ """Unit tests for temporal utilities.""" import datetime +from datetime import datetime as dt import unittest +from typing import List, Tuple, Optional + import iris import numpy as np +import pytest + +from iris.cube import Cube, CubeList from iris.exceptions import CoordinateNotFoundError from iris.tests import IrisTest @@ -45,8 +51,10 @@ from improver.utilities.temporal_interpolation import TemporalInterpolation -def _grid_params(spatial_grid, npoints): - # Set domain corner and grid spacing +def _grid_params(spatial_grid: str, npoints: int) -> Tuple[Tuple[float, float], float]: + """Set domain corner and grid spacing for lat-lon or equal area + projections.""" + domain_corner = None grid_spacing = None if spatial_grid == "latlon": @@ -58,584 +66,350 @@ def _grid_params(spatial_grid, npoints): return domain_corner, grid_spacing -class Test__init__(IrisTest): - - """Test the __init__ method.""" - - def test_raises_error_with_no_keyword_args(self): - """Test __init__ raises a ValueError if both interval_in_minutes - and times keywords are unset.""" - msg = "TemporalInterpolation: One of" - with self.assertRaisesRegex(ValueError, msg): - TemporalInterpolation() - - def test_raises_error_with_both_keyword_args(self): - """Test __init__ raises a ValueError if both interval_in_minutes - and times keywords are both set.""" - msg = "TemporalInterpolation: Only one of" - with self.assertRaisesRegex(ValueError, msg): - TemporalInterpolation( - interval_in_minutes=60, times=[datetime.datetime(2017, 11, 1, 9)] - ) - - def test_unknown_method(self): - """Test __init__ raises a ValueError if method unknown.""" - msg = "TemporalInterpolation: Unknown interpolation method" - with self.assertRaisesRegex(ValueError, msg): - TemporalInterpolation( - interval_in_minutes=60, interpolation_method="invalid" - ) - - -class Test__repr__(IrisTest): - - """Test the repr method.""" - - def test_basic(self): - """Test that the __repr__ returns the expected string.""" - result = str(TemporalInterpolation(interval_in_minutes=60)) - msg = ( - "" - ) - self.assertEqual(result, msg) - - def test_solar(self): - """Test that the __repr__ returns the expected string for solar.""" - result = str( - TemporalInterpolation(interval_in_minutes=60, interpolation_method="solar") - ) - msg = ( - "" - ) - self.assertEqual(result, msg) - - def test_daynight(self): - """Test that the __repr__ returns the expected string for daynight.""" - result = str( - TemporalInterpolation( - interval_in_minutes=60, interpolation_method="daynight" - ) - ) - msg = ( - "" - ) - self.assertEqual(result, msg) - - -class Test_construct_time_list(IrisTest): - - """Test construction of time lists suitable for iris interpolation using - this function.""" - - def setUp(self): - """Set up the test inputs.""" - self.time_0 = datetime.datetime(2017, 11, 1, 3) - self.time_1 = datetime.datetime(2017, 11, 1, 9) - self.times = [] - for i in range(4, 9): - self.times.append(datetime.datetime(2017, 11, 1, i)) - self.expected = [("time", list(self.times))] - - def test_return_type(self): - """Test that a list is returned.""" - - result = TemporalInterpolation(interval_in_minutes=60).construct_time_list( - self.time_0, self.time_1 - ) - self.assertIsInstance(result, list) - - def test_list_from_interval_in_minutes(self): - """Test generating a list between two times using the - interval_in_minutes keyword to define the spacing.""" - - result = TemporalInterpolation(interval_in_minutes=60).construct_time_list( - self.time_0, self.time_1 - ) - self.assertEqual(self.expected, result) - - def test_non_equally_divisible_interval_in_minutes(self): - """Test an exception is raised when trying to generate a list of times - that would not divide the time range equally.""" - - msg = "interval_in_minutes of" - with self.assertRaisesRegex(ValueError, msg): - TemporalInterpolation(interval_in_minutes=61).construct_time_list( - self.time_0, self.time_1 - ) - - def test_list_from_existing_list(self): - """Test generating an iris interpolation suitable list from an - existing list of datetime objects.""" - - result = TemporalInterpolation(times=self.times).construct_time_list( - self.time_0, self.time_1 - ) - self.assertEqual(self.expected, result) - - def test_time_list_out_of_bounds(self): - """Test an exception is raised when trying to generate a list of times - using a pre-existing list that includes times outside the range of the - initial and final times.""" - - self.times.append(datetime.datetime(2017, 11, 1, 10)) - - msg = "List of times falls outside the range given by" - with self.assertRaisesRegex(ValueError, msg): - TemporalInterpolation(times=self.times).construct_time_list( - self.time_0, self.time_1 - ) - - -class Test_enforce_time_coords_dtype(IrisTest): - +def diagnostic_cube(time: dt, frt: dt, data: np.ndarray, spatial_grid: str, realizations: Optional[List] = None) -> Cube: + """Return a diagnostic cube containing the provided data. + + Args: + time: + A datetime object that gives the validity time of the cube. + frt: + The forecast reference time for the cube. + data: + The data to be contained in the cube. + spatial_grid: + Whether this is a lat-lon or equal areas projection. + Returns: + A diagnostic cube for use in testing. + """ + npoints = data.shape[0] + domain_corner, grid_spacing = _grid_params(spatial_grid, npoints) + + if realizations: + data = np.stack([data] * len(realizations)) + + return set_up_variable_cube( + data, + time=time, + frt=frt, + spatial_grid=spatial_grid, + domain_corner=domain_corner, + grid_spacing=grid_spacing, + realizations=realizations, + ) + + +def multi_time_cube(times: List, data: np.ndarray, spatial_grid: str, bounds: bool = False, realizations: Optional[List] = None) -> Cube: + """Return a multi-time diagnostic cube containing the provided data. + + Args: + times: + A list of datetime objects that gives the validity times for + the cube. + data: + The data to be contained in the cube. If the cube is 3-D the + leading dimension should be the same size and the list of times + and will be sliced to associate each slice with each time. + spatial_grid: + Whether this is a lat-lon or equal areas projection. + bounds: + If True return time coordinates with time bounds. + Returns: + A diagnostic cube for use in testing. + """ + cubes = CubeList() + if data.ndim == 2: + data = np.stack([data] * len(times)) + + frt = sorted(times)[0] - (times[1] - times[0]) # Such that guess bounds is +ve + for time, data_slice in zip(times, data): + cubes.append(diagnostic_cube(time, frt, data_slice, spatial_grid, realizations=realizations)) + + cube = cubes.merge_cube() + + if bounds: + for crd in ["time", "forecast_period"]: + cube.coord(crd).guess_bounds(bound_position=1.0) + return cube + + +def non_standard_times(times: List, data: np.ndarray, spatial_grid: str, bounds: bool = False) -> Cube: + """Return a multi-time diagnostic cube containing the provided data. + The units of the time dimensions are made non-standards compliant. + + Args: + times: + A list of datetime objects that gives the validity times for + the cube. + data: + The data to be contained in the cube. If the cube is 3-D the + leading dimension should be the same size and the list of times + and will be sliced to associate each slice with each time. + spatial_grid: + Whether this is a lat-lon or equal areas projection. + bounds: + If True return time coordinates with time bounds. + Returns: + A diagnostic cube for use in testing. + """ + cube = multi_time_cube(times, data, spatial_grid, bounds=bounds) + + epoch = "hours since 1970-01-01 00:00:00" + for crd in ["time", "forecast_reference_time"]: + cube.coord(crd).convert_units(epoch) + cube.coord(crd).points = cube.coord(crd).points.astype(np.int32) + + cube.coord("forecast_period").convert_units("hours") + cube.coord("forecast_period").points.astype(np.float32) + + return cube + + +@pytest.fixture +def solar_expected(): + """Return the expected values for the solar interpolation tests.""" + return np.array( + [ + [0.02358028, 0.15887623, 0.2501732, 0.32049885, 0.3806127], + [0.0, 0.09494493, 0.21051247, 0.2947393, 0.36431003], + [0.0, 0.0, 0.11747278, 0.23689085, 0.32841164], + [0.0, 0.0, 0.0, 0.0, 0.15872595], + [0.0, 0.0, 0.0, 0.0, 0.0], + ] + ) + + +@pytest.fixture +def daynight_mask(): + """The mask matching the terminator position for the daynight + interpolation tests.""" + return np.array( + [ + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [0, 0, 1, 1, 1, 1, 1, 1, 1, 1], + [0, 0, 0, 1, 1, 1, 1, 1, 1, 1], + [0, 0, 0, 0, 1, 1, 1, 1, 1, 1], + [0, 0, 0, 0, 0, 0, 1, 1, 1, 1], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ] + ) + + +@pytest.mark.parametrize( + "kwargs,exception", + [ + ({}, "TemporalInterpolation: One of"), # No target times defined + ({"interval_in_minutes":60, "times":[datetime.datetime(2017, 11, 1, 9)]}, "TemporalInterpolation: Only one of"), # Two methods of defining targets used + ({"interval_in_minutes":60, "interpolation_method":"invalid"}, "TemporalInterpolation: Unknown interpolation method"), # Invalid interpolation method requested + ] +) +def test__init__(kwargs, exception): + """Test exceptions raised by the __init__ method.""" + with pytest.raises(ValueError, match=exception): + TemporalInterpolation(**kwargs) + + +@pytest.mark.parametrize( + "kwargs,exception", + [ + ({"interval_in_minutes": 60}, None), # Generate times between bounds using interval + ({"times": None}, None), # Use the expected times as the input + ({"times": datetime.datetime(2017, 11, 1, 10)}, "List of times falls outside the range given by"), # Use the expected times, plus another outside the range as the input + ({"interval_in_minutes":61}, "interval_in_minutes of"), # Use an invalid interval + ] +) +def test_construct_time_list(kwargs, exception): + """Test construction of target times using various inputs and testing + exceptions that can be raised.""" + time_0 = datetime.datetime(2017, 11, 1, 3) + time_1 = datetime.datetime(2017, 11, 1, 9) + + # Expected times are all those interpolated to and time_1. + times = [] + for i in range(4, 10): + times.append(datetime.datetime(2017, 11, 1, i)) + expected = [("time", list(times))] + + # If a times kwarg is supplied, populate the value with the default + # times plus any others specified in the kwarg. + try: + target_times = times.copy() + target_times.append(kwargs["times"]) if kwargs["times"] is not None else target_times + except KeyError: + pass + else: + kwargs["times"] = target_times + + plugin = TemporalInterpolation(**kwargs) + + # If an exception is provided as a kwarg test for this. + if exception is not None: + with pytest.raises(ValueError, match=exception): + plugin.construct_time_list(time_0, time_1) + else: + result = plugin.construct_time_list(time_0, time_1) + assert isinstance(result, list) + assert expected == result + + +@pytest.mark.parametrize("bounds", (False, True)) +def test_enforce_time_coords_dtype(bounds): """Test that the datatypes and units of the time, forecast_reference_time and forecast_period coordinates have been enforced.""" - def setUp(self): - """Set up the test inputs.""" - time_start = datetime.datetime(2017, 11, 1, 3) - time_mid = datetime.datetime(2017, 11, 1, 6) - time_end = datetime.datetime(2017, 11, 1, 9) - self.npoints = 10 - - domain_corner, grid_spacing = _grid_params("latlon", self.npoints) - - data_time_0 = np.ones((self.npoints, self.npoints), dtype=np.float32) - cube_time_0 = set_up_variable_cube( - data_time_0, - time=time_start, - frt=time_start, - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - cube_times = add_coordinate( - cube_time_0.copy(), - [time_start, time_mid, time_end], - "time", - is_datetime=True, - ) - # Convert units and datatypes, so that they are non-standard. - cube_times.coord("time").convert_units("hours since 1970-01-01 00:00:00") - cube_times.coord("time").points = cube_times.coord("time").points.astype( - np.int32 - ) - cube_times.coord("forecast_reference_time").convert_units( - "hours since 1970-01-01 00:00:00" - ) - cube_times.coord("forecast_reference_time").points = cube_times.coord( - "forecast_reference_time" - ).points.astype(np.int32) - cube_times.coord("forecast_period").convert_units("hours") - cube_times.coord("forecast_period").points.astype(np.float32) - self.cube = cube_times - - self.coord_dtypes = { - "time": np.int64, - "forecast_reference_time": np.int64, - "forecast_period": np.int32, - } - self.coord_units = { - "time": "seconds since 1970-01-01 00:00:00", - "forecast_reference_time": "seconds since 1970-01-01 00:00:00", - "forecast_period": "seconds", - } - - def test_check_points(self): - """Test that a cube is returned with the desired types for the points - of the time, forecast_reference_time and forecast_period - coordinates.""" - plugin = TemporalInterpolation(interval_in_minutes=60) - result = plugin.enforce_time_coords_dtype(self.cube.copy()) - self.assertIsInstance(result, iris.cube.Cube) - # All coordinates converted to the desired units and datatypes. - # Check time coordinate. - self.assertEqual(result.coord("time").points.dtype, self.coord_dtypes["time"]) - self.assertEqual(result.coord("time").units, self.coord_units["time"]) - # Check forecast_reference_time. - self.assertEqual( - result.coord("forecast_reference_time").points.dtype, - self.coord_dtypes["forecast_reference_time"], - ) - self.assertEqual( - result.coord("forecast_reference_time").units, - self.coord_units["forecast_reference_time"], - ) - # Check forecast_period. - self.assertEqual( - result.coord("forecast_period").points.dtype, - self.coord_dtypes["forecast_period"], - ) - self.assertEqual( - result.coord("forecast_period").units, self.coord_units["forecast_period"] - ) - - def test_check_bounds(self): - """Test that a cube is returned with the desired types when - the time and forecast_period coordinates have bounds.""" - plugin = TemporalInterpolation(interval_in_minutes=60) - cube = self.cube - # Use of guess_bounds converts datatype to float64. - cube.coord("time").guess_bounds() - cube.coord("forecast_period").guess_bounds() - - result = plugin.enforce_time_coords_dtype(cube.copy()) - self.assertIsInstance(result, iris.cube.Cube) - # All coordinates including bounds converted to the - # desired units and datatypes. - # Check time coordinate. - self.assertEqual(result.coord("time").points.dtype, self.coord_dtypes["time"]) - self.assertEqual(result.coord("time").bounds.dtype, self.coord_dtypes["time"]) - self.assertEqual(result.coord("time").units, self.coord_units["time"]) - # Check forecast_reference_time coordinate. - self.assertEqual( - result.coord("forecast_reference_time").points.dtype, - self.coord_dtypes["forecast_reference_time"], - ) - self.assertEqual( - result.coord("forecast_reference_time").units, - self.coord_units["forecast_reference_time"], - ) - # Check forecast_period coordinate. - self.assertEqual( - result.coord("forecast_period").points.dtype, - self.coord_dtypes["forecast_period"], - ) - self.assertEqual( - result.coord("forecast_period").bounds.dtype, - self.coord_dtypes["forecast_period"], - ) - self.assertEqual( - result.coord("forecast_period").units, self.coord_units["forecast_period"] - ) - - -class Test_calc_sin_phi(IrisTest): - - """Test Calculate sin of solar elevation.""" - - def test_sin_phi(self): - """Test that the function returns the values expected.""" - latitudes = np.array([50.0, 50.0, 50.0]) - longitudes = np.array([-5.0, 0.0, 5.0]) - dtval = datetime.datetime(2017, 1, 11, 8) - expected_array = np.array([-0.05481607, -0.00803911, 0.03659632]) - plugin = TemporalInterpolation( - interval_in_minutes=60, interpolation_method="solar" - ) - result = plugin.calc_sin_phi(dtval, latitudes, longitudes) - self.assertIsInstance(result, np.ndarray) - self.assertArrayAlmostEqual(result, expected_array) - - -class Test_calc_lats_lons(IrisTest): - - """Test Calculate lats and lons.""" - - def setUp(self): - """Set up the test inputs.""" - time_start = datetime.datetime(2017, 11, 1, 3) - time_mid = datetime.datetime(2017, 11, 1, 6) - time_end = datetime.datetime(2017, 11, 1, 9) - self.npoints = 3 - - equalarea_domain_corner, equalarea_grid_spacing = _grid_params( - "equalarea", self.npoints - ) - latlon_domain_corner, latlon_grid_spacing = _grid_params("latlon", self.npoints) - - data_time_0 = np.ones((self.npoints, self.npoints), dtype=np.float32) - cube_time_0 = set_up_variable_cube( - data_time_0, - time=time_start, - frt=time_start, - domain_corner=latlon_domain_corner, - grid_spacing=latlon_grid_spacing, - ) - self.cube = add_coordinate( - cube_time_0, [time_start, time_mid, time_end], "time", is_datetime=True - ) - cube_time_0_equalarea = set_up_variable_cube( - data_time_0, - time=time_start, - frt=time_start, - spatial_grid="equalarea", - domain_corner=equalarea_domain_corner, - grid_spacing=equalarea_grid_spacing, - ) - self.cube_equalarea = add_coordinate( - cube_time_0_equalarea, - [time_start, time_mid, time_end], - "time", - is_datetime=True, - ) - - def test_lat_lon(self): - """Test that the function returns the lats and lons expected.""" - expected_lons = np.array( - [[-20.0, 0.0, 20.0], [-20.0, 0.0, 20.0], [-20.0, 0.0, 20.0]] - ) - expected_lats = np.array( - [[40.0, 40.0, 40.0], [60.0, 60.0, 60.0], [80.0, 80.0, 80.0]] - ) - plugin = TemporalInterpolation( - interval_in_minutes=60, interpolation_method="solar" - ) - result_lats, result_lons = plugin.calc_lats_lons(self.cube) - self.assertIsInstance(result_lats, np.ndarray) - self.assertEqual(result_lats.shape, (3, 3)) - self.assertIsInstance(result_lons, np.ndarray) - self.assertEqual(result_lons.shape, (3, 3)) - self.assertArrayAlmostEqual(result_lats, expected_lats) - self.assertArrayAlmostEqual(result_lons, expected_lons) - - def test_x_y(self): - """Test that the function returns the lats and lons expected.""" - expected_lats = np.array( - [ - [53.84618597, 53.99730779, 53.93247526], - [56.82670954, 56.99111356, 56.9205672], - [59.8045105, 59.98499383, 59.90752513], - ] - ) - - expected_lons = np.array( - [ - [-8.58580705, -3.51660018, 1.56242662], - [-9.06131306, -3.59656346, 1.88105082], - [-9.63368459, -3.69298822, 2.26497216], - ] - ) - - plugin = TemporalInterpolation( - interval_in_minutes=60, interpolation_method="solar" - ) - result_lats, result_lons = plugin.calc_lats_lons(self.cube_equalarea) - self.assertIsInstance(result_lats, np.ndarray) - self.assertEqual(result_lats.shape, (3, 3)) - self.assertIsInstance(result_lons, np.ndarray) - self.assertEqual(result_lons.shape, (3, 3)) - self.assertArrayAlmostEqual(result_lats, expected_lats) - self.assertArrayAlmostEqual(result_lons, expected_lons) - - -class Test_solar_interpolation(IrisTest): - - """Test Solar interpolation.""" - - def setUp(self): - """Set up the test inputs spanning sunrise.""" - self.time_0 = datetime.datetime(2017, 11, 1, 6) - self.time_mid = datetime.datetime(2017, 11, 1, 8) - self.time_1 = datetime.datetime(2017, 11, 1, 10) - self.npoints = 5 - self.expected = np.array( - [ - [0.02358028, 0.15887623, 0.2501732, 0.32049885, 0.3806127], - [0.0, 0.09494493, 0.21051247, 0.2947393, 0.36431003], - [0.0, 0.0, 0.11747278, 0.23689085, 0.32841164], - [0.0, 0.0, 0.0, 0.0, 0.15872595], - [0.0, 0.0, 0.0, 0.0, 0.0], - ] - ) - - domain_corner, grid_spacing = _grid_params("latlon", self.npoints) - - data_time_0 = np.zeros((self.npoints, self.npoints), dtype=np.float32) - data_time_1 = np.ones((self.npoints, self.npoints), dtype=np.float32) - data_time_mid = np.ones((self.npoints, self.npoints), dtype=np.float32) - cube_time_0 = set_up_variable_cube( - data_time_0, - time=self.time_0, - frt=self.time_0, - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - cube_time_1 = set_up_variable_cube( - data_time_1, - time=self.time_1, - frt=self.time_0, - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - cubes = iris.cube.CubeList([cube_time_0, cube_time_1]) - self.cube = cubes.merge_cube() - interp_cube = set_up_variable_cube( - data_time_mid, - time=self.time_mid, - frt=self.time_0, - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - self.interpolated_cube = iris.util.new_axis(interp_cube, "time") - data_time_0_ens = np.zeros((3, self.npoints, self.npoints), dtype=np.float32) - data_time_1_ens = np.ones((3, self.npoints, self.npoints), dtype=np.float32) - data_time_mid_ens = np.ones((3, self.npoints, self.npoints), dtype=np.float32) - cube_time_0_ens = set_up_variable_cube( - data_time_0_ens, - time=self.time_0, - frt=self.time_0, - realizations=[0, 1, 2], - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - cube_time_1_ens = set_up_variable_cube( - data_time_1_ens, - time=self.time_1, - frt=self.time_0, - realizations=[0, 1, 2], - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - interp_cube_ens = set_up_variable_cube( - data_time_mid_ens, - time=self.time_mid, - frt=self.time_0, - realizations=[0, 1, 2], - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - self.interpolated_cube_ens = iris.util.new_axis(interp_cube_ens, "time") - cubes_ens = iris.cube.CubeList([cube_time_0_ens, cube_time_1_ens]) - self.cube_ens = cubes_ens.merge_cube() - - def test_return_type(self): - """Test that an iris cubelist is returned.""" - - plugin = TemporalInterpolation( - interpolation_method="solar", times=[self.time_mid] - ) - result = plugin.solar_interpolate(self.cube, self.interpolated_cube) - self.assertIsInstance(result, iris.cube.CubeList) - - def test_solar_interpolation(self): - """Test interpolating using solar method works correctly.""" - - expected_time = [1509523200] - expected_fp = 2 * 3600 - plugin = TemporalInterpolation( - interpolation_method="solar", times=[self.time_mid] - ) - (result,) = plugin.solar_interpolate(self.cube, self.interpolated_cube) - self.assertArrayAlmostEqual(result.data, self.expected) - self.assertArrayAlmostEqual(result.coord("time").points, expected_time) - self.assertAlmostEqual(result.coord("forecast_period").points[0], expected_fp) - - def test_solar_interpolation_shape(self): - """Test interpolating using solar method with len(shape) >= 3 - works correctly.""" - - expected_time = [1509523200] - expected_fp = 2 * 3600 - plugin = TemporalInterpolation( - interpolation_method="solar", times=[self.time_mid] - ) - (result,) = plugin.solar_interpolate(self.cube_ens, self.interpolated_cube_ens) - - self.assertArrayEqual(result.data.shape, (3, 5, 5)) - self.assertArrayAlmostEqual(result.data[0], self.expected) - self.assertArrayAlmostEqual(result.coord("time").points, expected_time) - self.assertAlmostEqual(result.coord("forecast_period").points[0], expected_fp) - - -class Test_daynight_interpolation(IrisTest): - - """Test daynight interpolation.""" - - def setUp(self): - """Set up the test inputs spanning sunrise.""" - self.time_0 = datetime.datetime(2017, 11, 1, 6) - self.time_mid = datetime.datetime(2017, 11, 1, 8) - self.npoints = 10 - self.daynight_mask = np.array( - [ - [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], - [0, 1, 1, 1, 1, 1, 1, 1, 1, 1], - [0, 1, 1, 1, 1, 1, 1, 1, 1, 1], - [0, 0, 1, 1, 1, 1, 1, 1, 1, 1], - [0, 0, 0, 1, 1, 1, 1, 1, 1, 1], - [0, 0, 0, 0, 1, 1, 1, 1, 1, 1], - [0, 0, 0, 0, 0, 0, 1, 1, 1, 1], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 1], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - ] - ) - - domain_corner, grid_spacing = _grid_params("latlon", self.npoints) - - data_time_mid = np.ones((self.npoints, self.npoints), dtype=np.float32) * 4 - - interp_cube = set_up_variable_cube( - data_time_mid, - time=self.time_mid, - frt=self.time_0, - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - self.interpolated_cube = iris.util.new_axis(interp_cube, "time") - - data_time_mid_ens = ( - np.ones((3, self.npoints, self.npoints), dtype=np.float32) * 4 - ) - interp_cube_ens = set_up_variable_cube( - data_time_mid_ens, - time=self.time_mid, - frt=self.time_0, - realizations=[0, 1, 2], - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - self.interpolated_cube_ens = iris.util.new_axis(interp_cube_ens, "time") - - def test_return_type(self): - """Test that an iris cubelist is returned.""" - - plugin = TemporalInterpolation( - interpolation_method="daynight", times=[self.time_mid] - ) - result = plugin.daynight_interpolate(self.interpolated_cube) - self.assertIsInstance(result, iris.cube.CubeList) - - def test_daynight_interpolation(self): - """Test interpolating to the a point where the daynight - mask is not all zero.""" - - expected_data = np.ones((self.npoints, self.npoints)) * 4 - index = np.where(self.daynight_mask == 0) - expected_data[index] = 0.0 - expected_time = [1509523200] - expected_fp = 2 * 3600 - plugin = TemporalInterpolation( - interpolation_method="daynight", times=[self.time_mid] - ) - (result,) = plugin.daynight_interpolate(self.interpolated_cube) - self.assertArrayAlmostEqual(expected_data, result.data) - self.assertArrayAlmostEqual(result.coord("time").points, expected_time) - self.assertAlmostEqual(result.coord("forecast_period").points[0], expected_fp) - - def test_daynight_interpolation_ens(self): - """Test interpolating to the a point where the daynight - mask is not all zero and the len(shape) of the cube > 2.""" - - expected_data_grid = np.ones((self.npoints, self.npoints)) * 4 - index = np.where(self.daynight_mask == 0) - expected_data_grid[index] = 0.0 - expected_data = np.repeat(expected_data_grid[np.newaxis, :, :], 3, axis=0) - expected_time = [1509523200] - expected_fp = 2 * 3600 - plugin = TemporalInterpolation( - interpolation_method="daynight", times=[self.time_mid] - ) - (result,) = plugin.daynight_interpolate(self.interpolated_cube_ens) - self.assertArrayAlmostEqual(expected_data, result.data) - self.assertArrayAlmostEqual(result.coord("time").points, expected_time) - self.assertAlmostEqual(result.coord("forecast_period").points[0], expected_fp) + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 6, 9]] + data = np.ones((10, 10), dtype=np.float32) + cube = non_standard_times(times, data, "latlon", bounds=bounds) + + expected_coord_dtypes = { + "time": np.int64, + "forecast_reference_time": np.int64, + "forecast_period": np.int32, + } + expected_coord_units = { + "time": "seconds since 1970-01-01 00:00:00", + "forecast_reference_time": "seconds since 1970-01-01 00:00:00", + "forecast_period": "seconds", + } + + plugin = TemporalInterpolation(interval_in_minutes=60) + result = plugin.enforce_time_coords_dtype(cube.copy()) + + assert isinstance(result, Cube) + + for crd, expected in expected_coord_dtypes.items(): + assert result.coord(crd).points.dtype == expected + try: + assert result.coord(crd).bounds.dtype == expected + except AttributeError: + pass + for crd, expected in expected_coord_units.items(): + assert result.coord(crd).units == expected + + +def test_sin_phi(): + """Test that the function returns the values expected.""" + latitudes = np.array([50.0, 50.0, 50.0]) + longitudes = np.array([-5.0, 0.0, 5.0]) + dtval = datetime.datetime(2017, 1, 11, 8) + expected_array = np.array([-0.05481607, -0.00803911, 0.03659632]) + plugin = TemporalInterpolation( + interval_in_minutes=60, interpolation_method="solar" + ) + result = plugin.calc_sin_phi(dtval, latitudes, longitudes) + assert isinstance(result, np.ndarray) + np.testing.assert_almost_equal(result, expected_array) + + +@pytest.mark.parametrize("spatial_grid,expected_lats,expected_lons", + [ + ( + "latlon", + np.array([[40.0, 40.0, 40.0], [60.0, 60.0, 60.0], [80.0, 80.0, 80.0]]), + np.array([[-20.0, 0.0, 20.0], [-20.0, 0.0, 20.0], [-20.0, 0.0, 20.0]]), + ), + ( + "equalarea", + np.array( + [ + [53.84618597, 53.99730779, 53.93247526], + [56.82670954, 56.99111356, 56.9205672], + [59.8045105, 59.98499383, 59.90752513], + ] + ), + np.array( + [ + [-8.58580705, -3.51660018, 1.56242662], + [-9.06131306, -3.59656346, 1.88105082], + [-9.63368459, -3.69298822, 2.26497216], + ] + ) + ), + ] +) +def test_calc_lats_lons(spatial_grid, expected_lats, expected_lons): + """Test that the function returns the lats and lons expected for a native + lat-lon projection and for an equal areas projection.""" + + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 6, 9]] + data = np.ones((3, 3), dtype=np.float32) + cube = multi_time_cube(times, data, spatial_grid) + + plugin = TemporalInterpolation( + interval_in_minutes=60, interpolation_method="solar" + ) + result_lats, result_lons = plugin.calc_lats_lons(cube) + assert isinstance(result_lats, np.ndarray) + assert result_lats.shape == (3, 3) + assert isinstance(result_lons, np.ndarray) + assert result_lons.shape == (3, 3) + np.testing.assert_almost_equal(result_lats, expected_lats) + np.testing.assert_almost_equal(result_lons, expected_lons) + + +@pytest.mark.parametrize("realizations", (None, [0, 1, 2])) +def test_solar_interpolation(solar_expected, realizations): + """Test interpolation using the solar method. Apply it to deterministic + and ensemble data.""" + + times = [datetime.datetime(2017, 11, 1, hour) for hour in [6, 8, 10]] + npoints = 5 + data = np.stack( + [ + np.zeros((npoints, npoints), dtype=np.float32), + np.ones((npoints, npoints), dtype=np.float32), + np.ones((npoints, npoints), dtype=np.float32), + ] + ) + cube = multi_time_cube(times, data, "latlon", realizations=realizations) + interpolated_cube = cube[1].copy() + interpolated_cube = iris.util.new_axis(interpolated_cube, "time") + cube = cube[0::2] + + plugin = TemporalInterpolation(interpolation_method="solar", times=[times[1]]) + + result = plugin.solar_interpolate(cube, interpolated_cube) + assert isinstance(result, CubeList) + result, = result + assert result.coord("time").points == 1509523200 + assert result.coord("forecast_period").points[0] == 3600 * 4 + if result.ndim == 2: + np.testing.assert_almost_equal(result.data, solar_expected) + else: + for dslice in result.data: + np.testing.assert_almost_equal(dslice, solar_expected) + + +@pytest.mark.parametrize("realizations", (None, [0, 1, 2])) +def test_daynight_interpolation(daynight_mask, realizations): + """Test daynight function applies a suitable mask to interpolated + data. In this test the day-night terminator crosses the domain to + ensure the impact is captured. A deterministic and ensemble input + are tested.""" + + frt = datetime.datetime(2017, 11, 1, 6) + time = datetime.datetime(2017, 11, 1, 8) + data = np.ones((10, 10), dtype=np.float32) * 4 + interpolated_cube = diagnostic_cube(time, frt, data, "latlon", realizations=realizations) + + expected = np.where(daynight_mask == 0, 0, data) + + plugin = TemporalInterpolation(interpolation_method="daynight", times=[time]) + result = plugin.daynight_interpolate(interpolated_cube) + assert isinstance(result, CubeList) + + result, = result + assert result.coord("time").points == 1509523200 + assert result.coord("forecast_period").points[0] == 7200 + + if result.ndim == 2: + np.testing.assert_almost_equal(result.data, expected) + else: + for dslice in result.data: + np.testing.assert_almost_equal(dslice, expected) class Test_process(IrisTest): From 56c579a133f89a22cb4febca2380986b865ff1de Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Thu, 15 Feb 2024 09:11:19 +0000 Subject: [PATCH 06/19] All existing unit tests for temporal interpolation converted to pytest. --- .../utilities/test_TemporalInterpolation.py | 405 +++++++----------- 1 file changed, 147 insertions(+), 258 deletions(-) diff --git a/improver_tests/utilities/test_TemporalInterpolation.py b/improver_tests/utilities/test_TemporalInterpolation.py index 8e9ed258a5..f28f3bd520 100644 --- a/improver_tests/utilities/test_TemporalInterpolation.py +++ b/improver_tests/utilities/test_TemporalInterpolation.py @@ -53,7 +53,17 @@ def _grid_params(spatial_grid: str, npoints: int) -> Tuple[Tuple[float, float], float]: """Set domain corner and grid spacing for lat-lon or equal area - projections.""" + projections. + + Args: + spatial_grid: + "latlon" or "equalarea" to determine the type of projection. + npoints: + The number of grid points to use in both x and y. + Returns: + A tuple containing a further tuple that includes the grid corner + coordinates, and a single value specifying the grid spacing. + """ domain_corner = None grid_spacing = None @@ -78,6 +88,9 @@ def diagnostic_cube(time: dt, frt: dt, data: np.ndarray, spatial_grid: str, real The data to be contained in the cube. spatial_grid: Whether this is a lat-lon or equal areas projection. + realizations: + An optional list of realizations identifiers. The length of this + list will determine how many realizations are created. Returns: A diagnostic cube for use in testing. """ @@ -113,6 +126,9 @@ def multi_time_cube(times: List, data: np.ndarray, spatial_grid: str, bounds: bo Whether this is a lat-lon or equal areas projection. bounds: If True return time coordinates with time bounds. + realizations: + An optional list of realizations identifiers. The length of this + list will determine how many realizations are created. Returns: A diagnostic cube for use in testing. """ @@ -178,8 +194,7 @@ def solar_expected(): ) -@pytest.fixture -def daynight_mask(): +def mask_values(): """The mask matching the terminator position for the daynight interpolation tests.""" return np.array( @@ -198,6 +213,12 @@ def daynight_mask(): ) +@pytest.fixture +def daynight_mask(): + """Fixture to return mask values.""" + return mask_values() + + @pytest.mark.parametrize( "kwargs,exception", [ @@ -291,7 +312,9 @@ def test_enforce_time_coords_dtype(bounds): def test_sin_phi(): - """Test that the function returns the values expected.""" + """Test that the function returns the values expected for solar + elevation.""" + latitudes = np.array([50.0, 50.0, 50.0]) longitudes = np.array([-5.0, 0.0, 5.0]) dtval = datetime.datetime(2017, 1, 11, 8) @@ -412,268 +435,134 @@ def test_daynight_interpolation(daynight_mask, realizations): np.testing.assert_almost_equal(dslice, expected) -class Test_process(IrisTest): - - """Test interpolation of cubes to intermediate times using the plugin.""" - - def setUp(self): - """Set up the test inputs.""" - self.time_0 = datetime.datetime(2017, 11, 1, 3) - self.time_extra = datetime.datetime(2017, 11, 1, 6) - self.time_1 = datetime.datetime(2017, 11, 1, 9) - self.npoints = 10 - - domain_corner, grid_spacing = _grid_params("latlon", self.npoints) - - data_time_0 = np.ones((self.npoints, self.npoints), dtype=np.float32) - data_time_1 = np.ones((self.npoints, self.npoints), dtype=np.float32) * 7 - self.cube_time_0 = set_up_variable_cube( - data_time_0, - time=self.time_0, - frt=self.time_0, - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - self.cube_time_1 = set_up_variable_cube( - data_time_1, - time=self.time_1, - frt=self.time_0, - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - - def test_return_type(self): - """Test that an iris cubelist is returned.""" - - result = TemporalInterpolation(interval_in_minutes=180).process( - self.cube_time_0, self.cube_time_1 - ) - self.assertIsInstance(result, iris.cube.CubeList) - - def test_wind_direction_interpolation_over_north(self): - """Test that wind directions are interpolated properly at the 0/360 - circular cross-over.""" - - data_time_0 = np.ones((self.npoints, self.npoints), dtype=np.float32) * 350 - data_time_1 = np.ones((self.npoints, self.npoints), dtype=np.float32) * 20 - domain_corner, grid_spacing = _grid_params("latlon", self.npoints) - cube_time_0 = set_up_variable_cube( - data_time_0, - name="wind_from_direction", - units="degrees", - time=self.time_0, - frt=self.time_0, - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - cube_time_1 = set_up_variable_cube( - data_time_1, - name="wind_from_direction", - units="degrees", - time=self.time_1, - frt=self.time_1, - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - expected_data = np.full((self.npoints, self.npoints), 5, dtype=np.float32) - (result,) = TemporalInterpolation(interval_in_minutes=180).process( - cube_time_0, cube_time_1 - ) - - self.assertArrayAlmostEqual(expected_data, result.data, decimal=4) - - def test_wind_direction_interpolation(self): - """Test that wind directions are interpolated properly when the interpolation - doesn't cross the 0/360 boundary.""" - - data_time_0 = np.ones((self.npoints, self.npoints), dtype=np.float32) * 40 - data_time_1 = np.ones((self.npoints, self.npoints), dtype=np.float32) * 60 - domain_corner, grid_spacing = _grid_params("latlon", self.npoints) - cube_time_0 = set_up_variable_cube( - data_time_0, - units="degrees", - time=self.time_0, - frt=self.time_0, - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - cube_time_1 = set_up_variable_cube( - data_time_1, - units="degrees", - time=self.time_1, - frt=self.time_1, - domain_corner=domain_corner, - grid_spacing=grid_spacing, - ) - expected_data = expected_data = np.full( - (self.npoints, self.npoints), 50, dtype=np.float32 - ) - (result,) = TemporalInterpolation(interval_in_minutes=180).process( - cube_time_0, cube_time_1 - ) - - self.assertArrayAlmostEqual(expected_data, result.data, decimal=4) - - def test_valid_single_interpolation(self): - """Test interpolating to the mid point of the time range. Expect the - data to be half way between, and the time coordinate should be at - 06Z November 11th 2017.""" - - expected_data = np.ones((self.npoints, self.npoints)) * 4 - expected_time = [1509516000] - expected_fp = 3 * 3600 - (result,) = TemporalInterpolation(interval_in_minutes=180).process( - self.cube_time_0, self.cube_time_1 - ) - - self.assertArrayAlmostEqual(expected_data, result.data) - self.assertArrayAlmostEqual(result.coord("time").points, expected_time) - self.assertAlmostEqual(result.coord("forecast_period").points[0], expected_fp) - - def test_valid_multiple_interpolations(self): - """Test interpolating to every hour between the two input cubes. - Check the data increments as expected and the time coordinates are also - set correctly. - - NB Interpolation in iris is prone to float precision errors of order - 10E-6, hence the need to use AlmostEqual below.""" - - result = TemporalInterpolation(interval_in_minutes=60).process( - self.cube_time_0, self.cube_time_1 - ) - for i, cube in enumerate(result): - expected_data = np.ones((self.npoints, self.npoints)) * i + 2 - expected_time = [1509508800 + i * 3600] - - self.assertArrayAlmostEqual(expected_data, cube.data) - self.assertArrayAlmostEqual( - cube.coord("time").points, expected_time, decimal=5 - ) - self.assertAlmostEqual( - cube.coord("forecast_period").points[0], (i + 1) * 3600 - ) +@pytest.mark.parametrize("bearings,expected_value", [ + ([350, 20], 5), + ([40, 60], 50), + ] +) +def test_process_wind_direction(bearings, expected_value): + """Test that wind directions are interpolated properly at the 0/360 + circular cross-over and away from this cross-over. The interpolated + values are returned, along with the cube corresponding to the later + input time. This later cube should be the same as the input.""" + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] + npoints = 10 + data = np.stack( + [ + np.ones((npoints, npoints), dtype=np.float32) * bearings[0], + np.ones((npoints, npoints), dtype=np.float32) * bearings[1], + ] + ) + cube = multi_time_cube(times, data, "latlon") + cube.rename("wind_from_direction") + cube.units = "degrees" - def test_valid_interpolation_from_given_list(self): - """Test interpolating to a point defined in a list between the two - input cube validity times. Check the data increments as expected and - the time coordinates are also set correctly. - - NB Interpolation in iris is prone to float precision errors of order - 10E-6, hence the need to use AlmostEqual below.""" - - (result,) = TemporalInterpolation(times=[self.time_extra]).process( - self.cube_time_0, self.cube_time_1 - ) - expected_data = np.ones((self.npoints, self.npoints)) * 4 - expected_time = [1509516000] - expected_fp = 3 * 3600 - - self.assertArrayAlmostEqual(expected_data, result.data) - self.assertArrayAlmostEqual( - result.coord("time").points, expected_time, decimal=5 - ) - self.assertAlmostEqual(result.coord("forecast_period").points[0], expected_fp) - self.assertEqual(result.data.dtype, np.float32) - self.assertEqual(str(result.coord("time").points.dtype), "int64") - self.assertEqual(str(result.coord("forecast_period").points.dtype), "int32") - - def test_solar_interpolation_from_given_list(self): - """Test solar interpolating to a point defined in a list - between the two input cube validity times. - Check the data increments as expected and - the time coordinates are also set correctly.""" - - plugin = TemporalInterpolation( - times=[self.time_extra], interpolation_method="solar" - ) - (result,) = plugin.process(self.cube_time_0, self.cube_time_1) - expected_time = [1509516000] - expected_fp = 3 * 3600 - - self.assertArrayAlmostEqual( - result.coord("time").points, expected_time, decimal=5 - ) - self.assertAlmostEqual(result.coord("forecast_period").points[0], expected_fp) - self.assertEqual(result.data.dtype, np.float32) - self.assertEqual(str(result.coord("time").points.dtype), "int64") - self.assertEqual(str(result.coord("forecast_period").points.dtype), "int32") - - def test_daynight_interpolation_from_given_list(self): - """Test daynight interpolating to a point defined in a list - between the two input cube validity times. - Check the data increments as expected and - the time coordinates are also set correctly.""" - - plugin = TemporalInterpolation( - times=[self.time_extra], interpolation_method="daynight" - ) - (result,) = plugin.process(self.cube_time_0, self.cube_time_1) - expected_data = np.zeros((self.npoints, self.npoints)) - expected_data[:2, 7:] = 4.0 - expected_data[2, 8:] = 4.0 - expected_data[3, 9] = 4.0 - expected_time = [1509516000] - expected_fp = 3 * 3600 - - self.assertArrayAlmostEqual(expected_data, result.data) - self.assertArrayAlmostEqual( - result.coord("time").points, expected_time, decimal=5 - ) - self.assertAlmostEqual(result.coord("forecast_period").points[0], expected_fp) - self.assertEqual(result.data.dtype, np.float32) - self.assertEqual(str(result.coord("time").points.dtype), "int64") - self.assertEqual(str(result.coord("forecast_period").points.dtype), "int32") - - def test_input_cube_without_time_coordinate(self): - """Test that an exception is raised if a cube is provided without a - time coordiate.""" - - self.cube_time_0.remove_coord("time") - - msg = "Cube provided to TemporalInterpolation " "contains no time coordinate" - with self.assertRaisesRegex(CoordinateNotFoundError, msg): - TemporalInterpolation(interval_in_minutes=180).process( - self.cube_time_0, self.cube_time_1 - ) + expected = np.full((npoints, npoints), expected_value, dtype=np.float32) + result = TemporalInterpolation(interval_in_minutes=180).process( + cube[0], cube[1] + ) - def test_input_cubes_in_incorrect_time_order(self): - """Test that an exception is raised if the cube representing the - initial time has a validity time that is after the cube representing - the final time.""" + assert isinstance(result, CubeList) + np.testing.assert_almost_equal(result[0].data, expected, decimal=4) + # Ideally the data at the later of the input cube times would be + # completely unchanged by the interpolation, but there are some small + # precision level changes. + np.testing.assert_almost_equal(result[1].data, cube[1].data, decimal=5) - msg = "TemporalInterpolation input cubes ordered incorrectly" - with self.assertRaisesRegex(ValueError, msg): - TemporalInterpolation(interval_in_minutes=180).process( - self.cube_time_1, self.cube_time_0 - ) - def test_input_cube_with_multiple_times(self): - """Test that an exception is raised if a cube is provided that has - multiple validity times, e.g. a multi-entried time dimension.""" - second_time = self.cube_time_0.copy() - second_time.coord("time").points = self.time_extra.timestamp() - cube = iris.cube.CubeList([self.cube_time_0, second_time]) - cube = cube.merge_cube() +@pytest.mark.parametrize( + "kwargs,offsets,expected", + [ + ({"interval_in_minutes": 180}, [3, 6], [4, 7]), + ({"interval_in_minutes": 60}, [1, 2, 3, 4, 5, 6], [2, 3, 4, 5, 6, 7]), + ({"times": [datetime.datetime(2017, 11, 1, 6)]}, [3, 6], [4, 7]), + ({"times": [datetime.datetime(2017, 11, 1, 8)], "interpolation_method": "daynight"}, [5], [6 * mask_values()]), + ({"times": [datetime.datetime(2017, 11, 1, 8)], "interpolation_method": "solar"}, [5], [None]) + ] +) +def test_process_interpolation(kwargs, offsets, expected): + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] + npoints = 10 + data = np.stack( + [ + np.ones((npoints, npoints), dtype=np.float32), + np.ones((npoints, npoints), dtype=np.float32) * 7, + ] + ) + cube = multi_time_cube(times, data, "latlon") - msg = "Cube provided to TemporalInterpolation contains multiple" - with self.assertRaisesRegex(ValueError, msg): - TemporalInterpolation(interval_in_minutes=180).process( - cube, self.cube_time_1 - ) + result = TemporalInterpolation(**kwargs).process(cube[0], cube[1]) - def test_input_cubelists_raises_exception(self): - """Test that providing cubelists instead of cubes raises an - exception.""" + for i, (offset, value) in enumerate(zip(offsets, expected)): + expected_data = np.full((npoints, npoints), value) + expected_time = 1509505200 + (offset * 3600) + expected_fp = (6 + offset) * 3600 - cubes = iris.cube.CubeList([self.cube_time_1]) + assert result[i].coord("time").points[0] == expected_time + assert result[i].coord("forecast_period").points[0] == expected_fp + assert result[i].coord("time").points.dtype == "int64" + assert result[i].coord("forecast_period").points.dtype == "int32" + if value is not None: + np.testing.assert_almost_equal(result[i].data, expected_data) - msg = "Inputs to TemporalInterpolation are not of type " - with self.assertRaisesRegex(TypeError, msg): - TemporalInterpolation(interval_in_minutes=180).process( - cubes, self.cube_time_0 - ) +def test_input_cube_without_time_coordinate(): + """Test that an exception is raised if a cube is provided without a + time coordiate.""" + + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] + data = np.ones((5, 5), dtype=np.float32) + cube = multi_time_cube(times, data, "latlon") + cube_0 = cube[0] + cube_1 = cube[1] + cube_1.remove_coord("time") + + msg = "Cube provided to TemporalInterpolation contains no time coordinate" + with pytest.raises(CoordinateNotFoundError, match=msg): + TemporalInterpolation(interval_in_minutes=180).process(cube_0, cube_1) + + +def test_input_cubes_in_incorrect_time_order(): + """Test that an exception is raised if the cube representing the + initial time has a validity time that is after the cube representing + the final time.""" + + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] + data = np.ones((5, 5), dtype=np.float32) + cube = multi_time_cube(times, data, "latlon") + cube_0 = cube[0] + cube_1 = cube[1] + + msg = "TemporalInterpolation input cubes ordered incorrectly" + with pytest.raises(ValueError, match=msg): + TemporalInterpolation(interval_in_minutes=180).process(cube_1, cube_0) + + +def test_input_cube_with_multiple_times(): + """Test that an exception is raised if a cube is provided that has + multiple validity times, e.g. a multi-entried time dimension.""" + + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 6, 9]] + data = np.ones((5, 5), dtype=np.float32) + cube = multi_time_cube(times, data, "latlon") + cube_0 = cube[0:2] + cube_1 = cube[2] + + msg = "Cube provided to TemporalInterpolation contains multiple" + with pytest.raises(ValueError, match=msg): + TemporalInterpolation(interval_in_minutes=60).process(cube_0, cube_1) + + +def test_input_cubelists_raises_exception(): + """Test that providing cubelists instead of cubes raises an + exception.""" + + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] + data = np.ones((5, 5), dtype=np.float32) + cube = multi_time_cube(times, data, "latlon") + cubes = CubeList([cube[0], cube[1]]) -if __name__ == "__main__": - unittest.main() + msg = "Inputs to TemporalInterpolation are not of type " + with pytest.raises(TypeError, match=msg): + TemporalInterpolation(interval_in_minutes=180).process(cubes, cube[1]) From 1b7a0cb4ddc60becacbac0caeb28d02a220e5ea5 Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Thu, 15 Feb 2024 09:17:40 +0000 Subject: [PATCH 07/19] Style fixes. --- improver/utilities/temporal_interpolation.py | 9 +- .../utilities/test_TemporalInterpolation.py | 124 +++++++++++------- 2 files changed, 82 insertions(+), 51 deletions(-) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index 7e13eef2cd..41ed070c93 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -73,7 +73,6 @@ def __init__( times: Optional[List[datetime]] = None, interpolation_method: str = "linear", accumulation: bool = False, - ) -> None: """ Initialise class. @@ -403,7 +402,7 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: raise TypeError(msg) if cube_t0.coord("time").bounds is not None: - if not self.interpolation_method in self.period_interpolation_methods: + if self.interpolation_method not in self.period_interpolation_methods: raise ValueError( "Period diagnostics can only be temporally interpolated " f"using these methods: {self.period_interpolation_methods}.\n" @@ -415,9 +414,7 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: (initial_time,) = iris_time_to_datetime(cube_t0.coord("time")) (final_time,) = iris_time_to_datetime(cube_t1.coord("time")) except CoordinateNotFoundError: - msg = ( - "Cube provided to TemporalInterpolation contains no time coordinate." - ) + msg = "Cube provided to TemporalInterpolation contains no time coordinate." raise CoordinateNotFoundError(msg) except ValueError: msg = ( @@ -478,7 +475,7 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: # allows for trends in the data to add some variation to the # different periods that are created. if self.accumulation: - time_coord, = interpolated_cube.coord_dims("time") + (time_coord,) = interpolated_cube.coord_dims("time") interpolated_total = np.sum(interpolated_cube.data, axis=time_coord) renormalisation = cube_t1.data / interpolated_total interpolated_cube.data *= renormalisation diff --git a/improver_tests/utilities/test_TemporalInterpolation.py b/improver_tests/utilities/test_TemporalInterpolation.py index f28f3bd520..42c1f0e93b 100644 --- a/improver_tests/utilities/test_TemporalInterpolation.py +++ b/improver_tests/utilities/test_TemporalInterpolation.py @@ -32,22 +32,15 @@ import datetime from datetime import datetime as dt -import unittest - -from typing import List, Tuple, Optional +from typing import List, Optional, Tuple import iris import numpy as np import pytest - from iris.cube import Cube, CubeList from iris.exceptions import CoordinateNotFoundError -from iris.tests import IrisTest -from improver.synthetic_data.set_up_test_cubes import ( - add_coordinate, - set_up_variable_cube, -) +from improver.synthetic_data.set_up_test_cubes import set_up_variable_cube from improver.utilities.temporal_interpolation import TemporalInterpolation @@ -76,7 +69,13 @@ def _grid_params(spatial_grid: str, npoints: int) -> Tuple[Tuple[float, float], return domain_corner, grid_spacing -def diagnostic_cube(time: dt, frt: dt, data: np.ndarray, spatial_grid: str, realizations: Optional[List] = None) -> Cube: +def diagnostic_cube( + time: dt, + frt: dt, + data: np.ndarray, + spatial_grid: str, + realizations: Optional[List] = None, +) -> Cube: """Return a diagnostic cube containing the provided data. Args: @@ -111,7 +110,13 @@ def diagnostic_cube(time: dt, frt: dt, data: np.ndarray, spatial_grid: str, real ) -def multi_time_cube(times: List, data: np.ndarray, spatial_grid: str, bounds: bool = False, realizations: Optional[List] = None) -> Cube: +def multi_time_cube( + times: List, + data: np.ndarray, + spatial_grid: str, + bounds: bool = False, + realizations: Optional[List] = None, +) -> Cube: """Return a multi-time diagnostic cube containing the provided data. Args: @@ -138,7 +143,11 @@ def multi_time_cube(times: List, data: np.ndarray, spatial_grid: str, bounds: bo frt = sorted(times)[0] - (times[1] - times[0]) # Such that guess bounds is +ve for time, data_slice in zip(times, data): - cubes.append(diagnostic_cube(time, frt, data_slice, spatial_grid, realizations=realizations)) + cubes.append( + diagnostic_cube( + time, frt, data_slice, spatial_grid, realizations=realizations + ) + ) cube = cubes.merge_cube() @@ -148,7 +157,9 @@ def multi_time_cube(times: List, data: np.ndarray, spatial_grid: str, bounds: bo return cube -def non_standard_times(times: List, data: np.ndarray, spatial_grid: str, bounds: bool = False) -> Cube: +def non_standard_times( + times: List, data: np.ndarray, spatial_grid: str, bounds: bool = False +) -> Cube: """Return a multi-time diagnostic cube containing the provided data. The units of the time dimensions are made non-standards compliant. @@ -223,9 +234,15 @@ def daynight_mask(): "kwargs,exception", [ ({}, "TemporalInterpolation: One of"), # No target times defined - ({"interval_in_minutes":60, "times":[datetime.datetime(2017, 11, 1, 9)]}, "TemporalInterpolation: Only one of"), # Two methods of defining targets used - ({"interval_in_minutes":60, "interpolation_method":"invalid"}, "TemporalInterpolation: Unknown interpolation method"), # Invalid interpolation method requested - ] + ( + {"interval_in_minutes": 60, "times": [datetime.datetime(2017, 11, 1, 9)]}, + "TemporalInterpolation: Only one of", + ), # Two methods of defining targets used + ( + {"interval_in_minutes": 60, "interpolation_method": "invalid"}, + "TemporalInterpolation: Unknown interpolation method", + ), # Invalid interpolation method requested + ], ) def test__init__(kwargs, exception): """Test exceptions raised by the __init__ method.""" @@ -236,11 +253,20 @@ def test__init__(kwargs, exception): @pytest.mark.parametrize( "kwargs,exception", [ - ({"interval_in_minutes": 60}, None), # Generate times between bounds using interval + ( + {"interval_in_minutes": 60}, + None, + ), # Generate times between bounds using interval ({"times": None}, None), # Use the expected times as the input - ({"times": datetime.datetime(2017, 11, 1, 10)}, "List of times falls outside the range given by"), # Use the expected times, plus another outside the range as the input - ({"interval_in_minutes":61}, "interval_in_minutes of"), # Use an invalid interval - ] + ( + {"times": datetime.datetime(2017, 11, 1, 10)}, + "List of times falls outside the range given by", + ), # Use the expected times, plus another outside the range as the input + ( + {"interval_in_minutes": 61}, + "interval_in_minutes of", + ), # Use an invalid interval + ], ) def test_construct_time_list(kwargs, exception): """Test construction of target times using various inputs and testing @@ -258,7 +284,9 @@ def test_construct_time_list(kwargs, exception): # times plus any others specified in the kwarg. try: target_times = times.copy() - target_times.append(kwargs["times"]) if kwargs["times"] is not None else target_times + target_times.append(kwargs["times"]) if kwargs[ + "times" + ] is not None else target_times except KeyError: pass else: @@ -319,15 +347,14 @@ def test_sin_phi(): longitudes = np.array([-5.0, 0.0, 5.0]) dtval = datetime.datetime(2017, 1, 11, 8) expected_array = np.array([-0.05481607, -0.00803911, 0.03659632]) - plugin = TemporalInterpolation( - interval_in_minutes=60, interpolation_method="solar" - ) + plugin = TemporalInterpolation(interval_in_minutes=60, interpolation_method="solar") result = plugin.calc_sin_phi(dtval, latitudes, longitudes) assert isinstance(result, np.ndarray) np.testing.assert_almost_equal(result, expected_array) -@pytest.mark.parametrize("spatial_grid,expected_lats,expected_lons", +@pytest.mark.parametrize( + "spatial_grid,expected_lats,expected_lons", [ ( "latlon", @@ -349,9 +376,9 @@ def test_sin_phi(): [-9.06131306, -3.59656346, 1.88105082], [-9.63368459, -3.69298822, 2.26497216], ] - ) + ), ), - ] + ], ) def test_calc_lats_lons(spatial_grid, expected_lats, expected_lons): """Test that the function returns the lats and lons expected for a native @@ -361,9 +388,7 @@ def test_calc_lats_lons(spatial_grid, expected_lats, expected_lons): data = np.ones((3, 3), dtype=np.float32) cube = multi_time_cube(times, data, spatial_grid) - plugin = TemporalInterpolation( - interval_in_minutes=60, interpolation_method="solar" - ) + plugin = TemporalInterpolation(interval_in_minutes=60, interpolation_method="solar") result_lats, result_lons = plugin.calc_lats_lons(cube) assert isinstance(result_lats, np.ndarray) assert result_lats.shape == (3, 3) @@ -396,7 +421,7 @@ def test_solar_interpolation(solar_expected, realizations): result = plugin.solar_interpolate(cube, interpolated_cube) assert isinstance(result, CubeList) - result, = result + (result,) = result assert result.coord("time").points == 1509523200 assert result.coord("forecast_period").points[0] == 3600 * 4 if result.ndim == 2: @@ -416,7 +441,9 @@ def test_daynight_interpolation(daynight_mask, realizations): frt = datetime.datetime(2017, 11, 1, 6) time = datetime.datetime(2017, 11, 1, 8) data = np.ones((10, 10), dtype=np.float32) * 4 - interpolated_cube = diagnostic_cube(time, frt, data, "latlon", realizations=realizations) + interpolated_cube = diagnostic_cube( + time, frt, data, "latlon", realizations=realizations + ) expected = np.where(daynight_mask == 0, 0, data) @@ -424,7 +451,7 @@ def test_daynight_interpolation(daynight_mask, realizations): result = plugin.daynight_interpolate(interpolated_cube) assert isinstance(result, CubeList) - result, = result + (result,) = result assert result.coord("time").points == 1509523200 assert result.coord("forecast_period").points[0] == 7200 @@ -435,11 +462,7 @@ def test_daynight_interpolation(daynight_mask, realizations): np.testing.assert_almost_equal(dslice, expected) -@pytest.mark.parametrize("bearings,expected_value", [ - ([350, 20], 5), - ([40, 60], 50), - ] -) +@pytest.mark.parametrize("bearings,expected_value", [([350, 20], 5), ([40, 60], 50)]) def test_process_wind_direction(bearings, expected_value): """Test that wind directions are interpolated properly at the 0/360 circular cross-over and away from this cross-over. The interpolated @@ -458,9 +481,7 @@ def test_process_wind_direction(bearings, expected_value): cube.units = "degrees" expected = np.full((npoints, npoints), expected_value, dtype=np.float32) - result = TemporalInterpolation(interval_in_minutes=180).process( - cube[0], cube[1] - ) + result = TemporalInterpolation(interval_in_minutes=180).process(cube[0], cube[1]) assert isinstance(result, CubeList) np.testing.assert_almost_equal(result[0].data, expected, decimal=4) @@ -470,16 +491,29 @@ def test_process_wind_direction(bearings, expected_value): np.testing.assert_almost_equal(result[1].data, cube[1].data, decimal=5) - @pytest.mark.parametrize( "kwargs,offsets,expected", [ ({"interval_in_minutes": 180}, [3, 6], [4, 7]), ({"interval_in_minutes": 60}, [1, 2, 3, 4, 5, 6], [2, 3, 4, 5, 6, 7]), ({"times": [datetime.datetime(2017, 11, 1, 6)]}, [3, 6], [4, 7]), - ({"times": [datetime.datetime(2017, 11, 1, 8)], "interpolation_method": "daynight"}, [5], [6 * mask_values()]), - ({"times": [datetime.datetime(2017, 11, 1, 8)], "interpolation_method": "solar"}, [5], [None]) - ] + ( + { + "times": [datetime.datetime(2017, 11, 1, 8)], + "interpolation_method": "daynight", + }, + [5], + [6 * mask_values()], + ), + ( + { + "times": [datetime.datetime(2017, 11, 1, 8)], + "interpolation_method": "solar", + }, + [5], + [None], + ), + ], ) def test_process_interpolation(kwargs, offsets, expected): times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] From 35dae5153f71f9c25c424b0516abe49e5523d235 Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Fri, 16 Feb 2024 11:55:51 +0000 Subject: [PATCH 08/19] Unit tests for new functionality added. Fixed a bug in the application of the daynight method to cubes that have been interpolated to multiple times, allowing for other coordinates as well. Style fixes. --- improver/cli/temporal_interpolate.py | 24 +- improver/utilities/temporal_interpolation.py | 213 +++++++++++-- .../utilities/test_TemporalInterpolation.py | 296 ++++++++++++++++-- 3 files changed, 472 insertions(+), 61 deletions(-) diff --git a/improver/cli/temporal_interpolate.py b/improver/cli/temporal_interpolate.py index 70e3abd3fd..3dc27087c4 100755 --- a/improver/cli/temporal_interpolate.py +++ b/improver/cli/temporal_interpolate.py @@ -43,6 +43,9 @@ def process( interval_in_mins: int = None, times: cli.comma_separated_list = None, interpolation_method="linear", + accumulation: bool = False, + maximum: bool = False, + minimum: bool = False, ): """Interpolate data between validity times. @@ -75,7 +78,23 @@ def process( solar interpolates using the solar elevation, daynight uses linear interpolation but sets night time points to 0.0 linear is linear interpolation. - + accumulation: + Set True if the diagnostic being temporally interpolated is a + period accumulation. The output will be renormalised to ensure + that the total across the period constructed from the shorter + intervals matches the total across the period from the coarser + intervals. Trends between adjacent input periods will be used + to provide variation across the interpolated periods. + maximum: + Set True if the diagnostic being temporally interpolated is a + period maximum. Trends between adjacent input periods will be used + to provide variation across the interpolated periods where these + are consistent with the inputs. + minimum: + Set True if the diagnostic being temporally interpolated is a + period minimum. Trends between adjacent input periods will be used + to provide variation across the interpolated periods where these + are consistent with the inputs. Returns: iris.cube.CubeList: A list of cubes interpolated to the desired times. The @@ -99,5 +118,8 @@ def process( interval_in_minutes=interval_in_mins, times=times, interpolation_method=interpolation_method, + accumulation=accumulation, + maximum=maximum, + minimum=minimum, )(start_cube, end_cube) return MergeCubes()(result) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index 41ed070c93..eae0df91f6 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -59,12 +59,37 @@ class TemporalInterpolation(BasePlugin): The plugin will return the interpolated times and the later of the two input times. This allows us to modify the input diagnostics if they - represent periods. The IMPROVER convention is that period diagnostics - have their time coordinate point at the end of the period. The later of - the two inputs therefore covers the period that has been broken down into - shorter periods by the interpolation and must itself be modified. The - result of this approach is that in a long run of lead-times, e.g. T+0 to - T+120 all the lead-times will be available except T+0. + represent accumulations. + + The IMPROVER convention is that period diagnostics have their time + coordinate point at the end of the period. The later of the two inputs + therefore covers the period that has been broken down into shorter periods + by the interpolation and, if working with accumulations, must itself be + modified. The result of this approach is that in a long run of + lead-times, e.g. T+0 to T+120 all the lead-times will be available except + T+0. + + If working with period maximums and minimums we cannot return values in + the new periods that do not adhere to the inputs. For example, we might + have a 3-hour maximum of 5 ms-1. The period before it might have a + maximum of 10 ms-1. However, on splitting the 3-hour period into 1-hour + periods, we cannot use the gradient between the two original period + maximums to produce a wind speed maximum in any of the interpolated + e.g. 1-hour periods that is above 5 ms-1 as we already know this to be the + maximum over the full 3-hours. We could use the trend to produce a lower + maximum in the interpolated 1-hour periods as this is not an unreasonable + interpretation of the gradient information and is consistent with the + original period maximum. A similar argument can be made about modifying + the minimums. + + We could use the cell methods to determine whether we are working with + accumulations, maximums, or minimums. This should be denoted as a cell + method associated with the time coordinate, e.g. for an accumulation it + would be `time: sum`, whilst a maximum would have `time: max`. However + we cannot guarantee these cell methods are present and not assumed. As + such the interpolation of periods here relies on the user supplying a + suitable keyword argument that denotes the type of period being + processed. """ def __init__( @@ -73,6 +98,8 @@ def __init__( times: Optional[List[datetime]] = None, interpolation_method: str = "linear", accumulation: bool = False, + maximum: bool = False, + minimum: bool = False, ) -> None: """ Initialise class. @@ -98,10 +125,22 @@ def __init__( that the total across the period constructed from the shorter intervals matches the total across the period from the coarser intervals. + maximum: + Set True if the diagnostic being temporally interpolated is a + period maximum. Trends between adjacent input periods will be used + to provide variation across the interpolated periods where these + are consistent with the inputs. + minimum: + Set True if the diagnostic being temporally interpolated is a + period minimum. Trends between adjacent input periods will be used + to provide variation across the interpolated periods where these + are consistent with the inputs. Raises: ValueError: If neither interval_in_minutes nor times are set. + ValueError: If both interval_in_minutes and times are both set. ValueError: If interpolation method not in known list. + ValueError: If multiple period diagnostic kwargs are set True. """ if interval_in_minutes is None and times is None: raise ValueError( @@ -126,7 +165,15 @@ def __init__( ) self.interpolation_method = interpolation_method self.period_inputs = False + if np.sum([accumulation, maximum, minimum]) > 1: + raise ValueError( + "Only one type of period diagnostics may be specified: " + f"accumulation = {accumulation}, maximum = {maximum}, " + f"minimum = {minimum}" + ) self.accumulation = accumulation + self.maximum = maximum + self.minimum = minimum def construct_time_list( self, initial_time: datetime, final_time: datetime @@ -364,9 +411,107 @@ def daynight_interpolate(interpolated_cube: Cube) -> CubeList: daynightplugin = DayNightMask() daynight_mask = daynightplugin(interpolated_cube) index = daynight_mask.data == daynightplugin.night - interpolated_cube.data[..., index] = 0.0 + + # Reshape the time, y, x mask to match the input which may include addtional + # dimensions, such as realization. + dropped_crds = [ + crd + for crd in interpolated_cube.coords(dim_coords=True) + if crd not in daynight_mask.coords(dim_coords=True) + ] + if dropped_crds: + cslices = interpolated_cube.slices_over(dropped_crds) + masked_data = CubeList() + for cslice in cslices: + cslice.data[index] = 0.0 + masked_data.append(cslice) + interpolated_cube = masked_data.merge_cube() + else: + interpolated_cube.data[index] = 0.0 + return iris.cube.CubeList(list(interpolated_cube.slices_over("time"))) + @staticmethod + def add_bounds(cube_t0: Cube, interpolated_cube: Cube): + """Calcualte bounds using the interpolated times and the time + taken from cube_t0. This function is used rather than iris's guess + bounds method as we want to use the earlier time cube to inform + the lowest bound. The interpolated_cube `crd` is modified in + place. + + Args: + cube_t0: + The input cube corresponding to the earlier time. + interpolated_cube: + The cube containing the interpolated times, which includes + the data corresponding to the time of the later of the two + input cubes. + + Raises: + CoordinateNotFoundError: if time or forecast_period coordinates + are not present on the input cubes. + """ + for crd in ["time", "forecast_period"]: + try: + interpolated_times = np.concatenate( + [cube_t0.coord(crd).points, interpolated_cube.coord(crd).points] + ) + except CoordinateNotFoundError: + raise CoordinateNotFoundError( + f"Period diagnostic cube is missing expected coordinate: {crd}" + ) + bounds = np.diff(interpolated_times) + all_bounds = [] + for bound, time in zip(bounds, interpolated_times[1:]): + all_bounds.append([time - bound, time]) + interpolated_cube.coord(crd).bounds = all_bounds + + @staticmethod + def _calculate_accumulation( + cube_t0: Cube, period_reference: Cube, interpolated_cube: Cube + ): + """If the input is an accumulation we use the trapezium rule to + calculate a new accumulation for each output period from the rates + we converted the accumulations to prior to interpolating. We then + renormalise to ensure the total accumulation across the period is + unchanged by expressing it as a series of shorter periods. + + The interpolated cube is modified in place. + + Args: + cube_t0: + The input cube corresponding to the earlier time. + period_refrence: + The input cube corresponding to the later time, with the + values prior to conversion to rates. + interpolated_cube: + The cube containing the interpolated times, which includes + the data corresponding to the time of the later of the two + input cubes. + """ + # Calculate an average rate for the period from the edges + accumulation_edges = [cube_t0, *interpolated_cube.slices_over("time")] + period_rates = np.array( + [ + (a.data + b.data) / 2 + for a, b in zip(accumulation_edges[:-1], accumulation_edges[1:]) + ] + ) + interpolated_cube.data = period_rates + + # Multiply the averate rate by the length of each period to get a new + # accumulation. + new_periods = np.diff(interpolated_cube.coord("forecast_period").bounds) + new_periods = new_periods[:, np.newaxis] + interpolated_cube.data = np.multiply(new_periods, interpolated_cube.data) + + # Renormalise the total of the new periods to ensure it matches the + # total expressed in the longer original period. + (time_coord,) = interpolated_cube.coord_dims("time") + interpolated_total = np.sum(interpolated_cube.data, axis=time_coord) + renormalisation = period_reference.data / interpolated_total + interpolated_cube.data *= renormalisation + def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: """ Interpolate data to intermediate times between validity times of @@ -385,6 +530,10 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: Raises: TypeError: If cube_t0 and cube_t1 are not of type iris.cube.Cube. + ValueError: A period diagnostic is being interpolated with a method + not found in the period_interpolation_methods list. + ValueError: A period diagnostic is being processed without the + type of period being specified. CoordinateNotFoundError: The input cubes contain no time coordinate. ValueError: Cubes contain multiple validity times. @@ -408,6 +557,12 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: f"using these methods: {self.period_interpolation_methods}.\n" f"Currently selected method is: {self.interpolation_method}." ) + if not any([self.accumulation, self.maximum, self.minimum]): + raise ValueError( + "A type of period must be specified when interpolating a " + "period diagnostic. This may be a period maximum, minimum" + " or an accumulation." + ) self.period_inputs = True try: @@ -441,6 +596,13 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: cube_t0.data = WindDirection.deg_to_complex(cube_t0.data) cube_t1.data = WindDirection.deg_to_complex(cube_t1.data) + # Convert accumulations into rates to allow interpolation using trends + # in the data and to accommodate non-uniform output intervals. + if self.accumulation: + cube_t0.data /= np.diff(cube_t0.coord("forecast_period").bounds[0])[0] + period_reference = cube_t1.copy() + cube_t1.data /= np.diff(cube_t1.coord("forecast_period").bounds[0])[0] + cubes = iris.cube.CubeList([cube_t0, cube_t1]) cube = MergeCubes()(cubes) @@ -453,32 +615,21 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: if self.period_inputs: # Add bounds to the time coordinates of the interpolated outputs # if the inputs were period diagnostics. - for crd in ["time", "forecast_period"]: - interpolated_cube.coord(crd).guess_bounds(bound_position=1.0) - - # Check the interpolated output bounds cover the expected period. - lowest_bound = interpolated_cube.coord(crd).bounds[0][0] - highest_bound = interpolated_cube.coord(crd).bounds[-1][-1] - lowest_match = lowest_bound == cube_t1.coord(crd).bounds[0][0] - highest_match = highest_bound == cube_t1.coord(crd).bounds[-1][-1] - if not lowest_match or not highest_match: - raise ValueError( - "The interpolated periods do not cover the entire " - "period of the input." - ) + self.add_bounds(cube_t0, interpolated_cube) - # If the input is an accumulation the total must be renormalised - # to avoid double counting. The input cube the contains an - # accumulation spanning the whole period is adjusted down to - # represent only a suitable fraction of the total time. The - # use of interpolation, rather than just splitting up the period, - # allows for trends in the data to add some variation to the - # different periods that are created. + # Apply suitable if self.accumulation: - (time_coord,) = interpolated_cube.coord_dims("time") - interpolated_total = np.sum(interpolated_cube.data, axis=time_coord) - renormalisation = cube_t1.data / interpolated_total - interpolated_cube.data *= renormalisation + self._calculate_accumulation( + cube_t0, period_reference, interpolated_cube + ) + elif self.maximum: + interpolated_cube.data = np.minimum( + cube_t1.data, interpolated_cube.data + ) + elif self.minimum: + interpolated_cube.data = np.maximum( + cube_t1.data, interpolated_cube.data + ) self.enforce_time_coords_dtype(interpolated_cube) interpolated_cubes = iris.cube.CubeList() diff --git a/improver_tests/utilities/test_TemporalInterpolation.py b/improver_tests/utilities/test_TemporalInterpolation.py index 42c1f0e93b..cb0eed9c29 100644 --- a/improver_tests/utilities/test_TemporalInterpolation.py +++ b/improver_tests/utilities/test_TemporalInterpolation.py @@ -141,7 +141,7 @@ def multi_time_cube( if data.ndim == 2: data = np.stack([data] * len(times)) - frt = sorted(times)[0] - (times[1] - times[0]) # Such that guess bounds is +ve + frt = sorted(times)[0] - (times[1] - times[0]) # Such that guess bounds are +ve for time, data_slice in zip(times, data): cubes.append( diagnostic_cube( @@ -210,16 +210,30 @@ def mask_values(): interpolation tests.""" return np.array( [ - [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], - [0, 1, 1, 1, 1, 1, 1, 1, 1, 1], - [0, 1, 1, 1, 1, 1, 1, 1, 1, 1], - [0, 0, 1, 1, 1, 1, 1, 1, 1, 1], - [0, 0, 0, 1, 1, 1, 1, 1, 1, 1], - [0, 0, 0, 0, 1, 1, 1, 1, 1, 1], - [0, 0, 0, 0, 0, 0, 1, 1, 1, 1], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 1], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [0, 0, 1, 1, 1, 1, 1, 1, 1, 1], + [0, 0, 0, 1, 1, 1, 1, 1, 1, 1], + [0, 0, 0, 0, 1, 1, 1, 1, 1, 1], + [0, 0, 0, 0, 0, 0, 1, 1, 1, 1], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ], + [ + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [0, 0, 1, 1, 1, 1, 1, 1, 1, 1], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ], ] ) @@ -242,6 +256,10 @@ def daynight_mask(): {"interval_in_minutes": 60, "interpolation_method": "invalid"}, "TemporalInterpolation: Unknown interpolation method", ), # Invalid interpolation method requested + ( + {"interval_in_minutes": 60, "maximum": True, "minimum": True}, + "Only one type of period diagnostics may be specified:", + ), # Invalid interpolation method requested ], ) def test__init__(kwargs, exception): @@ -432,34 +450,49 @@ def test_solar_interpolation(solar_expected, realizations): @pytest.mark.parametrize("realizations", (None, [0, 1, 2])) -def test_daynight_interpolation(daynight_mask, realizations): +@pytest.mark.parametrize( + "interp_times,expected_times,expected_fps", + [ + ([8], [1509523200], [7200]), # Interpolate to a single time. + ([8, 10], [1509523200, 1509530400], [7200, 14400]), # Interpolate to two times. + ], +) +def test_daynight_interpolation( + daynight_mask, realizations, interp_times, expected_times, expected_fps +): """Test daynight function applies a suitable mask to interpolated data. In this test the day-night terminator crosses the domain to ensure the impact is captured. A deterministic and ensemble input are tested.""" - frt = datetime.datetime(2017, 11, 1, 6) - time = datetime.datetime(2017, 11, 1, 8) + times = [datetime.datetime(2017, 11, 1, hour) for hour in interp_times] data = np.ones((10, 10), dtype=np.float32) * 4 - interpolated_cube = diagnostic_cube( - time, frt, data, "latlon", realizations=realizations - ) - - expected = np.where(daynight_mask == 0, 0, data) + if len(times) > 1: + interpolated_cube = multi_time_cube( + times, data, "latlon", realizations=realizations + ) + else: + frt = datetime.datetime(2017, 11, 1, 6) + interpolated_cube = diagnostic_cube( + times, frt, data, "latlon", realizations=realizations + ) - plugin = TemporalInterpolation(interpolation_method="daynight", times=[time]) + plugin = TemporalInterpolation(interpolation_method="daynight", times=[times]) result = plugin.daynight_interpolate(interpolated_cube) + assert isinstance(result, CubeList) + for index, cube in enumerate(result): + expected = np.where(daynight_mask[index] == 0, 0, data) - (result,) = result - assert result.coord("time").points == 1509523200 - assert result.coord("forecast_period").points[0] == 7200 + assert cube.coord("time").points == expected_times[index] + assert cube.coord("forecast_period").points[0] == expected_fps[index] - if result.ndim == 2: - np.testing.assert_almost_equal(result.data, expected) - else: - for dslice in result.data: - np.testing.assert_almost_equal(dslice, expected) + if cube.coords("realization"): + cslices = cube.slices_over("realization") + else: + cslices = [cube] + for cslice in cslices: + np.testing.assert_almost_equal(cslice.data, expected) @pytest.mark.parametrize("bearings,expected_value", [([350, 20], 5), ([40, 60], 50)]) @@ -503,7 +536,7 @@ def test_process_wind_direction(bearings, expected_value): "interpolation_method": "daynight", }, [5], - [6 * mask_values()], + [6 * mask_values()[0]], ), ( { @@ -516,6 +549,10 @@ def test_process_wind_direction(bearings, expected_value): ], ) def test_process_interpolation(kwargs, offsets, expected): + """Test the process method with a variety of kwargs, selecting different + interpolation methods and output times. Check that the returned times and + data are as expected.""" + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] npoints = 10 data = np.stack( @@ -600,3 +637,204 @@ def test_input_cubelists_raises_exception(): msg = "Inputs to TemporalInterpolation are not of type " with pytest.raises(TypeError, match=msg): TemporalInterpolation(interval_in_minutes=180).process(cubes, cube[1]) + + +def test_invalid_method_for_period_exception(): + """Test that providing a period diagnostic and attempting to apply an + unsuitable interpolation method raises an exception.""" + + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] + data = np.ones((5, 5), dtype=np.float32) + cube = multi_time_cube(times, data, "latlon", bounds=True) + + msg = "Period diagnostics can only be temporally interpolated" + with pytest.raises(ValueError, match=msg): + TemporalInterpolation( + interval_in_minutes=180, interpolation_method="solar" + ).process(cube[0], cube[1]) + + +def test_period_without_chosen_type_exception(): + """Test that providing a period diagnostic but not specifying a type from + the available minimum, maximum, or accumulation raises an exception.""" + + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] + data = np.ones((5, 5), dtype=np.float32) + cube = multi_time_cube(times, data, "latlon", bounds=True) + + msg = "A type of period must be specified when interpolating" + with pytest.raises(ValueError, match=msg): + TemporalInterpolation(interval_in_minutes=180).process(cube[0], cube[1]) + + +@pytest.mark.parametrize( + "input_times,expected_time_bounds,expected_fp_bounds", + [ + ( + [3, 6, 9], # Forecast reference time ends up as 0 AM. + [[1509505200, 1509516000], [1509516000, 1509526800]], # 3-6, 6-9 AM + [[10800, 21600], [21600, 32400]], # T+3 - T+6, T+6 - T+9 + ), + ( + [3, 4, 6, 9], # Forecast reference time ends up as 2 AM. + [ + [1509505200, 1509508800], + [1509508800, 1509516000], + [1509516000, 1509526800], + ], # 3-4, 4-6, 6-9 AM + [ + [3600, 7200], + [7200, 14400], + [14400, 25200], + ], # T+1 - T+2, T+2 - T+4, T+4 - T+7 + ), + ( + [3, 4, 9], # Forecast reference time ends up as 2 AM. + [[1509505200, 1509508800], [1509508800, 1509526800]], # 3-4, 4-9 AM + [[3600, 7200], [7200, 25200]], # T+1 - T+2, T+2 - T+7 + ), + ], +) +def test_add_bounds(input_times, expected_time_bounds, expected_fp_bounds): + """Test the add bounds method creates the expected bounds for interpolated + data with different interpolated time intervals.""" + + times = [datetime.datetime(2017, 11, 1, hour) for hour in input_times] + + data = np.ones((5, 5), dtype=np.float32) + cube = multi_time_cube(times, data, "latlon") + # The first of the input times is used to represent the earlier of the + # input cubes. The other input time represent the interpolated times. + cube_t0 = cube[0] + interpolated_cube = cube[1:].copy() + + # Note the interval_in_minutes defined here is not used but required. + TemporalInterpolation(interval_in_minutes=60).add_bounds(cube_t0, interpolated_cube) + + assert (interpolated_cube.coord("time").bounds == expected_time_bounds).all() + assert ( + interpolated_cube.coord("forecast_period").bounds == expected_fp_bounds + ).all() + + +@pytest.mark.parametrize( + "kwargs,values,offsets,expected", + [ + # Equal adjacent accumulations, divided into equal shorter periods. + ( + {"interval_in_minutes": 180, "accumulation": True}, + [5, 5], + [3, 6], + [2.5, 2.5], + ), + # Equal adjacent e.g. period maxes, shorter periods have the same max. + ({"interval_in_minutes": 180, "maximum": True}, [5, 5], [3, 6], [5, 5],), + # Equal adjacent e.g. period minimums, shorter periods have the same + # min. + ({"interval_in_minutes": 180, "minimum": True}, [5, 5], [3, 6], [5, 5],), + # Trend of increasing accumulations with time, which is reflected + # in the shorter periods generated. + ( + {"interval_in_minutes": 180, "accumulation": True}, + [3, 9], + [3, 6], + [3.375, 5.625], + ), + # Trend of increasing maxes with time, which is reflected in the + # shorter periods generated. + ({"interval_in_minutes": 180, "maximum": True}, [3, 9], [3, 6], [6, 9],), + # Later input period minimum is 9, expect all new periods to be >= 9 + ({"interval_in_minutes": 180, "minimum": True}, [3, 9], [3, 6], [9, 9],), + # Trend of increasing accumulations with time, which is reflected + # in the shorter periods generated. + ( + {"interval_in_minutes": 120, "accumulation": True}, + [0, 9], + [2, 4, 6], + [1, 3, 5], + ), + # Trend of increasing maxes with time, which is reflected in the + # shorter periods generated. + ({"interval_in_minutes": 120, "maximum": True}, [0, 9], [2, 4, 6], [3, 6, 9],), + # Trend of increasing maxes with time, which is reflected in the + # shorter periods generated. + ({"interval_in_minutes": 120, "minimum": True}, [0, 9], [2, 4, 6], [9, 9, 9],), + # Later input period is 0, expect all new periods to be 0 + ( + {"interval_in_minutes": 120, "accumulation": True}, + [9, 0], + [2, 4, 6], + [0, 0, 0], + ), + # Later input period max is 0, expect all new periods to be 0 + ({"interval_in_minutes": 120, "maximum": True}, [9, 0], [2, 4, 6], [0, 0, 0],), + # Later input period minimum is 0, expect all new periods to be >= 0 + ({"interval_in_minutes": 120, "minimum": True}, [9, 0], [2, 4, 6], [6, 3, 0],), + # Equal adjacent accumulations, divided into unequal shorter periods. + ( + {"times": [datetime.datetime(2017, 11, 1, 4)], "accumulation": True}, + [6, 6], + [1, 6], + [1, 5], + ), + # Equal adjacent e.g. period maxes, unequal shorter periods have the + # same max. + ( + {"times": [datetime.datetime(2017, 11, 1, 4)], "maximum": True}, + [6, 6], + [1, 6], + [6, 6], + ), + # Equal adjacent e.g. period minimums, unequal shorter periods have + # the same min. + ( + {"times": [datetime.datetime(2017, 11, 1, 4)], "minimum": True}, + [6, 6], + [1, 6], + [6, 6], + ), + # Trend of increasing accumulations with time, which is reflected + # in the unequal shorter periods generated. + ( + {"times": [datetime.datetime(2017, 11, 1, 4)], "accumulation": True}, + [0, 9], + [1, 6], + [0.25, 8.75], + ), + # Trend of decreasing accumulations with time, which is reflected + # in the unequal shorter periods generated. + ( + {"times": [datetime.datetime(2017, 11, 1, 4)], "accumulation": True}, + [12, 3], + [1, 6], + [0.75, 2.25], + ), + ], +) +def test_process_periods(kwargs, values, offsets, expected): + """Test the process method when applied to period diagnostics, some + accumlations and some not.""" + + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] + npoints = 10 + data = np.stack( + [ + np.full((npoints, npoints), values[0], dtype=np.float32), + np.full((npoints, npoints), values[1], dtype=np.float32), + ] + ) + cube = multi_time_cube(times, data, "latlon", bounds=True) + + result = TemporalInterpolation(**kwargs).process(cube[0], cube[1]) + + for i, (offset, value) in enumerate(zip(offsets, expected)): + expected_data = np.full((npoints, npoints), value) + expected_time = 1509505200 + (offset * 3600) + expected_fp = (6 + offset) * 3600 + + assert result[i].coord("time").points[0] == expected_time + assert result[i].coord("forecast_period").points[0] == expected_fp + assert result[i].coord("time").points.dtype == "int64" + assert result[i].coord("forecast_period").points.dtype == "int32" + if value is not None: + np.testing.assert_almost_equal(result[i].data, expected_data) From c35a863518b65605b8ee19b925db279aef69b11d Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Fri, 16 Feb 2024 15:45:46 +0000 Subject: [PATCH 09/19] Enforce float type on renormalised accumulations. --- improver/utilities/temporal_interpolation.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index eae0df91f6..b17f60055d 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -40,6 +40,7 @@ from numpy import ndarray from improver import BasePlugin +from improver.metadata.constants import FLOAT_DTYPE from improver.metadata.constants.time_types import TIME_COORDS from improver.utilities.cube_manipulation import MergeCubes from improver.utilities.round import round_close @@ -511,6 +512,7 @@ def _calculate_accumulation( interpolated_total = np.sum(interpolated_cube.data, axis=time_coord) renormalisation = period_reference.data / interpolated_total interpolated_cube.data *= renormalisation + interpolated_cube.data = interpolated_cube.data.astype(FLOAT_DTYPE) def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: """ From c7f749a7d54ce283b76e0c0fcd44ce7ca5a27d93 Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Mon, 19 Feb 2024 12:32:20 +0000 Subject: [PATCH 10/19] Add acceptance tests for temporal interpolation of period data. --- improver_tests/acceptance/SHA256SUMS | 31 +++++++++++-- .../acceptance/test_temporal_interpolate.py | 46 +++++++++++++++++++ 2 files changed, 73 insertions(+), 4 deletions(-) diff --git a/improver_tests/acceptance/SHA256SUMS b/improver_tests/acceptance/SHA256SUMS index 52fc6fcee5..a5085d1fa0 100644 --- a/improver_tests/acceptance/SHA256SUMS +++ b/improver_tests/acceptance/SHA256SUMS @@ -632,6 +632,20 @@ aa54f6292c8e948eb02c3f127d8240252b1e626db330bee26d00fa29c0edb60e ./neighbour-fi 478ca4750aae458380c6abe730f9bdfd664e4007e5a298ae59d01a4356f99ece ./neighbour-finding/outputs/nearest_minimum_dz_uk_kgo.nc da22d11a9597890fa31f578e21dbb274b4c0fc5b7a8a9f2df5d7d510c2f83b70 ./neighbour-finding/outputs/nearest_uk_kgo.nc 936c9fe815274d2bd9d18f5cc79719e15ae7a4424ee8121694ea701815b7b82b ./neighbour-finding/outputs/nearest_uk_kgo_some_unset_wmo_ids_unique_ids.nc +4a1be5d407fcaf353e6960b64331cbafc3c67ff25ee51d2d19328804f7f6da17 ./new_temporal_interpolate/accumulation/20240217T1900Z-PT0033H00M-precipitation_accumulation-PT03H.nc +eeb4ed8ec35950edd9472f08c6feeacc10e26fb22fb541b996c5275a8be1f872 ./new_temporal_interpolate/accumulation/20240217T2200Z-PT0036H00M-precipitation_accumulation-PT03H.nc +eafaa596ac4d54774ee3e2890ec63bc3c8da7e9f6c32198fe930b085aff715ac ./new_temporal_interpolate/accumulation/kgo.nc +7111e3bf02103c91e38ad7b3fa752b6ffd93c76938f5bad1e7742ed4fcf12f08 ./new_temporal_interpolate/basic/20190116T0900Z-PT0033H00M-temperature_at_screen_level.nc +46c3d7d58af0dffbbdfd8a9e37094a446420ceafbf5d9aadb1e9667bf0317cb6 ./new_temporal_interpolate/basic/20190116T1200Z-PT0036H00M-temperature_at_screen_level.nc +d9a76a7fd4f8a614895072d0897eea8aaae6c23e53d8af91ca51d0fcef75664d ./new_temporal_interpolate/basic/kgo.nc +d2e4b9ca9dd1dcf0eca81c52e817e8297a1653abba77d0cbdcf09da330aa5206 ./new_temporal_interpolate/basic/kgo_t1.nc +ca0be5b3cc30779a21cc2f7131531dbd4da03aaef3b06418d8ba6749dc27e055 ./new_temporal_interpolate/period/20240217T0300Z-PT0016H00M-wind_gust_at_10m_max-PT03H.nc +7efb8363fc7094642351aa7a48a15fbf66795877c50afc00cb2f37f254828a5a ./new_temporal_interpolate/period/20240217T0600Z-PT0019H00M-wind_gust_at_10m_max-PT03H.nc +3e1b0074d7e4303800957dbfe14e28624b543a76b46c63077834c69e37a10bd3 ./new_temporal_interpolate/period/kgo.nc +eb6f7c3f646c4c51a0964b9a19367f43d6e3762ff5523b982cfaf7bf2610f091 ./new_temporal_interpolate/uv/20181220T0900Z-PT0021H00M-uv_index.nc +e3b8f51a0be52c4fead55f95c0e3da29ee3d93f92deed26314e60ad43e8fd5ef ./new_temporal_interpolate/uv/20181220T1200Z-PT0024H00M-uv_index.nc +b3fde693b3a8e144cb8f9ee9ff23c51ef92701858667cff850b2a49986bacaab ./new_temporal_interpolate/uv/kgo_t1.nc +1065ae1f25e6bc6df8d02e61c1f8ef92ab3dae679595d5165bd94d9c740adb2c ./new_temporal_interpolate/uv/kgo_t1_daynight.nc 83fa4bfc407a5ecd74f3d7de807e3d9d567e1052bbaa366e77cf75a80951494f ./normalise-to-reference/percentile/input_rain_rate.nc ef00f9b76df85efe1a77d4c81f735789337b1e882daf13bee9f28ea1fda9c810 ./normalise-to-reference/percentile/input_sleet_rate.nc 4974aa5b23458fab845c520ea3ca3b546062e5d045dbc51ec91ecef92490ab5b ./normalise-to-reference/percentile/input_snow_rate.nc @@ -758,6 +772,9 @@ f862eb29b6ec18916d25e11080feed5e1820ff1410d9cae4d0922b0d9aa0d058 ./resolve-wind ce3ecb08d6cc3569d84fc6c17502018582a5eab3f7d6411eb3240fd4563bf74f ./shower-condition-probability/cloud_input.nc e9086965416e3a4c02df5644eef55449acac46842399749b553f20c3b065ed08 ./shower-condition-probability/convection_input.nc 47d18297114668157bf509d991b195c5d7ca8b5ead35cc833514d2452bd644da ./shower-condition-probability/kgo.nc +d5916af8c067620094fe6383e42b0fafe8debe10f2089d6a5d9c7a698ab2fc1a ./shrunk/period/v2_20240217T0300Z-PT0016H00M-wind_gust_at_10m_max-PT03H.nc +2c8a64a8250e301e89a4a9cc804f8274ca042e4bfec7dc168a5c9b5e1c5ca375 ./shrunk/period/v2_20240217T0600Z-PT0019H00M-wind_gust_at_10m_max-PT03H.nc +967857a367f2a30657e4a7176566101f41862fbcb4069f35763eedba305b0f7e ./shrunk/period/v2_kgo.nc 77d35d36e52a65aa782e969e791f79c5558be8b0cf838fdc00bef3b4e75e6723 ./sleet_probability/basic/half_prob_snow_falling_level.nc e592f74bec62c53fe906d6f1c1d87f4acec5f2660457c0db98ffc00c64066663 ./sleet_probability/basic/kgo.nc efff8720d452591b6531a279bad6c0168db2336fa497c1f1087ecfcb9321af92 ./sleet_probability/basic/tenth_prob_snow_falling_level.nc @@ -820,14 +837,20 @@ c0b4643cf8fb0de1688c4f432ae5abb3c898582f5b41e996ebee03f679965314 ./temp-lapse-r 78c40fb411b63bd8c4cfef38f163b6f5df30679873211b2d7ba690ddeb9bdecb ./temp-lapse-rate/realizations/enukx_orography.nc 186f5d2cda6708e09459b69b92d04791a12c938f2aad9518c0c5d0a5b65b1e94 ./temp-lapse-rate/realizations/enukx_temperature.nc 8185d38dde574099b138b6cacf8b2aaee624cd0e19b98979616a8f675c317d84 ./temp-lapse-rate/realizations/kgo.nc +40612597709a889fa9b90f5ee6feb8c5bf893268db727db0084d4be1eee4a485 ./temporal-interpolate/accumulation/20240217T1900Z-PT0033H00M-precipitation_accumulation-PT03H.nc +173fa6e5548674af8ae3b48b4c305e8b522187d7b7386251c1af7cf600b369bd ./temporal-interpolate/accumulation/20240217T2200Z-PT0036H00M-precipitation_accumulation-PT03H.nc +e53201fb2d1226225546c5b7f06400ed4802aeb38484407c796e97480f614bfe ./temporal-interpolate/accumulation/kgo.nc 7111e3bf02103c91e38ad7b3fa752b6ffd93c76938f5bad1e7742ed4fcf12f08 ./temporal-interpolate/basic/20190116T0900Z-PT0033H00M-temperature_at_screen_level.nc 46c3d7d58af0dffbbdfd8a9e37094a446420ceafbf5d9aadb1e9667bf0317cb6 ./temporal-interpolate/basic/20190116T1200Z-PT0036H00M-temperature_at_screen_level.nc -42881311798b72302196546d9d106928833ccc28cd87867831d51265022f6227 ./temporal-interpolate/basic/kgo.nc -205f894969fdc5d8295eb44228e64fbacb74f668ec7f02373d1e5f56e064ce2b ./temporal-interpolate/basic/kgo_t1.nc +d9a76a7fd4f8a614895072d0897eea8aaae6c23e53d8af91ca51d0fcef75664d ./temporal-interpolate/basic/kgo.nc +d2e4b9ca9dd1dcf0eca81c52e817e8297a1653abba77d0cbdcf09da330aa5206 ./temporal-interpolate/basic/kgo_t1.nc +d5916af8c067620094fe6383e42b0fafe8debe10f2089d6a5d9c7a698ab2fc1a ./temporal-interpolate/period/20240217T0300Z-PT0016H00M-wind_gust_at_10m_max-PT03H.nc +2c8a64a8250e301e89a4a9cc804f8274ca042e4bfec7dc168a5c9b5e1c5ca375 ./temporal-interpolate/period/20240217T0600Z-PT0019H00M-wind_gust_at_10m_max-PT03H.nc +db7ea71aef580b1d49406e442bd2a3783255a9ccb990e3827fbc0d422bb7f4c8 ./temporal-interpolate/period/kgo.nc eb6f7c3f646c4c51a0964b9a19367f43d6e3762ff5523b982cfaf7bf2610f091 ./temporal-interpolate/uv/20181220T0900Z-PT0021H00M-uv_index.nc e3b8f51a0be52c4fead55f95c0e3da29ee3d93f92deed26314e60ad43e8fd5ef ./temporal-interpolate/uv/20181220T1200Z-PT0024H00M-uv_index.nc -71e137948f4de612e8c5b48134a42f171b2be5d9dd81ac9d56688e943e7f5c3f ./temporal-interpolate/uv/kgo_t1.nc -68a8f0f2c4f62ef8c2d3a597bf3c0de69e769d4920ad7f4129dbcf9e0f95ef2f ./temporal-interpolate/uv/kgo_t1_daynight.nc +b3fde693b3a8e144cb8f9ee9ff23c51ef92701858667cff850b2a49986bacaab ./temporal-interpolate/uv/kgo_t1.nc +1065ae1f25e6bc6df8d02e61c1f8ef92ab3dae679595d5165bd94d9c740adb2c ./temporal-interpolate/uv/kgo_t1_daynight.nc ac93ed67c9947547e5879af6faaa329fede18afd822c720ac3afcb18fa41077a ./threshold/basic/input.nc eb3fdc9400401ec47d95961553aed452abcbd91891d0fbca106b3a05131adaa9 ./threshold/basic/kgo.nc 6b50fa16b663869b3e3fbff36197603886ff7383b2df2a8ba92579bcc9461a16 ./threshold/below_threshold/kgo.nc diff --git a/improver_tests/acceptance/test_temporal_interpolate.py b/improver_tests/acceptance/test_temporal_interpolate.py index ff27839c9b..124e46b731 100644 --- a/improver_tests/acceptance/test_temporal_interpolate.py +++ b/improver_tests/acceptance/test_temporal_interpolate.py @@ -112,3 +112,49 @@ def test_daynight(tmp_path): ] run_cli(args) acc.compare(output_path, kgo_path) + + +def test_period_max(tmp_path): + """Test interpolation of an period maximum.""" + kgo_dir = acc.kgo_root() / "temporal-interpolate/period" + kgo_path = kgo_dir / "kgo.nc" + input_paths = [ + kgo_dir / f"20240217T{v:04}Z-PT{l:04}H00M-wind_gust_at_10m_max-PT03H.nc" + for v, l in ((300, 16), (600, 19)) + ] + output_path = tmp_path / "output.nc" + args = [ + *input_paths, + "--times", + "20240217T0430Z", + "--interpolation-method", + "linear", + "--maximum", + "--output", + output_path, + ] + run_cli(args) + acc.compare(output_path, kgo_path) + + +def test_accumulation(tmp_path): + """Test interpolation of an accumulation.""" + kgo_dir = acc.kgo_root() / "temporal-interpolate/accumulation" + kgo_path = kgo_dir / "kgo.nc" + input_paths = [ + kgo_dir / f"20240217T{v:04}Z-PT{l:04}H00M-precipitation_accumulation-PT03H.nc" + for v, l in ((1900, 33), (2200, 36)) + ] + output_path = tmp_path / "output.nc" + args = [ + *input_paths, + "--times", + "20240217T2000Z,20240217T2100Z", + "--interpolation-method", + "linear", + "--accumulation", + "--output", + output_path, + ] + run_cli(args) + acc.compare(output_path, kgo_path) From b12be5c9bc540fd5d3b5aaa8ecc2acd69099198a Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Mon, 19 Feb 2024 13:45:50 +0000 Subject: [PATCH 11/19] Add an additional test to check for overlapping period diagnostics. --- improver/utilities/temporal_interpolation.py | 16 ++++++++++++++++ .../utilities/test_TemporalInterpolation.py | 19 +++++++++++++++++++ 2 files changed, 35 insertions(+) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index b17f60055d..6ee373f45b 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -536,6 +536,7 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: not found in the period_interpolation_methods list. ValueError: A period diagnostic is being processed without the type of period being specified. + ValueError: Period diagnostics with overlapping periods. CoordinateNotFoundError: The input cubes contain no time coordinate. ValueError: Cubes contain multiple validity times. @@ -565,6 +566,21 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: "period diagnostic. This may be a period maximum, minimum" " or an accumulation." ) + cube_interval = ( + cube_t1.coord("time").points[0] - cube_t0.coord("time").points[0] + ) + (period,) = np.diff(cube_t1.coord("time").bounds[0]) + if not cube_interval == period: + raise ValueError( + "The diagnostic provided represents the period " + f"{period / 3600} hours. The interval between the " + f"diagnostics is {cube_interval / 3600} hours. Temporal " + "interpolation can only be applied to a period " + "diagnostic provided at intervals that match the " + "diagnostic period such that all points in time are " + "captured by only one of the inputs and do not overlap." + ) + self.period_inputs = True try: diff --git a/improver_tests/utilities/test_TemporalInterpolation.py b/improver_tests/utilities/test_TemporalInterpolation.py index cb0eed9c29..36cb8889bb 100644 --- a/improver_tests/utilities/test_TemporalInterpolation.py +++ b/improver_tests/utilities/test_TemporalInterpolation.py @@ -667,6 +667,25 @@ def test_period_without_chosen_type_exception(): TemporalInterpolation(interval_in_minutes=180).process(cube[0], cube[1]) +def test_period_unequal_to_interval_exception(): + """Test that providing a period diagnostic where the represented + periods overlap raises an exception.""" + + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] + data = np.ones((5, 5), dtype=np.float32) + cube = multi_time_cube(times, data, "latlon", bounds=True) + for crd in ["time", "forecast_period"]: + bounds = cube.coord(crd).bounds + bounds = [[lower, upper + 3600] for lower, upper in bounds] + cube.coord(crd).bounds = bounds + + msg = "The diagnostic provided represents the period" + with pytest.raises(ValueError, match=msg): + TemporalInterpolation(interval_in_minutes=180, accumulation=True).process( + cube[0], cube[1] + ) + + @pytest.mark.parametrize( "input_times,expected_time_bounds,expected_fp_bounds", [ From 2bddb65dfc5a9c5e18aeca8742402814b410fab2 Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Mon, 19 Feb 2024 15:10:53 +0000 Subject: [PATCH 12/19] Add pass through if only target is the trailing input cube. --- improver/utilities/temporal_interpolation.py | 16 ++++++--- .../utilities/test_TemporalInterpolation.py | 34 +++++++++++++++++++ 2 files changed, 46 insertions(+), 4 deletions(-) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index 6ee373f45b..157c58f57f 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -342,7 +342,7 @@ def solar_interpolate(self, diag_cube: Cube, interpolated_cube: Cube) -> CubeLis A list of cubes interpolated to the desired times. """ - interpolated_cubes = iris.cube.CubeList() + interpolated_cubes = CubeList() (lats, lons) = self.calc_lats_lons(diag_cube) prev_data = diag_cube[0].data next_data = diag_cube[1].data @@ -430,7 +430,7 @@ def daynight_interpolate(interpolated_cube: Cube) -> CubeList: else: interpolated_cube.data[index] = 0.0 - return iris.cube.CubeList(list(interpolated_cube.slices_over("time"))) + return CubeList(list(interpolated_cube.slices_over("time"))) @staticmethod def add_bounds(cube_t0: Cube, interpolated_cube: Cube): @@ -606,6 +606,14 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: time_list = self.construct_time_list(initial_time, final_time) + # If the target output time is the same as time at which the + # trailing input is valid, just return it unchanged. + if ( + len(time_list[0][1]) == 1 + and time_list[0][1][0] == cube_t1.coord("time").cell(0).point + ): + return CubeList([cube_t1]) + # If the units of the two cubes are degrees, assume we are dealing with # directions. Convert the directions to complex numbers so # interpolations (esp. the 0/360 wraparound) are handled in a sane @@ -621,7 +629,7 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: period_reference = cube_t1.copy() cube_t1.data /= np.diff(cube_t1.coord("forecast_period").bounds[0])[0] - cubes = iris.cube.CubeList([cube_t0, cube_t1]) + cubes = CubeList([cube_t0, cube_t1]) cube = MergeCubes()(cubes) interpolated_cube = cube.interpolate(time_list, iris.analysis.Linear()) @@ -650,7 +658,7 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: ) self.enforce_time_coords_dtype(interpolated_cube) - interpolated_cubes = iris.cube.CubeList() + interpolated_cubes = CubeList() if self.interpolation_method == "solar": interpolated_cubes = self.solar_interpolate(cube, interpolated_cube) elif self.interpolation_method == "daynight": diff --git a/improver_tests/utilities/test_TemporalInterpolation.py b/improver_tests/utilities/test_TemporalInterpolation.py index 36cb8889bb..50047d7bdd 100644 --- a/improver_tests/utilities/test_TemporalInterpolation.py +++ b/improver_tests/utilities/test_TemporalInterpolation.py @@ -857,3 +857,37 @@ def test_process_periods(kwargs, values, offsets, expected): assert result[i].coord("forecast_period").points.dtype == "int32" if value is not None: np.testing.assert_almost_equal(result[i].data, expected_data) + + +@pytest.mark.parametrize( + "kwargs", + ( + [ + {"interval_in_minutes": 360, "accumulation": True}, + {"times": [datetime.datetime(2017, 11, 1, 9)], "accumulation": True}, + ] + ), +) +def test_process_return_input(kwargs): + """Test the process method returns an unmodified cube when the + target time is identical to that of the trailing input.""" + + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] + npoints = 5 + data = np.stack( + [ + np.full((npoints, npoints), 3, dtype=np.float32), + np.full((npoints, npoints), 4, dtype=np.float32), + ] + ) + cube = multi_time_cube(times, data, "latlon", bounds=True) + + # Slice here to keep memory addresses consistent when passed in to + # plugin. + cube_0 = cube[0] + cube_1 = cube[1] + result = TemporalInterpolation(**kwargs).process(cube_0, cube_1) + + # assert that the object returned is the same one in memory that was + # passed in. + assert result[0] is cube_1 From c762517412a59efcd0cfcab01f12d6794d285fcd Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Mon, 19 Feb 2024 17:10:14 +0000 Subject: [PATCH 13/19] Fix up multi-realization accumulation handling and add test coverage. --- improver/utilities/temporal_interpolation.py | 5 +++-- .../utilities/test_TemporalInterpolation.py | 16 +++++++++++----- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index 157c58f57f..078bceb701 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -500,10 +500,11 @@ def _calculate_accumulation( ) interpolated_cube.data = period_rates - # Multiply the averate rate by the length of each period to get a new + # Multiply the average rate by the length of each period to get a new # accumulation. new_periods = np.diff(interpolated_cube.coord("forecast_period").bounds) - new_periods = new_periods[:, np.newaxis] + for _ in range(interpolated_cube.ndim - new_periods.ndim): + new_periods = np.expand_dims(new_periods, axis=1) interpolated_cube.data = np.multiply(new_periods, interpolated_cube.data) # Renormalise the total of the new periods to ensure it matches the diff --git a/improver_tests/utilities/test_TemporalInterpolation.py b/improver_tests/utilities/test_TemporalInterpolation.py index 50047d7bdd..22103d9651 100644 --- a/improver_tests/utilities/test_TemporalInterpolation.py +++ b/improver_tests/utilities/test_TemporalInterpolation.py @@ -736,6 +736,7 @@ def test_add_bounds(input_times, expected_time_bounds, expected_fp_bounds): ).all() +@pytest.mark.parametrize("realizations", (None, [0, 1, 2])) @pytest.mark.parametrize( "kwargs,values,offsets,expected", [ @@ -830,24 +831,29 @@ def test_add_bounds(input_times, expected_time_bounds, expected_fp_bounds): ), ], ) -def test_process_periods(kwargs, values, offsets, expected): +def test_process_periods(kwargs, values, offsets, expected, realizations): """Test the process method when applied to period diagnostics, some - accumlations and some not.""" + accumlations and some not. Test with and without multiple realizations.""" times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] - npoints = 10 + npoints = 5 data = np.stack( [ np.full((npoints, npoints), values[0], dtype=np.float32), np.full((npoints, npoints), values[1], dtype=np.float32), ] ) - cube = multi_time_cube(times, data, "latlon", bounds=True) + cube = multi_time_cube( + times, data, "latlon", bounds=True, realizations=realizations + ) result = TemporalInterpolation(**kwargs).process(cube[0], cube[1]) for i, (offset, value) in enumerate(zip(offsets, expected)): - expected_data = np.full((npoints, npoints), value) + if realizations: + expected_data = np.full((len(realizations), npoints, npoints), value) + else: + expected_data = np.full((npoints, npoints), value) expected_time = 1509505200 + (offset * 3600) expected_fp = (6 + offset) * 3600 From cf1fc70f3bae6cffbae52d9539cdcd5d010a2bfe Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Tue, 20 Feb 2024 08:26:21 +0000 Subject: [PATCH 14/19] Improve comment. --- improver/utilities/temporal_interpolation.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index 078bceb701..34d8eaa204 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -624,7 +624,10 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: cube_t1.data = WindDirection.deg_to_complex(cube_t1.data) # Convert accumulations into rates to allow interpolation using trends - # in the data and to accommodate non-uniform output intervals. + # in the data and to accommodate non-uniform output intervals. This also + # accommodates cube_t0 and cube_t1 representing different periods of + # accumulation, for example where the forecast period interval changes + # in an NWP model's output. if self.accumulation: cube_t0.data /= np.diff(cube_t0.coord("forecast_period").bounds[0])[0] period_reference = cube_t1.copy() From 31eca75dfe3ce173fc0d23d7592b3ceebf6cc5e8 Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Tue, 20 Feb 2024 09:22:17 +0000 Subject: [PATCH 15/19] Add unit tests to cover unequal length inputs when interpolating accumulation diagnostics. --- .../utilities/test_TemporalInterpolation.py | 88 +++++++++++++++++++ 1 file changed, 88 insertions(+) diff --git a/improver_tests/utilities/test_TemporalInterpolation.py b/improver_tests/utilities/test_TemporalInterpolation.py index 22103d9651..9f6d945a9e 100644 --- a/improver_tests/utilities/test_TemporalInterpolation.py +++ b/improver_tests/utilities/test_TemporalInterpolation.py @@ -897,3 +897,91 @@ def test_process_return_input(kwargs): # assert that the object returned is the same one in memory that was # passed in. assert result[0] is cube_1 + + +@pytest.mark.parametrize("realizations", (None, [0, 1, 2])) +@pytest.mark.parametrize( + "kwargs,values,offsets,expected", + [ + # Unequal input periods and accumulations give effective rates of + # 1 mm/hr and 2 mm/hr at the start and end of the period. This gives + # a gradient of 1/6 mm/hr across the period which results in the + # expected 3-hour accumulations returned across the period. + ({"interval_in_minutes": 180, "accumulation": True}, [3, 12], [3, 6], [5, 7],), + # Unequal input periods and accumulations give effective rates of + # 2 mm/hr and 1 mm/hr at the start and end of the period. This gives + # a gradient of -1/6 mm/hr across the period which results in the + # expected 3-hour accumulations returned across the period. + ( + {"interval_in_minutes": 180, "accumulation": True}, + [6, 6], + [3, 6], + [3.5, 2.5], + ), + # Unequal input periods and accumulations give a consistent effective + # rate of 1 mm/hr across the the period. This results in equal + # accumulations across the two returned 3-hour periods. + ({"interval_in_minutes": 180, "accumulation": True}, [3, 6], [3, 6], [3, 3],), + # Unequal input periods and accumulations give a consistent effective + # rate of 1 mm/hr across the the period. The unequal output periods + # split the total accumulation as expected. + ( + {"times": [datetime.datetime(2017, 11, 1, 4)], "accumulation": True}, + [3, 6], + [1, 6], + [1, 5], + ), + # Unequal input periods and accumulations give effective rates of + # 1 mm/hr and 1.5 mm/hr at the start and end of the period. This gives + # a gradient of 1/12 mm/hr across the period. The unequal output periods + # divide up the total accumulation in line with this as expected. + ( + {"times": [datetime.datetime(2017, 11, 1, 5)], "accumulation": True}, + [3, 9], + [2, 6], + [2.6, 6.4], + ), + ], +) +def test_process_accumulation_unequal_inputs( + kwargs, values, offsets, expected, realizations +): + """Test that the expected values are returned when the accumulation inputs + are of different periods. The accumulations are converted to rates using + each input cube's period prior to interpolation which allos for this.""" + + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] + npoints = 5 + data = np.stack( + [ + np.full((npoints, npoints), values[0], dtype=np.float32), + np.full((npoints, npoints), values[1], dtype=np.float32), + ] + ) + cube = multi_time_cube( + times, data, "latlon", bounds=True, realizations=realizations + ) + cube_0 = cube[0] + cube_1 = cube[1] + + for crd in ["time", "forecast_period"]: + bounds = cube_0.coord(crd).bounds + bounds = [[lower + 10800, upper] for lower, upper in bounds] + cube_0.coord(crd).bounds = bounds + + result = TemporalInterpolation(**kwargs).process(cube_0, cube_1) + + for i, (offset, value) in enumerate(zip(offsets, expected)): + if realizations: + expected_data = np.full((len(realizations), npoints, npoints), value) + else: + expected_data = np.full((npoints, npoints), value) + expected_time = 1509505200 + (offset * 3600) + expected_fp = (6 + offset) * 3600 + + assert result[i].coord("time").points[0] == expected_time + assert result[i].coord("forecast_period").points[0] == expected_fp + assert result[i].coord("time").points.dtype == "int64" + assert result[i].coord("forecast_period").points.dtype == "int32" + if value is not None: + np.testing.assert_almost_equal(result[i].data, expected_data) From f4836386b9acf1f325b8ef671ce58d35c5235648 Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Wed, 21 Feb 2024 11:45:37 +0000 Subject: [PATCH 16/19] Modify arguments from maximum and minimum to max and min to align with other CLIs like combine. --- improver/cli/temporal_interpolate.py | 12 ++++----- improver/utilities/temporal_interpolation.py | 26 +++++++++---------- improver_tests/acceptance/SHA256SUMS | 17 ------------ .../acceptance/test_temporal_interpolate.py | 2 +- .../utilities/test_TemporalInterpolation.py | 24 ++++++++--------- 5 files changed, 32 insertions(+), 49 deletions(-) diff --git a/improver/cli/temporal_interpolate.py b/improver/cli/temporal_interpolate.py index 3dc27087c4..168b68dedc 100755 --- a/improver/cli/temporal_interpolate.py +++ b/improver/cli/temporal_interpolate.py @@ -44,8 +44,8 @@ def process( times: cli.comma_separated_list = None, interpolation_method="linear", accumulation: bool = False, - maximum: bool = False, - minimum: bool = False, + max: bool = False, + min: bool = False, ): """Interpolate data between validity times. @@ -85,12 +85,12 @@ def process( intervals matches the total across the period from the coarser intervals. Trends between adjacent input periods will be used to provide variation across the interpolated periods. - maximum: + max: Set True if the diagnostic being temporally interpolated is a period maximum. Trends between adjacent input periods will be used to provide variation across the interpolated periods where these are consistent with the inputs. - minimum: + min: Set True if the diagnostic being temporally interpolated is a period minimum. Trends between adjacent input periods will be used to provide variation across the interpolated periods where these @@ -119,7 +119,7 @@ def process( times=times, interpolation_method=interpolation_method, accumulation=accumulation, - maximum=maximum, - minimum=minimum, + max=max, + min=min, )(start_cube, end_cube) return MergeCubes()(result) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index 34d8eaa204..3225f86758 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -99,8 +99,8 @@ def __init__( times: Optional[List[datetime]] = None, interpolation_method: str = "linear", accumulation: bool = False, - maximum: bool = False, - minimum: bool = False, + max: bool = False, + min: bool = False, ) -> None: """ Initialise class. @@ -126,12 +126,12 @@ def __init__( that the total across the period constructed from the shorter intervals matches the total across the period from the coarser intervals. - maximum: + max: Set True if the diagnostic being temporally interpolated is a period maximum. Trends between adjacent input periods will be used to provide variation across the interpolated periods where these are consistent with the inputs. - minimum: + min: Set True if the diagnostic being temporally interpolated is a period minimum. Trends between adjacent input periods will be used to provide variation across the interpolated periods where these @@ -166,15 +166,15 @@ def __init__( ) self.interpolation_method = interpolation_method self.period_inputs = False - if np.sum([accumulation, maximum, minimum]) > 1: + if np.sum([accumulation, max, min]) > 1: raise ValueError( "Only one type of period diagnostics may be specified: " - f"accumulation = {accumulation}, maximum = {maximum}, " - f"minimum = {minimum}" + f"accumulation = {accumulation}, max = {max}, " + f"min = {min}" ) self.accumulation = accumulation - self.maximum = maximum - self.minimum = minimum + self.max = max + self.min = min def construct_time_list( self, initial_time: datetime, final_time: datetime @@ -561,10 +561,10 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: f"using these methods: {self.period_interpolation_methods}.\n" f"Currently selected method is: {self.interpolation_method}." ) - if not any([self.accumulation, self.maximum, self.minimum]): + if not any([self.accumulation, self.max, self.min]): raise ValueError( "A type of period must be specified when interpolating a " - "period diagnostic. This may be a period maximum, minimum" + "period diagnostic. This may be a period max, min" " or an accumulation." ) cube_interval = ( @@ -652,11 +652,11 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: self._calculate_accumulation( cube_t0, period_reference, interpolated_cube ) - elif self.maximum: + elif self.max: interpolated_cube.data = np.minimum( cube_t1.data, interpolated_cube.data ) - elif self.minimum: + elif self.min: interpolated_cube.data = np.maximum( cube_t1.data, interpolated_cube.data ) diff --git a/improver_tests/acceptance/SHA256SUMS b/improver_tests/acceptance/SHA256SUMS index a5085d1fa0..82028532d4 100644 --- a/improver_tests/acceptance/SHA256SUMS +++ b/improver_tests/acceptance/SHA256SUMS @@ -632,20 +632,6 @@ aa54f6292c8e948eb02c3f127d8240252b1e626db330bee26d00fa29c0edb60e ./neighbour-fi 478ca4750aae458380c6abe730f9bdfd664e4007e5a298ae59d01a4356f99ece ./neighbour-finding/outputs/nearest_minimum_dz_uk_kgo.nc da22d11a9597890fa31f578e21dbb274b4c0fc5b7a8a9f2df5d7d510c2f83b70 ./neighbour-finding/outputs/nearest_uk_kgo.nc 936c9fe815274d2bd9d18f5cc79719e15ae7a4424ee8121694ea701815b7b82b ./neighbour-finding/outputs/nearest_uk_kgo_some_unset_wmo_ids_unique_ids.nc -4a1be5d407fcaf353e6960b64331cbafc3c67ff25ee51d2d19328804f7f6da17 ./new_temporal_interpolate/accumulation/20240217T1900Z-PT0033H00M-precipitation_accumulation-PT03H.nc -eeb4ed8ec35950edd9472f08c6feeacc10e26fb22fb541b996c5275a8be1f872 ./new_temporal_interpolate/accumulation/20240217T2200Z-PT0036H00M-precipitation_accumulation-PT03H.nc -eafaa596ac4d54774ee3e2890ec63bc3c8da7e9f6c32198fe930b085aff715ac ./new_temporal_interpolate/accumulation/kgo.nc -7111e3bf02103c91e38ad7b3fa752b6ffd93c76938f5bad1e7742ed4fcf12f08 ./new_temporal_interpolate/basic/20190116T0900Z-PT0033H00M-temperature_at_screen_level.nc -46c3d7d58af0dffbbdfd8a9e37094a446420ceafbf5d9aadb1e9667bf0317cb6 ./new_temporal_interpolate/basic/20190116T1200Z-PT0036H00M-temperature_at_screen_level.nc -d9a76a7fd4f8a614895072d0897eea8aaae6c23e53d8af91ca51d0fcef75664d ./new_temporal_interpolate/basic/kgo.nc -d2e4b9ca9dd1dcf0eca81c52e817e8297a1653abba77d0cbdcf09da330aa5206 ./new_temporal_interpolate/basic/kgo_t1.nc -ca0be5b3cc30779a21cc2f7131531dbd4da03aaef3b06418d8ba6749dc27e055 ./new_temporal_interpolate/period/20240217T0300Z-PT0016H00M-wind_gust_at_10m_max-PT03H.nc -7efb8363fc7094642351aa7a48a15fbf66795877c50afc00cb2f37f254828a5a ./new_temporal_interpolate/period/20240217T0600Z-PT0019H00M-wind_gust_at_10m_max-PT03H.nc -3e1b0074d7e4303800957dbfe14e28624b543a76b46c63077834c69e37a10bd3 ./new_temporal_interpolate/period/kgo.nc -eb6f7c3f646c4c51a0964b9a19367f43d6e3762ff5523b982cfaf7bf2610f091 ./new_temporal_interpolate/uv/20181220T0900Z-PT0021H00M-uv_index.nc -e3b8f51a0be52c4fead55f95c0e3da29ee3d93f92deed26314e60ad43e8fd5ef ./new_temporal_interpolate/uv/20181220T1200Z-PT0024H00M-uv_index.nc -b3fde693b3a8e144cb8f9ee9ff23c51ef92701858667cff850b2a49986bacaab ./new_temporal_interpolate/uv/kgo_t1.nc -1065ae1f25e6bc6df8d02e61c1f8ef92ab3dae679595d5165bd94d9c740adb2c ./new_temporal_interpolate/uv/kgo_t1_daynight.nc 83fa4bfc407a5ecd74f3d7de807e3d9d567e1052bbaa366e77cf75a80951494f ./normalise-to-reference/percentile/input_rain_rate.nc ef00f9b76df85efe1a77d4c81f735789337b1e882daf13bee9f28ea1fda9c810 ./normalise-to-reference/percentile/input_sleet_rate.nc 4974aa5b23458fab845c520ea3ca3b546062e5d045dbc51ec91ecef92490ab5b ./normalise-to-reference/percentile/input_snow_rate.nc @@ -772,9 +758,6 @@ f862eb29b6ec18916d25e11080feed5e1820ff1410d9cae4d0922b0d9aa0d058 ./resolve-wind ce3ecb08d6cc3569d84fc6c17502018582a5eab3f7d6411eb3240fd4563bf74f ./shower-condition-probability/cloud_input.nc e9086965416e3a4c02df5644eef55449acac46842399749b553f20c3b065ed08 ./shower-condition-probability/convection_input.nc 47d18297114668157bf509d991b195c5d7ca8b5ead35cc833514d2452bd644da ./shower-condition-probability/kgo.nc -d5916af8c067620094fe6383e42b0fafe8debe10f2089d6a5d9c7a698ab2fc1a ./shrunk/period/v2_20240217T0300Z-PT0016H00M-wind_gust_at_10m_max-PT03H.nc -2c8a64a8250e301e89a4a9cc804f8274ca042e4bfec7dc168a5c9b5e1c5ca375 ./shrunk/period/v2_20240217T0600Z-PT0019H00M-wind_gust_at_10m_max-PT03H.nc -967857a367f2a30657e4a7176566101f41862fbcb4069f35763eedba305b0f7e ./shrunk/period/v2_kgo.nc 77d35d36e52a65aa782e969e791f79c5558be8b0cf838fdc00bef3b4e75e6723 ./sleet_probability/basic/half_prob_snow_falling_level.nc e592f74bec62c53fe906d6f1c1d87f4acec5f2660457c0db98ffc00c64066663 ./sleet_probability/basic/kgo.nc efff8720d452591b6531a279bad6c0168db2336fa497c1f1087ecfcb9321af92 ./sleet_probability/basic/tenth_prob_snow_falling_level.nc diff --git a/improver_tests/acceptance/test_temporal_interpolate.py b/improver_tests/acceptance/test_temporal_interpolate.py index 124e46b731..b8c1c5c172 100644 --- a/improver_tests/acceptance/test_temporal_interpolate.py +++ b/improver_tests/acceptance/test_temporal_interpolate.py @@ -129,7 +129,7 @@ def test_period_max(tmp_path): "20240217T0430Z", "--interpolation-method", "linear", - "--maximum", + "--max", "--output", output_path, ] diff --git a/improver_tests/utilities/test_TemporalInterpolation.py b/improver_tests/utilities/test_TemporalInterpolation.py index 9f6d945a9e..f61ca9d6b3 100644 --- a/improver_tests/utilities/test_TemporalInterpolation.py +++ b/improver_tests/utilities/test_TemporalInterpolation.py @@ -257,7 +257,7 @@ def daynight_mask(): "TemporalInterpolation: Unknown interpolation method", ), # Invalid interpolation method requested ( - {"interval_in_minutes": 60, "maximum": True, "minimum": True}, + {"interval_in_minutes": 60, "max": True, "min": True}, "Only one type of period diagnostics may be specified:", ), # Invalid interpolation method requested ], @@ -656,7 +656,7 @@ def test_invalid_method_for_period_exception(): def test_period_without_chosen_type_exception(): """Test that providing a period diagnostic but not specifying a type from - the available minimum, maximum, or accumulation raises an exception.""" + the available min, max, or accumulation raises an exception.""" times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] data = np.ones((5, 5), dtype=np.float32) @@ -748,10 +748,10 @@ def test_add_bounds(input_times, expected_time_bounds, expected_fp_bounds): [2.5, 2.5], ), # Equal adjacent e.g. period maxes, shorter periods have the same max. - ({"interval_in_minutes": 180, "maximum": True}, [5, 5], [3, 6], [5, 5],), + ({"interval_in_minutes": 180, "max": True}, [5, 5], [3, 6], [5, 5],), # Equal adjacent e.g. period minimums, shorter periods have the same # min. - ({"interval_in_minutes": 180, "minimum": True}, [5, 5], [3, 6], [5, 5],), + ({"interval_in_minutes": 180, "min": True}, [5, 5], [3, 6], [5, 5],), # Trend of increasing accumulations with time, which is reflected # in the shorter periods generated. ( @@ -762,9 +762,9 @@ def test_add_bounds(input_times, expected_time_bounds, expected_fp_bounds): ), # Trend of increasing maxes with time, which is reflected in the # shorter periods generated. - ({"interval_in_minutes": 180, "maximum": True}, [3, 9], [3, 6], [6, 9],), + ({"interval_in_minutes": 180, "max": True}, [3, 9], [3, 6], [6, 9],), # Later input period minimum is 9, expect all new periods to be >= 9 - ({"interval_in_minutes": 180, "minimum": True}, [3, 9], [3, 6], [9, 9],), + ({"interval_in_minutes": 180, "min": True}, [3, 9], [3, 6], [9, 9],), # Trend of increasing accumulations with time, which is reflected # in the shorter periods generated. ( @@ -775,10 +775,10 @@ def test_add_bounds(input_times, expected_time_bounds, expected_fp_bounds): ), # Trend of increasing maxes with time, which is reflected in the # shorter periods generated. - ({"interval_in_minutes": 120, "maximum": True}, [0, 9], [2, 4, 6], [3, 6, 9],), + ({"interval_in_minutes": 120, "max": True}, [0, 9], [2, 4, 6], [3, 6, 9],), # Trend of increasing maxes with time, which is reflected in the # shorter periods generated. - ({"interval_in_minutes": 120, "minimum": True}, [0, 9], [2, 4, 6], [9, 9, 9],), + ({"interval_in_minutes": 120, "min": True}, [0, 9], [2, 4, 6], [9, 9, 9],), # Later input period is 0, expect all new periods to be 0 ( {"interval_in_minutes": 120, "accumulation": True}, @@ -787,9 +787,9 @@ def test_add_bounds(input_times, expected_time_bounds, expected_fp_bounds): [0, 0, 0], ), # Later input period max is 0, expect all new periods to be 0 - ({"interval_in_minutes": 120, "maximum": True}, [9, 0], [2, 4, 6], [0, 0, 0],), + ({"interval_in_minutes": 120, "max": True}, [9, 0], [2, 4, 6], [0, 0, 0],), # Later input period minimum is 0, expect all new periods to be >= 0 - ({"interval_in_minutes": 120, "minimum": True}, [9, 0], [2, 4, 6], [6, 3, 0],), + ({"interval_in_minutes": 120, "min": True}, [9, 0], [2, 4, 6], [6, 3, 0],), # Equal adjacent accumulations, divided into unequal shorter periods. ( {"times": [datetime.datetime(2017, 11, 1, 4)], "accumulation": True}, @@ -800,7 +800,7 @@ def test_add_bounds(input_times, expected_time_bounds, expected_fp_bounds): # Equal adjacent e.g. period maxes, unequal shorter periods have the # same max. ( - {"times": [datetime.datetime(2017, 11, 1, 4)], "maximum": True}, + {"times": [datetime.datetime(2017, 11, 1, 4)], "max": True}, [6, 6], [1, 6], [6, 6], @@ -808,7 +808,7 @@ def test_add_bounds(input_times, expected_time_bounds, expected_fp_bounds): # Equal adjacent e.g. period minimums, unequal shorter periods have # the same min. ( - {"times": [datetime.datetime(2017, 11, 1, 4)], "minimum": True}, + {"times": [datetime.datetime(2017, 11, 1, 4)], "min": True}, [6, 6], [1, 6], [6, 6], From e78b07da6a7b45c9ca0b3e98eed12ad21afbd95b Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Tue, 27 Feb 2024 15:37:49 +0000 Subject: [PATCH 17/19] Review changes. --- improver/utilities/temporal_interpolation.py | 17 +++++++++++------ .../utilities/test_TemporalInterpolation.py | 14 +++++++------- 2 files changed, 18 insertions(+), 13 deletions(-) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index 3225f86758..362e4a21ec 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -461,10 +461,9 @@ def add_bounds(cube_t0: Cube, interpolated_cube: Cube): raise CoordinateNotFoundError( f"Period diagnostic cube is missing expected coordinate: {crd}" ) - bounds = np.diff(interpolated_times) all_bounds = [] - for bound, time in zip(bounds, interpolated_times[1:]): - all_bounds.append([time - bound, time]) + for start, end in zip(interpolated_times[:-1], interpolated_times[1:]): + all_bounds.append([start, end]) interpolated_cube.coord(crd).bounds = all_bounds @staticmethod @@ -482,7 +481,7 @@ def _calculate_accumulation( Args: cube_t0: The input cube corresponding to the earlier time. - period_refrence: + period_reference: The input cube corresponding to the later time, with the values prior to conversion to rates. interpolated_cube: @@ -607,7 +606,7 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: time_list = self.construct_time_list(initial_time, final_time) - # If the target output time is the same as time at which the + # If the target output time is the same as the time at which the # trailing input is valid, just return it unchanged. if ( len(time_list[0][1]) == 1 @@ -647,7 +646,13 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: # if the inputs were period diagnostics. self.add_bounds(cube_t0, interpolated_cube) - # Apply suitable + # Apply suitable constraints to the returned values. + # - accumulations are renormalised to ensure the period total is + # unchanged when broken into shorter periods. + # - period maximums are enforced to not exceed the original + # maximum that occurred across the whole longer period. + # - period minimums are enforced to not be below the original + # minimum that occurred across the whole longer period. if self.accumulation: self._calculate_accumulation( cube_t0, period_reference, interpolated_cube diff --git a/improver_tests/utilities/test_TemporalInterpolation.py b/improver_tests/utilities/test_TemporalInterpolation.py index f61ca9d6b3..6f9beb2ed2 100644 --- a/improver_tests/utilities/test_TemporalInterpolation.py +++ b/improver_tests/utilities/test_TemporalInterpolation.py @@ -125,7 +125,7 @@ def multi_time_cube( the cube. data: The data to be contained in the cube. If the cube is 3-D the - leading dimension should be the same size and the list of times + leading dimension should be the same size as the list of times and will be sliced to associate each slice with each time. spatial_grid: Whether this is a lat-lon or equal areas projection. @@ -612,7 +612,7 @@ def test_input_cubes_in_incorrect_time_order(): def test_input_cube_with_multiple_times(): """Test that an exception is raised if a cube is provided that has - multiple validity times, e.g. a multi-entried time dimension.""" + multiple validity times (a multi-entried time dimension).""" times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 6, 9]] data = np.ones((5, 5), dtype=np.float32) @@ -747,9 +747,9 @@ def test_add_bounds(input_times, expected_time_bounds, expected_fp_bounds): [3, 6], [2.5, 2.5], ), - # Equal adjacent e.g. period maxes, shorter periods have the same max. + # Equal adjacent period maxes, shorter periods have the same max. ({"interval_in_minutes": 180, "max": True}, [5, 5], [3, 6], [5, 5],), - # Equal adjacent e.g. period minimums, shorter periods have the same + # Equal adjacent period minimums, shorter periods have the same # min. ({"interval_in_minutes": 180, "min": True}, [5, 5], [3, 6], [5, 5],), # Trend of increasing accumulations with time, which is reflected @@ -797,7 +797,7 @@ def test_add_bounds(input_times, expected_time_bounds, expected_fp_bounds): [1, 6], [1, 5], ), - # Equal adjacent e.g. period maxes, unequal shorter periods have the + # Equal adjacent period maxes, unequal shorter periods have the # same max. ( {"times": [datetime.datetime(2017, 11, 1, 4)], "max": True}, @@ -805,7 +805,7 @@ def test_add_bounds(input_times, expected_time_bounds, expected_fp_bounds): [1, 6], [6, 6], ), - # Equal adjacent e.g. period minimums, unequal shorter periods have + # Equal adjacent period minimums, unequal shorter periods have # the same min. ( {"times": [datetime.datetime(2017, 11, 1, 4)], "min": True}, @@ -948,7 +948,7 @@ def test_process_accumulation_unequal_inputs( ): """Test that the expected values are returned when the accumulation inputs are of different periods. The accumulations are converted to rates using - each input cube's period prior to interpolation which allos for this.""" + each input cube's period prior to interpolation which allows for this.""" times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] npoints = 5 From 37a0fd20968bc8a4019434b530903d9417ec1212 Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Mon, 18 Mar 2024 17:04:43 +0000 Subject: [PATCH 18/19] Second review changes. --- improver/utilities/temporal_interpolation.py | 127 +++++++++++------- .../utilities/test_TemporalInterpolation.py | 43 ++++-- 2 files changed, 108 insertions(+), 62 deletions(-) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index 362e4a21ec..8084676619 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -72,25 +72,40 @@ class TemporalInterpolation(BasePlugin): If working with period maximums and minimums we cannot return values in the new periods that do not adhere to the inputs. For example, we might - have a 3-hour maximum of 5 ms-1. The period before it might have a - maximum of 10 ms-1. However, on splitting the 3-hour period into 1-hour - periods, we cannot use the gradient between the two original period - maximums to produce a wind speed maximum in any of the interpolated - e.g. 1-hour periods that is above 5 ms-1 as we already know this to be the - maximum over the full 3-hours. We could use the trend to produce a lower - maximum in the interpolated 1-hour periods as this is not an unreasonable - interpretation of the gradient information and is consistent with the - original period maximum. A similar argument can be made about modifying - the minimums. + have a 3-hour maximum of 5 ms-1 between 03-06Z. The period before it might + have a maximum of 11 ms-1. Upon splitting the 3-hour period into 1-hour + periods the gradient might give us the following results: + + Inputs: 00-03Z: 11 ms-1, 03-06Z: 5 ms-1 + Outputs: 03-04Z: 9 ms-1, 04-05Z: 7 ms-1, 05-06Z: 5ms-1 + + However these outputs are not in agreement with the original 3-hour period + maximum of 5 ms-1 over the period 03-06Z. We enforce the maximum from the + original period which results in: + + Inputs: 00-03Z: 10 ms-1, 03-06Z: 5 ms-1 + Outputs: 03-04Z: 5 ms-1, 04-05Z: 5 ms-1, 05-06Z: 5ms-1 + + If instead the preceding period maximum was 2 ms-1 we would use the trend + to produce lower maximums in the interpolated 1-hour periods, becoming: + + Inputs: 00-03Z: 2 ms-1, 03-06Z: 5 ms-1 + Outputs: 03-04Z: 3 ms-1, 04-05Z: 4 ms-1, 05-06Z: 5ms-1 + + This interpretation of the gradient information is retained in the output + as it is consistent with the original period maximum of 5 ms-1 between + 03-06Z. As such we can impart increasing trends into maximums over periods + but not decreasing trends. The counter argument can be made when + interpolating minimums in periods, allowing us only to introduce + decreasing trends for these. We could use the cell methods to determine whether we are working with accumulations, maximums, or minimums. This should be denoted as a cell method associated with the time coordinate, e.g. for an accumulation it would be `time: sum`, whilst a maximum would have `time: max`. However - we cannot guarantee these cell methods are present and not assumed. As - such the interpolation of periods here relies on the user supplying a - suitable keyword argument that denotes the type of period being - processed. + we cannot guarantee these cell methods are present. As such the + interpolation of periods here relies on the user supplying a suitable + keyword argument that denotes the type of period being processed. """ def __init__( @@ -142,6 +157,8 @@ def __init__( ValueError: If both interval_in_minutes and times are both set. ValueError: If interpolation method not in known list. ValueError: If multiple period diagnostic kwargs are set True. + ValueError: A period diagnostic is being interpolated with a method + not found in the period_interpolation_methods list. """ if interval_in_minutes is None and times is None: raise ValueError( @@ -158,7 +175,6 @@ def __init__( self.interval_in_minutes = interval_in_minutes self.times = times known_interpolation_methods = ["linear", "solar", "daynight"] - self.period_interpolation_methods = ["linear"] if interpolation_method not in known_interpolation_methods: raise ValueError( "TemporalInterpolation: Unknown interpolation " @@ -175,6 +191,16 @@ def __init__( self.accumulation = accumulation self.max = max self.min = min + if any([accumulation, max, min]): + self.period_inputs = True + + period_interpolation_methods = ["linear"] + if self.interpolation_method not in period_interpolation_methods: + raise ValueError( + "Period diagnostics can only be temporally interpolated " + f"using these methods: {period_interpolation_methods}.\n" + f"Currently selected method is: {self.interpolation_method}." + ) def construct_time_list( self, initial_time: datetime, final_time: datetime @@ -532,10 +558,10 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: Raises: TypeError: If cube_t0 and cube_t1 are not of type iris.cube.Cube. - ValueError: A period diagnostic is being interpolated with a method - not found in the period_interpolation_methods list. - ValueError: A period diagnostic is being processed without the - type of period being specified. + ValueError: A mix of instantaneous and period diagnostics have + been used as inputs. + ValueError: A period type has been declared but inputs are not + period diagnostics. ValueError: Period diagnostics with overlapping periods. CoordinateNotFoundError: The input cubes contain no time coordinate. @@ -553,36 +579,6 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: ) raise TypeError(msg) - if cube_t0.coord("time").bounds is not None: - if self.interpolation_method not in self.period_interpolation_methods: - raise ValueError( - "Period diagnostics can only be temporally interpolated " - f"using these methods: {self.period_interpolation_methods}.\n" - f"Currently selected method is: {self.interpolation_method}." - ) - if not any([self.accumulation, self.max, self.min]): - raise ValueError( - "A type of period must be specified when interpolating a " - "period diagnostic. This may be a period max, min" - " or an accumulation." - ) - cube_interval = ( - cube_t1.coord("time").points[0] - cube_t0.coord("time").points[0] - ) - (period,) = np.diff(cube_t1.coord("time").bounds[0]) - if not cube_interval == period: - raise ValueError( - "The diagnostic provided represents the period " - f"{period / 3600} hours. The interval between the " - f"diagnostics is {cube_interval / 3600} hours. Temporal " - "interpolation can only be applied to a period " - "diagnostic provided at intervals that match the " - "diagnostic period such that all points in time are " - "captured by only one of the inputs and do not overlap." - ) - - self.period_inputs = True - try: (initial_time,) = iris_time_to_datetime(cube_t0.coord("time")) (final_time,) = iris_time_to_datetime(cube_t1.coord("time")) @@ -604,6 +600,39 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: "time." ) + cube_t0_bounds = cube_t0.coord("time").has_bounds() + cube_t1_bounds = cube_t1.coord("time").has_bounds() + if cube_t0_bounds + cube_t1_bounds == 1: + raise ValueError( + "Period and non-period diagnostics cannot be combined for" + " temporal interpolation." + ) + + if self.period_inputs: + # Declaring period type requires the inputs be period diagnostics. + if not cube_t0_bounds: + raise ValueError( + "A period method has been declared for temporal " + "interpolation (max, min, or accumulation). Period " + "diagnostics must be provided. The input cubes have no " + "time bounds." + ) + + cube_interval = ( + cube_t1.coord("time").points[0] - cube_t0.coord("time").points[0] + ) + (period,) = np.diff(cube_t1.coord("time").bounds[0]) + if not cube_interval == period: + raise ValueError( + "The diagnostic provided represents the period " + f"{period / 3600} hours. The interval between the " + f"diagnostics is {cube_interval / 3600} hours. Temporal " + "interpolation can only be applied to a period " + "diagnostic provided at intervals that match the " + "diagnostic period such that all points in time are " + "captured by only one of the inputs and do not overlap." + ) + time_list = self.construct_time_list(initial_time, final_time) # If the target output time is the same as the time at which the diff --git a/improver_tests/utilities/test_TemporalInterpolation.py b/improver_tests/utilities/test_TemporalInterpolation.py index 6f9beb2ed2..43fc3eb3ba 100644 --- a/improver_tests/utilities/test_TemporalInterpolation.py +++ b/improver_tests/utilities/test_TemporalInterpolation.py @@ -260,6 +260,10 @@ def daynight_mask(): {"interval_in_minutes": 60, "max": True, "min": True}, "Only one type of period diagnostics may be specified:", ), # Invalid interpolation method requested + ( + {"interval_in_minutes": 60, "max": True, "interpolation_method": "solar"}, + "Period diagnostics can only be temporally interpolated", + ), # Invalid interpolation method requested ], ) def test__init__(kwargs, exception): @@ -639,32 +643,45 @@ def test_input_cubelists_raises_exception(): TemporalInterpolation(interval_in_minutes=180).process(cubes, cube[1]) -def test_invalid_method_for_period_exception(): - """Test that providing a period diagnostic and attempting to apply an - unsuitable interpolation method raises an exception.""" +def test_mix_instantaneous_and_period(): + """Test that providing one instantaneous and one period diagnostic raises + an exception.""" times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] data = np.ones((5, 5), dtype=np.float32) cube = multi_time_cube(times, data, "latlon", bounds=True) - msg = "Period diagnostics can only be temporally interpolated" + cube_0 = cube[0] + cube_1 = cube[1] + for crd in ["time", "forecast_period"]: + cube_0.coord(crd).bounds = None + + msg = "Period and non-period diagnostics cannot be combined" with pytest.raises(ValueError, match=msg): - TemporalInterpolation( - interval_in_minutes=180, interpolation_method="solar" - ).process(cube[0], cube[1]) + TemporalInterpolation(interval_in_minutes=180).process(cube_0, cube_1) -def test_period_without_chosen_type_exception(): - """Test that providing a period diagnostic but not specifying a type from - the available min, max, or accumulation raises an exception.""" +@pytest.mark.parametrize( + "kwargs", + ( + [ + {"interval_in_minutes": 180, "accumulation": True}, + {"interval_in_minutes": 180, "max": True}, + {"interval_in_minutes": 180, "min": True}, + ] + ), +) +def test_period_method_non_period_diagnostics(kwargs): + """Test that declaring a period type for the interpolation and then + passing in non-period diagnostics raises an exception.""" times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] data = np.ones((5, 5), dtype=np.float32) - cube = multi_time_cube(times, data, "latlon", bounds=True) + cube = multi_time_cube(times, data, "latlon") - msg = "A type of period must be specified when interpolating" + msg = "A period method has been declared for temporal" with pytest.raises(ValueError, match=msg): - TemporalInterpolation(interval_in_minutes=180).process(cube[0], cube[1]) + TemporalInterpolation(**kwargs).process(cube[0], cube[1]) def test_period_unequal_to_interval_exception(): From 85694875d0bb22f2eeca8be9fc6105affaf21f37 Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Tue, 19 Mar 2024 16:36:50 +0000 Subject: [PATCH 19/19] Additional check within plugin. Added bounds checking to unit tests. --- improver/utilities/temporal_interpolation.py | 6 +++ .../utilities/test_TemporalInterpolation.py | 37 +++++++++++++++++++ 2 files changed, 43 insertions(+) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index 8084676619..c747ae5ee7 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -608,6 +608,12 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: " temporal interpolation." ) + if cube_t0_bounds and not self.period_inputs: + raise ValueError( + "Interpolation of period diagnostics should be done using " + "the appropriate period specifier (accumulation, min or max)." + ) + if self.period_inputs: # Declaring period type requires the inputs be period diagnostics. if not cube_t0_bounds: diff --git a/improver_tests/utilities/test_TemporalInterpolation.py b/improver_tests/utilities/test_TemporalInterpolation.py index 43fc3eb3ba..be00c42841 100644 --- a/improver_tests/utilities/test_TemporalInterpolation.py +++ b/improver_tests/utilities/test_TemporalInterpolation.py @@ -661,6 +661,19 @@ def test_mix_instantaneous_and_period(): TemporalInterpolation(interval_in_minutes=180).process(cube_0, cube_1) +def test_period_diagnostic_no_period_type(): + """Test that providing a period diagnostic without declaring the type of + period to be processed raises an exception.""" + + times = [datetime.datetime(2017, 11, 1, hour) for hour in [3, 9]] + data = np.ones((5, 5), dtype=np.float32) + cube = multi_time_cube(times, data, "latlon", bounds=True) + + msg = "Interpolation of period diagnostics should be done using" + with pytest.raises(ValueError, match=msg): + TemporalInterpolation(interval_in_minutes=180).process(cube[0], cube[1]) + + @pytest.mark.parametrize( "kwargs", ( @@ -872,10 +885,22 @@ def test_process_periods(kwargs, values, offsets, expected, realizations): else: expected_data = np.full((npoints, npoints), value) expected_time = 1509505200 + (offset * 3600) + expected_lower_bound_time = 1509505200 + [0, *offsets][i] * 3600 + expected_upper_bound_time = expected_time expected_fp = (6 + offset) * 3600 + expected_lower_bound_fp = (6 + [0, *offsets][i]) * 3600 + expected_upper_bound_fp = expected_fp assert result[i].coord("time").points[0] == expected_time + np.testing.assert_array_equal( + result[i].coord("time").bounds, + [[expected_lower_bound_time, expected_upper_bound_time]], + ) assert result[i].coord("forecast_period").points[0] == expected_fp + np.testing.assert_array_equal( + result[i].coord("forecast_period").bounds, + [[expected_lower_bound_fp, expected_upper_bound_fp]], + ) assert result[i].coord("time").points.dtype == "int64" assert result[i].coord("forecast_period").points.dtype == "int32" if value is not None: @@ -994,10 +1019,22 @@ def test_process_accumulation_unequal_inputs( else: expected_data = np.full((npoints, npoints), value) expected_time = 1509505200 + (offset * 3600) + expected_lower_bound_time = 1509505200 + [0, *offsets][i] * 3600 + expected_upper_bound_time = expected_time expected_fp = (6 + offset) * 3600 + expected_lower_bound_fp = (6 + [0, *offsets][i]) * 3600 + expected_upper_bound_fp = expected_fp assert result[i].coord("time").points[0] == expected_time + np.testing.assert_array_equal( + result[i].coord("time").bounds, + [[expected_lower_bound_time, expected_upper_bound_time]], + ) assert result[i].coord("forecast_period").points[0] == expected_fp + np.testing.assert_array_equal( + result[i].coord("forecast_period").bounds, + [[expected_lower_bound_fp, expected_upper_bound_fp]], + ) assert result[i].coord("time").points.dtype == "int64" assert result[i].coord("forecast_period").points.dtype == "int32" if value is not None: