From 1c3b584e17557d48a88d01f8e4c589377839f2c6 Mon Sep 17 00:00:00 2001 From: gavinevans Date: Mon, 16 Oct 2023 23:27:21 +0100 Subject: [PATCH] Wind direction changes (#1958) * Simplify how angles are converted from (-180, 180] to [0, 360) The simpler method takes about 65% the CPU time as the original method (tested on a massive array of angles). Furthermore, the original method introduces a tiny error for a very few angles near 0/360 degrees. * Implement "proper" blending of wind directions. We blend wind direction data using the IMPROVER wind direction functionality. We identify wind direction data by it having units of degrees. This identification of wind directions could be generalised to any variabe which has units that can be converted to degrees, not just degrees. For example, if the units were radians. But I've not seen things like radians used in meteorology, so just restricting the code to variables with units of "degrees" is OK for now. * Add proper handling of wind direction to temporal interpolation. * Fix a temporal interpolation bug. This was introduced by me at the previous commit. * Some tests that wind direction interpolation works. Test 1: 350->20 degrees, expect 5 degrees. Test 2: 40->60 degrees, expect 50 degrees. * Put in unit tests to test blending of wind directions We assume that wind directions can be identified by the units attribute being "degrees". If this is the case, we use the special wind direction blending functionality to handle blending. * Move direction blending into the correct place. Initially we had direction blending in the wrong function. It worked, as long as that function was called. If it wasn't then direction blending didn't happen. * Change contributor name from Tim to Timothy. * Remove extra line. * Fix .mailmap * Some suggestions from the reviewer. * For some reason bilinear.py got messed around. * Covert to float32 * simplify code * Add keyword. --------- Co-authored-by: Timothy Hume Co-authored-by: Belinda Trotta --- .mailmap | 2 +- CONTRIBUTING.md | 2 +- improver/blending/weighted_blend.py | 15 ++++- improver/utilities/temporal_interpolation.py | 15 +++++ improver/wind_calculations/wind_direction.py | 7 +- .../test_WeightedBlendAcrossWholeDimension.py | 50 +++++++++++++++ .../utilities/test_TemporalInterpolation.py | 64 +++++++++++++++++++ 7 files changed, 146 insertions(+), 9 deletions(-) diff --git a/.mailmap b/.mailmap index 4d0a95f97d..4e357a68bb 100644 --- a/.mailmap +++ b/.mailmap @@ -31,7 +31,7 @@ Paul Abernethy <87321907+ShafiatDewan@users.noreply.github.com> Shubhendra Singh Chauhan Simon Jackson -Tim Hume <57921955+thbom001@users.noreply.github.com> <57921955+thbom001@users.noreply.github.com> +Timothy Hume <57921955+thbom001@users.noreply.github.com> <57921955+thbom001@users.noreply.github.com> Tomasz Trzeciak Tomasz Trzeciak Tom Gale diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 8a3740c2c8..7f41be145f 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -51,7 +51,7 @@ below: - Ben Hooper (Met Office, UK) - Aaron Hopkinson (Met Office, UK) - Kathryn Howard (Met Office, UK) - - Tim Hume (Bureau of Meteorology, Australia) + - Timothy Hume (Bureau of Meteorology, Australia) - Katharine Hurst (Met Office, UK) - Simon Jackson (Met Office, UK) - Caroline Jones (Met Office, UK) diff --git a/improver/blending/weighted_blend.py b/improver/blending/weighted_blend.py index a01da70f0f..2fe4f17348 100644 --- a/improver/blending/weighted_blend.py +++ b/improver/blending/weighted_blend.py @@ -55,6 +55,7 @@ get_dim_coord_names, sort_coord_in_cube, ) +from improver.wind_calculations.wind_direction import WindDirection class MergeCubesForWeightedBlending(BasePlugin): @@ -640,7 +641,8 @@ def percentile_weighted_mean(self, cube: Cube, weights: Optional[Cube]) -> Cube: def weighted_mean(self, cube: Cube, weights: Optional[Cube]) -> Cube: """ - Blend data using a weighted mean using the weights provided. + Blend data using a weighted mean using the weights provided. Circular + data identified with a unit of "degrees" are blended appropriately. Args: cube: @@ -653,6 +655,11 @@ def weighted_mean(self, cube: Cube, weights: Optional[Cube]) -> Cube: The cube with values blended over self.blend_coord, with suitable weightings applied. """ + + # If units are degrees, convert degrees to complex numbers. + if cube.units == "degrees": + cube.data = WindDirection.deg_to_complex(cube.data) + weights_array = self.get_weights_array(cube, weights) slice_dim = 1 @@ -676,6 +683,10 @@ def weighted_mean(self, cube: Cube, weights: Optional[Cube]) -> Cube: result = result_slices.merge_cube() if allow_slicing else result_slices[0] + # If units are degrees, convert complex numbers back to degrees. + if cube.units == "degrees": + result.data = WindDirection.complex_to_deg(result.data) + return result def process(self, cube: Cube, weights: Optional[Cube] = None) -> Cube: @@ -732,7 +743,7 @@ def process(self, cube: Cube, weights: Optional[Cube] = None) -> Cube: # Check that the time coordinate is single valued if required. self.check_compatible_time_points(cube) - # Do blending and update metadata + # Do blending and update metadata. if self.check_percentile_coord(cube): enforce_coordinate_ordering(cube, [self.blend_coord, "percentile"]) result = self.percentile_weighted_mean(cube, weights) diff --git a/improver/utilities/temporal_interpolation.py b/improver/utilities/temporal_interpolation.py index 9f8e2d2bc1..d93220a153 100644 --- a/improver/utilities/temporal_interpolation.py +++ b/improver/utilities/temporal_interpolation.py @@ -46,6 +46,7 @@ from improver.utilities.solar import DayNightMask, calc_solar_elevation from improver.utilities.spatial import lat_lon_determine, transform_grid_to_lat_lon from improver.utilities.temporal import iris_time_to_datetime +from improver.wind_calculations.wind_direction import WindDirection class TemporalInterpolation(BasePlugin): @@ -411,10 +412,24 @@ def process(self, cube_t0: Cube, cube_t1: Cube) -> CubeList: ) time_list = self.construct_time_list(initial_time, final_time) + + # 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 + # fashion. + if cube_t0.units == "degrees" and cube_t1.units == "degrees": + cube_t0.data = WindDirection.deg_to_complex(cube_t0.data) + cube_t1.data = WindDirection.deg_to_complex(cube_t1.data) + cubes = iris.cube.CubeList([cube_t0, cube_t1]) cube = MergeCubes()(cubes) interpolated_cube = cube.interpolate(time_list, iris.analysis.Linear()) + if cube_t0.units == "degrees" and cube_t1.units == "degrees": + interpolated_cube.data = WindDirection.complex_to_deg( + interpolated_cube.data + ) + self.enforce_time_coords_dtype(interpolated_cube) interpolated_cubes = iris.cube.CubeList() if self.interpolation_method == "solar": diff --git a/improver/wind_calculations/wind_direction.py b/improver/wind_calculations/wind_direction.py index 97da7c1e25..d92c927edc 100644 --- a/improver/wind_calculations/wind_direction.py +++ b/improver/wind_calculations/wind_direction.py @@ -184,11 +184,8 @@ def complex_to_deg(complex_in: ndarray) -> ndarray: angle = np.angle(complex_in, deg=True) - # If angle negative value - add 360 degrees. - angle = np.where(angle < 0, angle + 360, angle) - - # If angle == 360 - set to zero degrees. - angle = np.where(np.isclose(angle, 360, atol=0.001), 0.0, angle) + # Convert angles so they are in the range [0, 360) + angle = np.mod(np.float32(angle), 360) # We don't need 64 bit precision. angle = angle.astype(np.float32) diff --git a/improver_tests/blending/weighted_blend/test_WeightedBlendAcrossWholeDimension.py b/improver_tests/blending/weighted_blend/test_WeightedBlendAcrossWholeDimension.py index 4cf52ca27b..2a96d82987 100644 --- a/improver_tests/blending/weighted_blend/test_WeightedBlendAcrossWholeDimension.py +++ b/improver_tests/blending/weighted_blend/test_WeightedBlendAcrossWholeDimension.py @@ -474,6 +474,56 @@ def test_without_weights(self): self.assertIsInstance(result, iris.cube.Cube) self.assertArrayAlmostEqual(result.data, expected) + def test_wind_directions(self): + """Test function when a wind direction data cube is provided, and + the directions don't cross the 0/360° boundary.""" + frt_points = [ + datetime(2015, 11, 19, 0), + datetime(2015, 11, 19, 1), + ] + cube = set_up_variable_cube( + np.zeros((2, 2), dtype=np.float32), + name="wind_from_direction", + units="degrees", + time=datetime(2015, 11, 19, 2), + frt=datetime(2015, 11, 19, 0), + standard_grid_metadata="gl_det", + attributes={"title": "Operational ENGL Model Forecast"}, + ) + cube = add_coordinate( + cube, frt_points, "forecast_reference_time", is_datetime=True + ) + cube.data[0] = 10.0 + cube.data[1] = 30.0 + expected = np.full((2, 2), 20.0) + result = self.plugin.weighted_mean(cube, weights=None) + self.assertArrayAlmostEqual(result.data, expected, decimal=4) + + def test_wind_directions_over_north(self): + """Test function when a wind direction data cube is provided, and + the directions cross the 0/360° boundary.""" + frt_points = [ + datetime(2015, 11, 19, 0), + datetime(2015, 11, 19, 1), + ] + cube = set_up_variable_cube( + np.zeros((2, 2), dtype=np.float32), + name="wind_direction", + units="degrees", + time=datetime(2015, 11, 19, 2), + frt=datetime(2015, 11, 19, 0), + standard_grid_metadata="gl_det", + attributes={"title": "Operational ENGL Model Forecast"}, + ) + cube = add_coordinate( + cube, frt_points, "forecast_reference_time", is_datetime=True + ) + cube.data[0] = 350.0 + cube.data[1] = 30.0 + expected = np.full((2, 2), 10.0) + result = self.plugin.weighted_mean(cube, weights=None) + self.assertArrayAlmostEqual(result.data, expected, decimal=4) + def test_collapse_dims_with_weights(self): """Test function matches when the blend coordinate is first or second.""" # Create a new axis. diff --git a/improver_tests/utilities/test_TemporalInterpolation.py b/improver_tests/utilities/test_TemporalInterpolation.py index 0e21408e44..c272995483 100644 --- a/improver_tests/utilities/test_TemporalInterpolation.py +++ b/improver_tests/utilities/test_TemporalInterpolation.py @@ -676,6 +676,70 @@ def test_return_type(self): ) 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