Tests for rotate_grid_vectors

lib/iris/tests/__init__.py

@@ -642,7 +642,7 @@ def assertMaskedArrayAlmostEqual(self, a, b, decimal=6, strict=False):
         self._assertMaskedArray(np.testing.assert_array_almost_equal, a, b,
                                 strict, decimal=decimal)
 
-    def assertArrayAllClose(self, a, b, rtol=1.0e-7, atol=0.0, **kwargs):
+    def assertArrayAllClose(self, a, b, rtol=1.0e-7, atol=1.0e-8, **kwargs):
         """
         Check arrays are equal, within given relative + absolute tolerances.
 
@@ -656,20 +656,22 @@ def assertArrayAllClose(self, a, b, rtol=1.0e-7, atol=0.0, **kwargs):
         * rtol, atol (float):
             Relative and absolute tolerances to apply.
 
-        Any additional kwargs are passed to numpy.testing.assert_allclose.
+        Any additional kwargs are passed to numpy.allclose.
 
         Performs pointwise toleranced comparison, and raises an assertion if
         the two are not equal 'near enough'.
-        For full details see underlying routine numpy.testing.assert_allclose.
+        For full details see underlying routine numpy.allclose.
 
         """
-        ok = np.allclose(a, b, atol=atol, rtol=rtol, equal_nan=True)
+        ok = np.allclose(a, b, rtol=rtol, atol=atol, **kwargs)
         if not ok:
             # Calculate errors above a pointwise tolerance : The method is
             # taken from "numpy.core.numeric.isclose".
+            a, b = np.broadcast_arrays(a, b)
             errors = (np.abs(a-b) - atol + rtol * np.abs(b))
             worst_inds = np.unravel_index(np.argmax(errors.flat), errors.shape)
-            # Build a more useful message than from np.testing.assert_allclose.
+
+            # Also raise a more useful message than np.testing.assert_allclose.
             msg = ('\nARRAY CHECK FAILED "assertArrayAllClose" :'
                    '\n  with shapes={} {}, atol={}, rtol={}'
                    '\n  worst at element {} :  a={}  b={}'

lib/iris/tests/unit/analysis/cartography/test_rotate_grid_vectors.py

@@ -0,0 +1,150 @@
+# (C) British Crown Copyright 2018, Met Office
+#
+# This file is part of Iris.
+#
+# Iris is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# Iris is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with Iris.  If not, see <http://www.gnu.org/licenses/>.
+"""
+Unit tests for the function
+:func:`iris.analysis.cartography.gridcell_angles`.
+
+"""
+from __future__ import (absolute_import, division, print_function)
+from six.moves import (filter, input, map, range, zip)  # noqa
+
+# Import iris.tests first so that some things can be initialised before
+# importing anything else.
+import iris.tests as tests
+
+from mock import Mock
+from mock import call as mock_call
+import numpy as np
+
+from iris.cube import Cube
+from iris.coords import AuxCoord
+
+from iris.analysis.cartography import rotate_grid_vectors
+from iris.tests.stock import sample_2d_latlons
+
+
+class TestRotateGridVectors(tests.IrisTest):
+    def _check_angles_calculation(self, angles_in_degrees=True,
+                                  nan_angles_mask=None):
+        # Check basic maths on a 2d latlon grid.
+        u_cube = sample_2d_latlons(regional=True, transformed=True)
+        u_cube.units = 'ms-1'
+        u_cube.rename('dx')
+        u_cube.data[...] = 0
+        v_cube = u_cube.copy()
+        v_cube.name('dy')
+
+        # Define 6 different vectors, repeated in each data row.
+        in_vu = np.array([(0, 1), (2, -1), (-1, -1), (-3, 1), (2, 0), (0, 0)])
+        in_angs = np.rad2deg(np.arctan2(in_vu[..., 0], in_vu[..., 1]))
+        in_mags = np.sqrt(np.sum(in_vu * in_vu, axis=1))
+        v_cube.data[...] = in_vu[..., 0]
+        u_cube.data[...] = in_vu[..., 1]
+
+        # Define 5 different test rotation angles, one for each data row.
+        rotation_angles = np.array([0., -45., 135, -140., 90.])
+        ang_cube_data = np.broadcast_to(rotation_angles[:, None],
+                                        u_cube.shape)
+        ang_cube = u_cube.copy()
+        if angles_in_degrees:
+            ang_cube.units = 'degrees'
+        else:
+            ang_cube.units = 'radians'
+            ang_cube_data = np.deg2rad(ang_cube_data)
+        ang_cube.data[:] = ang_cube_data
+
+        if nan_angles_mask is not None:
+            ang_cube.data[nan_angles_mask] = np.nan
+
+        # Rotate all vectors by all the given angles.
+        result = rotate_grid_vectors(u_cube, v_cube, ang_cube)
+        out_u, out_v = [cube.data for cube in result]
+
+        # Check that vector magnitudes were unchanged.
+        out_mags = np.sqrt(out_u * out_u + out_v * out_v)
+        expect_mags = in_mags[None, :]
+        self.assertArrayAllClose(out_mags, expect_mags)
+
+        # Check that vector angles are all as expected.
+        out_angs = np.rad2deg(np.arctan2(out_v, out_u))
+        expect_angs = in_angs[None, :] + rotation_angles[:, None]
+        ang_diffs = out_angs - expect_angs
+        # Fix for null vectors, and +/-360 differences.
+        ang_diffs[np.abs(out_mags) < 0.001] = 0.0
+        ang_diffs = ang_diffs % 360.0
+        # Check that any differences are very small.
+        self.assertArrayAllClose(ang_diffs, 0.0)
+
+        # Check that results are always masked arrays, masked at NaN angles.
+        self.assertTrue(np.ma.isMaskedArray(out_u))
+        self.assertTrue(np.ma.isMaskedArray(out_v))
+        if nan_angles_mask is not None:
+            self.assertArrayEqual(out_u.mask, nan_angles_mask)
+            self.assertArrayEqual(out_v.mask, nan_angles_mask)
+
+    def test_angles_calculation(self):
+        self._check_angles_calculation()
+
+    def test_angles_in_radians(self):
+        self._check_angles_calculation(angles_in_degrees=True)
+
+    def test_angles_from_grid(self):
+        # Check it will gets angles from 'u_cube', and pass any kwargs on to
+        # the angles routine.
+        u_cube = sample_2d_latlons(regional=True, transformed=True)
+        u_cube = u_cube[:2, :3]
+        u_cube.units = 'ms-1'
+        u_cube.rename('dx')
+        u_cube.data[...] = 1.0
+        v_cube = u_cube.copy()
+        v_cube.name('dy')
+        v_cube.data[...] = 0.0
+
+        # Setup a fake angles result from the inner call to 'gridcell_angles'.
+        angles_result_data = np.array([[0.0, 90.0, 180.0],
+                                       [-180.0, -90.0, 270.0]])
+        angles_result_cube = Cube(angles_result_data, units='degrees')
+        angles_kwargs = {'this': 2}
+        angles_call_patch = self.patch(
+                'iris.analysis._grid_angles.gridcell_angles',
+                Mock(return_value=angles_result_cube))
+
+        # Call the routine.
+        result = rotate_grid_vectors(u_cube, v_cube,
+                                     grid_angles_kwargs=angles_kwargs)
+
+        self.assertEqual(angles_call_patch.call_args_list,
+                         [mock_call(u_cube, this=2)])
+
+        out_u, out_v = [cube.data for cube in result]
+        # Records what results should be for the various n*90deg rotations.
+        expect_u = np.array([[1.0, 0.0, -1.0],
+                             [-1.0, 0.0, 0.0]])
+        expect_v = np.array([[0.0, 1.0, 0.0],
+                             [0.0, -1.0, -1.0]])
+        # Check results are as expected.
+        self.assertArrayAllClose(out_u, expect_u)
+        self.assertArrayAllClose(out_v, expect_v)
+
+    def test_nan_vectors(self):
+        bad_angle_points = np.zeros((5, 6), dtype=bool)
+        bad_angle_points[2,3] = True
+        self._check_angles_calculation(nan_angles_mask=bad_angle_points)
+
+
+if __name__ == "__main__":
+    tests.main()