added masking of erroneous values

lib/iris/analysis/cartography.py

@@ -858,6 +858,45 @@ def _transform_distance_vectors(src_crs, x, y, u_dist, v_dist, tgt_crs):
     return u2_dist, v2_dist
 
 
+def _transform_distance_vectors_tolerance_mask(src_crs, x, y, tgt_crs):
+    """
+    Return a mask that can be applied to data array to mask elements
+    where the magnitude of vectors are not preserved due to numerical
+    errors introduced by the tranformation between coordinate systems.
+
+    Args:
+    * src_crs (`cartopy.crs.Projection`):
+        The source coordinate reference systems.
+    * x, y (array):
+        Locations of each vector defined in 'src_crs'.
+    * tgt_crs (`cartopy.crs.Projection`):
+        The target coordinate reference systems.
+
+    Returns:
+        2d boolean array that is the same shape as x and y.
+
+    """
+    if x.shape != y.shape:
+        raise ValueError('Arrays do not have matching shapes. '
+                         'x.shape is {}, y.shape is {}.'.format(x.shape,
+                                                                y.shape))
+    ones = np.ones(x.shape)
+    zeros = np.zeros(x.shape)
+    u_one_t, v_zero_t = _transform_distance_vectors(src_crs, x, y,
+                                                    ones, zeros, tgt_crs)
+    u_zero_t, v_one_t = _transform_distance_vectors(src_crs, x, y,
+                                                    zeros, ones, tgt_crs)
+    # Squared magnitudes should be equal to one within acceptable tolerance.
+    # A value of atol=2e-3 is used, which corresponds to a change in magnitude
+    # of approximately 0.1%.
+    sqmag_1_0 = u_one_t**2 + v_zero_t**2
+    sqmag_0_1 = u_zero_t**2 + v_one_t**2
+    mask = np.logical_not(
+        np.logical_and(np.isclose(sqmag_1_0, ones, atol=2e-3),
+                       np.isclose(sqmag_0_1, ones, atol=2e-3)))
+    return mask
+
+
 def rotate_winds(u_cube, v_cube, target_cs):
     """
     Transform wind vectors to a different coordinate system.
@@ -902,12 +941,13 @@ def rotate_winds(u_cube, v_cube, target_cs):
         The names of the output cubes are those of the inputs, prefixed with
         'transformed\_' (e.g. 'transformed_x_wind').
 
-    .. note::
+    .. warning::
 
         Conversion between rotated-pole and non-rotated systems can be
         expressed analytically.  However, this function always uses a numerical
-        approach.
-
+        approach. In locations where this numerical approach does not preserve
+        magnitude to an accuracy of 0.1%, the corresponding elements of the
+        returned cubes will be masked.
 
     """
     # Check u_cube and v_cube have the same shape. We iterate through
@@ -975,12 +1015,27 @@ def rotate_winds(u_cube, v_cube, target_cs):
     else:
         dims = x_dims
 
+    # Transpose x, y 2d arrays to match the order in cube's data
+    # array so that x, y and the sliced data all line up.
+    if dims[0] > dims[1]:
+        x = x.transpose()
+        y = y.transpose()
+
     # Create resulting cubes.
     ut_cube = u_cube.copy()
     vt_cube = v_cube.copy()
     ut_cube.rename('transformed_{}'.format(u_cube.name()))
     vt_cube.rename('transformed_{}'.format(v_cube.name()))
 
+    # Calculate mask based on preservation of magnitude.
+    mask = _transform_distance_vectors_tolerance_mask(src_crs, x, y,
+                                                      target_crs)
+    apply_mask = mask.any()
+    if apply_mask:
+        # Make masked arrays to accept masking.
+        ut_cube.data = ma.asanyarray(ut_cube.data)
+        vt_cube.data = ma.asanyarray(vt_cube.data)
+
     # Project vectors with u, v components one horiz slice at a time and
     # insert into the resulting cubes.
     shape = list(u_cube.shape)
@@ -995,13 +1050,26 @@ def rotate_winds(u_cube, v_cube, target_cs):
         u = u_cube.data[index]
         v = v_cube.data[index]
         ut, vt = _transform_distance_vectors(src_crs, x, y, u, v, target_crs)
+        if apply_mask:
+            ut = ma.asanyarray(ut)
+            ut[mask] = ma.masked
+            vt = ma.asanyarray(vt)
+            vt[mask] = ma.masked
         ut_cube.data[index] = ut
         vt_cube.data[index] = vt
 
     # Calculate new coords of locations in target coordinate system.
     xyz_tran = target_crs.transform_points(src_crs, x, y)
     xt = xyz_tran[..., 0].reshape(x.shape)
     yt = xyz_tran[..., 1].reshape(y.shape)
+
+    # Transpose xt, yt 2d arrays to match the dim order
+    # of the original x an y arrays - i.e. undo the earlier
+    # transpose (if applied).
+    if dims[0] > dims[1]:
+        xt = xt.transpose()
+        yt = yt.transpose()
+
     xt_coord = iris.coords.AuxCoord(xt,
                                     standard_name='projection_x_coordinate',
                                     coord_system=target_cs)

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

@@ -26,6 +26,7 @@
 import iris.tests as tests
 
 import numpy as np
+import numpy.ma as ma
 
 import cartopy.crs as ccrs
 from iris.analysis.cartography import rotate_winds, unrotate_pole
@@ -261,6 +262,48 @@ def test_new_coords(self):
         self.assertEqual(ut.coord('projection_y_coordinate'), expected_y)
         self.assertEqual(vt.coord('projection_y_coordinate'), expected_y)
 
+    def test_new_coords_transposed(self):
+        u, v = self._uv_cubes_limited_extent()
+        # Transpose cubes so that cube is in xy order rather than the
+        # typical yx order of meshgrid.
+        u.transpose()
+        v.transpose()
+        x = u.coord('grid_longitude').points
+        y = u.coord('grid_latitude').points
+        x2d, y2d = np.meshgrid(x, y)
+        src_crs = ccrs.RotatedPole(pole_longitude=177.5, pole_latitude=37.5)
+        tgt_crs = ccrs.OSGB()
+        xyz_tran = tgt_crs.transform_points(src_crs, x2d, y2d)
+
+        ut, vt = rotate_winds(u, v, iris.coord_systems.OSGB())
+
+        points = xyz_tran[..., 0].reshape(x2d.shape)
+        expected_x = AuxCoord(points,
+                              standard_name='projection_x_coordinate',
+                              units='m',
+                              coord_system=iris.coord_systems.OSGB())
+        self.assertEqual(ut.coord('projection_x_coordinate'), expected_x)
+        self.assertEqual(vt.coord('projection_x_coordinate'), expected_x)
+
+        points = xyz_tran[..., 1].reshape(y2d.shape)
+        expected_y = AuxCoord(points,
+                              standard_name='projection_y_coordinate',
+                              units='m',
+                              coord_system=iris.coord_systems.OSGB())
+        self.assertEqual(ut.coord('projection_y_coordinate'), expected_y)
+        self.assertEqual(vt.coord('projection_y_coordinate'), expected_y)
+        # Check dim mapping for 2d coords is yx.
+        expected_dims = (u.coord_dims('grid_latitude') +
+                         u.coord_dims('grid_longitude'))
+        self.assertEqual(ut.coord_dims('projection_x_coordinate'),
+                         expected_dims)
+        self.assertEqual(ut.coord_dims('projection_y_coordinate'),
+                         expected_dims)
+        self.assertEqual(vt.coord_dims('projection_x_coordinate'),
+                         expected_dims)
+        self.assertEqual(vt.coord_dims('projection_y_coordinate'),
+                         expected_dims)
+
     def test_orig_coords(self):
         u, v = self._uv_cubes_limited_extent()
         ut, vt = rotate_winds(u, v, iris.coord_systems.OSGB())
@@ -310,6 +353,66 @@ def test_nd_data(self):
         self.assertArrayAlmostEqual(ut.data, expected_ut_data)
         self.assertArrayAlmostEqual(vt.data, expected_vt_data)
 
+    def test_transposed(self):
+        # Test case where the coordinates are not ordered yx in the cube.
+        u, v = self._uv_cubes_limited_extent()
+        # Slice out 4 points that lie in and outside OSGB extent.
+        u = u[1:3, 3:5]
+        v = v[1:3, 3:5]
+        # Transpose cubes (in-place)
+        u.transpose()
+        v.transpose()
+        ut, vt = rotate_winds(u, v, iris.coord_systems.OSGB())
+        # Values precalculated and checked.
+        expected_ut_data = np.array([[0.16285514,  0.35323639],
+                                     [1.82650698,  2.62455840]]).T
+        expected_vt_data = np.array([[19.88979966,  19.01921346],
+                                     [19.88018847,  19.01424281]]).T
+        # Compare u and v data values against previously calculated values.
+        self.assertArrayAllClose(ut.data, expected_ut_data, rtol=1e-5)
+        self.assertArrayAllClose(vt.data, expected_vt_data, rtol=1e-5)
+
+
+class TestMasking(tests.IrisTest):
+    def test_rotated_to_osgb(self):
+        # Rotated Pole data with large extent.
+        x = np.linspace(311.9, 391.1, 10)
+        y = np.linspace(-23.6, 24.8, 8)
+        u, v = uv_cubes(x, y)
+        ut, vt = rotate_winds(u, v, iris.coord_systems.OSGB())
+
+        # Ensure cells with discrepancies in magnitude are masked.
+        self.assertTrue(ma.isMaskedArray(ut.data))
+        self.assertTrue(ma.isMaskedArray(vt.data))
+
+        # Snapshot of mask with fixed tolerance of atol=2e-3
+        expected_mask = np.array([[1, 1, 1, 0, 0, 0, 0, 0, 0, 1],
+                                  [1, 1, 1, 0, 0, 0, 0, 0, 0, 1],
+                                  [1, 1, 1, 1, 0, 0, 0, 0, 1, 1],
+                                  [1, 1, 1, 1, 0, 0, 0, 0, 1, 1],
+                                  [1, 1, 1, 1, 0, 0, 0, 0, 1, 1],
+                                  [1, 1, 1, 1, 1, 0, 0, 1, 1, 1],
+                                  [1, 1, 1, 1, 1, 0, 0, 1, 1, 1],
+                                  [1, 1, 1, 1, 1, 0, 0, 1, 1, 1]], np.bool)
+        self.assertArrayEqual(expected_mask, ut.data.mask)
+        self.assertArrayEqual(expected_mask, vt.data.mask)
+
+        # Check unmasked values have sufficiently small error in mag.
+        expected_mag = np.sqrt(u.data**2 + v.data**2)
+        # Use underlying data to ignore mask in calculation.
+        res_mag = np.sqrt(ut.data.data**2 + vt.data.data**2)
+        # Calculate percentage error (note there are no zero magnitudes
+        # so we can divide safely).
+        anom = 100.0 * np.abs(res_mag - expected_mag) / expected_mag
+        self.assertTrue(anom[~ut.data.mask].max() < 0.1)
+
+    def test_rotated_to_unrotated(self):
+        # Suffiently accurate so that no mask is introduced.
+        u, v = uv_cubes()
+        ut, vt = rotate_winds(u, v, iris.coord_systems.GeogCS(6371229))
+        self.assertFalse(ma.isMaskedArray(ut.data))
+        self.assertFalse(ma.isMaskedArray(vt.data))
+
 
 class TestRoundTrip(tests.IrisTest):
     def test_rotated_to_unrotated(self):