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 <[email protected]>
Co-authored-by: Belinda Trotta <[email protected]>

.mailmap

@@ -31,7 +31,7 @@ Paul Abernethy <[email protected]> <[email protected]
 Shafiat Dewan <[email protected]> <[email protected]>
 Shubhendra Singh Chauhan <[email protected]> <[email protected]>
 Simon Jackson <[email protected]> <[email protected]>
-Tim Hume <[email protected]> <[email protected]>
+Timothy Hume <[email protected]> <[email protected]>
 Tomasz Trzeciak <[email protected]> <[email protected]>
 Tomasz Trzeciak <[email protected]> <[email protected]>
 Tom Gale <[email protected]> <[email protected]>

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)

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)

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":

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)

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.

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