From f9359f3593bdaac1c6753ceca65e40efadc43ff9 Mon Sep 17 00:00:00 2001 From: Ed Campbell Date: Fri, 6 Feb 2015 16:59:51 +0000 Subject: [PATCH 1/3] added masking of erroneous values --- lib/iris/analysis/cartography.py | 56 ++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/lib/iris/analysis/cartography.py b/lib/iris/analysis/cartography.py index 82035d850f..2e3ef2bd07 100644 --- a/lib/iris/analysis/cartography.py +++ b/lib/iris/analysis/cartography.py @@ -858,6 +858,43 @@ 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. + 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=1e-3), + np.isclose(sqmag_0_1, ones, atol=1e-3))) + return mask + + def rotate_winds(u_cube, v_cube, target_cs): """ Transform wind vectors to a different coordinate system. @@ -975,12 +1012,26 @@ def rotate_winds(u_cube, v_cube, target_cs): else: dims = x_dims + # Transpose x, y 2d arrays to match order in cube's data. + 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,6 +1046,11 @@ 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 From 706caccb7c8ba255589a52b92ffd604f4b82d056 Mon Sep 17 00:00:00 2001 From: Ed Campbell Date: Wed, 11 Feb 2015 08:20:24 +0000 Subject: [PATCH 2/3] review actions --- lib/iris/analysis/cartography.py | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/lib/iris/analysis/cartography.py b/lib/iris/analysis/cartography.py index 2e3ef2bd07..3b449fb70f 100644 --- a/lib/iris/analysis/cartography.py +++ b/lib/iris/analysis/cartography.py @@ -887,6 +887,8 @@ def _transform_distance_vectors_tolerance_mask(src_crs, x, y, 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=1e-3 is used, which corresponds to a change in magnitude + # of approximately 0.05%. 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( @@ -939,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.05%, the corresponding elements of the + returned cubes will be masked. """ # Check u_cube and v_cube have the same shape. We iterate through @@ -1012,7 +1015,8 @@ def rotate_winds(u_cube, v_cube, target_cs): else: dims = x_dims - # Transpose x, y 2d arrays to match order in cube's data. + # 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() @@ -1023,7 +1027,7 @@ def rotate_winds(u_cube, v_cube, target_cs): ut_cube.rename('transformed_{}'.format(u_cube.name())) vt_cube.rename('transformed_{}'.format(v_cube.name())) - # Calculate mask based on preservation of magnitude + # Calculate mask based on preservation of magnitude. mask = _transform_distance_vectors_tolerance_mask(src_crs, x, y, target_crs) apply_mask = mask.any() @@ -1058,6 +1062,14 @@ def rotate_winds(u_cube, v_cube, target_cs): 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) From 6588dd546f4e29ab7b1a35015dd9c3d82a4ad6b8 Mon Sep 17 00:00:00 2001 From: Ed Campbell Date: Wed, 11 Feb 2015 15:17:12 +0000 Subject: [PATCH 3/3] increased tolerance and added testing --- lib/iris/analysis/cartography.py | 10 +- .../analysis/cartography/test_rotate_winds.py | 103 ++++++++++++++++++ 2 files changed, 108 insertions(+), 5 deletions(-) diff --git a/lib/iris/analysis/cartography.py b/lib/iris/analysis/cartography.py index 3b449fb70f..a7447907ba 100644 --- a/lib/iris/analysis/cartography.py +++ b/lib/iris/analysis/cartography.py @@ -887,13 +887,13 @@ def _transform_distance_vectors_tolerance_mask(src_crs, x, y, 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=1e-3 is used, which corresponds to a change in magnitude - # of approximately 0.05%. + # 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=1e-3), - np.isclose(sqmag_0_1, ones, atol=1e-3))) + np.logical_and(np.isclose(sqmag_1_0, ones, atol=2e-3), + np.isclose(sqmag_0_1, ones, atol=2e-3))) return mask @@ -946,7 +946,7 @@ def rotate_winds(u_cube, v_cube, target_cs): Conversion between rotated-pole and non-rotated systems can be expressed analytically. However, this function always uses a numerical approach. In locations where this numerical approach does not preserve - magnitude to an accuracy of 0.05%, the corresponding elements of the + magnitude to an accuracy of 0.1%, the corresponding elements of the returned cubes will be masked. """ diff --git a/lib/iris/tests/unit/analysis/cartography/test_rotate_winds.py b/lib/iris/tests/unit/analysis/cartography/test_rotate_winds.py index 173a0cf6f1..ae109e896a 100644 --- a/lib/iris/tests/unit/analysis/cartography/test_rotate_winds.py +++ b/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):