Public contiguity checking

lib/iris/coords.py

@@ -1013,7 +1013,7 @@ def _sanity_check_bounds(self):
                              'are not defined for coordinates with more than '
                              '2 dimensions.'.format(self.name()))
 
-    def _discontiguity_in_bounds(self, rtol=1e-05, atol=1e-08):
+    def _discontiguity_in_bounds(self, rtol=1e-5, atol=1e-8):
         """
         Checks that the bounds of the coordinate are contiguous.
 
@@ -1064,9 +1064,14 @@ def mod360_adjust(compare_axis):
                         modification = (index.astype(int) * 360) * sign
                         upper_bounds -= modification
 
-                diffs_along_axis = np.abs(upper_bounds - lower_bounds)
-                contiguous_along_axis = np.allclose(upper_bounds, lower_bounds,
-                                                    rtol=rtol, atol=atol)
+                diffs_between_cells = np.abs(upper_bounds - lower_bounds)
+                cell_size = lower_bounds - upper_bounds
+                diffs_along_axis = diffs_between_cells > (atol +
+                                                          rtol * cell_size)
+
+                points_close_enough = diffs_along_axis <= (atol +
+                                                           rtol * cell_size)
+                contiguous_along_axis = np.all(points_close_enough)
                 return diffs_along_axis, contiguous_along_axis
 
             diffs_along_x, match_cell_x1 = mod360_adjust(compare_axis='x')

lib/iris/plot.py

@@ -104,7 +104,7 @@ def get_span(coord):
             raise ValueError('The coordinate {!r} has {} dimensions.'
                              'Cell-based plotting is only supported for'
                              'coordinates with one or two dimensions.'
-                             .format(coord.name()), len(span))
+                             .format(coord.name(), len(span)))
 
     # Check the combination of coordinates spans enough (ndims) data
     # dimensions.
@@ -271,7 +271,7 @@ def _invert_yaxis(v_coord, axes=None):
             axes.invert_yaxis()
 
 
-def _check_bounds_contiguity_and_mask(coord, data, atol=None):
+def _check_bounds_contiguity_and_mask(coord, data, atol=None, rtol=None):
     """
     Checks that any discontiguities in the bounds of the given coordinate only
     occur where the data is masked.
@@ -297,6 +297,7 @@ def _check_bounds_contiguity_and_mask(coord, data, atol=None):
             tolerance.
 
     """
+    kwargs = {}
     data_is_masked = hasattr(data, 'mask')
     if data_is_masked:
         # When checking the location of the discontiguities, we check against
@@ -317,15 +318,16 @@ def _check_bounds_contiguity_and_mask(coord, data, atol=None):
 
     elif coord.ndim == 2:
         if atol:
-            kwargs = {'atol': atol}
-        else:
-            kwargs = {}
+            kwargs['atol'] = atol
+        if rtol:
+            kwargs['rtol'] = rtol
         contiguous, diffs = coord._discontiguity_in_bounds(**kwargs)
 
         if not contiguous and data_is_masked:
             diffs_along_x, diffs_along_y = diffs
 
-            # Check along both dimensions.
+            # Check along both dimensions that any discontiguous
+            # points are correctly masked.
             not_masked_at_discontiguity_along_x = np.any(
                 np.logical_and(mask_invert[:, :-1], diffs_along_x))
 
@@ -341,12 +343,20 @@ def _check_bounds_contiguity_and_mask(coord, data, atol=None):
         if not data_is_masked:
             raise ValueError('The bounds of the {} coordinate are not '
                              'contiguous. Not able to create a suitable grid'
-                             'to plot.'.format(coord.name()))
+                             'to plot. You can use '
+                             'iris.util.find_discontiguities() to identify '
+                             'discontiguities in your x and y coordinate '
+                             'bounds arrays.'.format(coord.name()))
         if not_masked_at_discontiguity:
             raise ValueError('The bounds of the {} coordinate are not '
                              'contiguous and data is not masked where the '
                              'discontiguity occurs. Not able to create a '
-                             'suitable grid to plot.'.format(
+                             'suitable grid to plot. You can use '
+                             'iris.util.find_discontiguities() to identify '
+                             'discontiguities in your x and y coordinate '
+                             'bounds arrays, and then mask them with '
+                             'iris.util.mask_cube()'
+                             ''.format(
                               coord.name()))
 
 

lib/iris/tests/stock/_stock_2d_latlons.py

@@ -302,17 +302,17 @@ def sample_cube(xargs, yargs):
     return cube
 
 
-def make_bounds_discontiguous_at_point(cube, at_iy, at_ix):
+def make_bounds_discontiguous_at_point(cube, at_iy, at_ix, in_y=False):
     """
-    Meddle with the XY grid bounds of a cube to make the grid discontiguous.
+    Meddle with the XY grid bounds of a 2D cube to make the grid discontiguous.
 
     Changes the points and bounds of a single gridcell, so that it becomes
-    discontinuous with the next gridcell to its right.
-    Also masks the cube data at the given point.
+    discontinuous with an adjacent gridcell : either the one to its right, or
+    the one above it (if 'in_y' is True).
 
-    The cube must have bounded 2d 'x' and 'y' coordinates.
+    Also masks the cube data at the given point.
 
-    TODO: add a switch to make a discontiguity in the *y*-direction instead ?
+    The cube must be 2-dimensional and have bounded 2d 'x' and 'y' coordinates.
 
     """
     x_coord = cube.coord(axis='x')
@@ -326,11 +326,21 @@ def adjust_coord(coord):
         pts, bds = coord.points, coord.bounds
         # Fetch the 4 bounds (bottom-left, bottom-right, top-right, top-left)
         bds_bl, bds_br, bds_tr, bds_tl = bds[at_iy, at_ix]
-        # Make a discontinuity "at" (iy, ix), by moving the right-hand edge of
-        # the cell to the midpoint of the existing left+right bounds.
-        new_bds_br = 0.5 * (bds_bl + bds_br)
-        new_bds_tr = 0.5 * (bds_tl + bds_tr)
-        bds_br, bds_tr = new_bds_br, new_bds_tr
+        if not in_y:
+            # Make a discontinuity "at" (iy, ix), by moving the right-hand edge
+            # of the cell to the midpoint of the existing left+right bounds.
+            new_bds_br = 0.5 * (bds_bl + bds_br)
+            new_bds_tr = 0.5 * (bds_tl + bds_tr)
+            bds_br, bds_tr = new_bds_br, new_bds_tr
+        else:
+            # Same but in the 'grid y direction' :
+            # Make a discontinuity "at" (iy, ix), by moving the **top** edge of
+            # the cell to the midpoint of the existing **top+bottom** bounds.
+            new_bds_tl = 0.5 * (bds_bl + bds_tl)
+            new_bds_tr = 0.5 * (bds_br + bds_tr)
+            bds_tl, bds_tr = new_bds_tl, new_bds_tr
+
+        # Write in the new bounds (all 4 corners).
         bds[at_iy, at_ix] = [bds_bl, bds_br, bds_tr, bds_tl]
         # Also reset the cell midpoint to the middle of the 4 new corners,
         # in case having a midpoint outside the corners might cause a problem.
@@ -341,10 +351,15 @@ def adjust_coord(coord):
 
     adjust_coord(x_coord)
     adjust_coord(y_coord)
-    # Also mask the relevant data point.
-    data = cube.data  # N.B. fetch all the data.
-    if not ma.isMaskedArray(data):
-        # Promote to masked array, to avoid converting mask to NaN.
-        data = ma.masked_array(data)
-    data[at_iy, at_ix] = ma.masked
-    cube.data = data
+
+    # Check which dimensions are spanned by each coordinate.
+    for coord in (x_coord, y_coord):
+        span = set(cube.coord_dims(coord))
+        if not span:
+            msg = 'The coordinate {!r} doesn\'t span a data dimension.'
+            raise ValueError(msg.format(coord.name()))
+
+        masked_data = ma.masked_array(cube.data)
+        masked_data[at_iy, at_ix] = ma.masked
+
+    cube.data = masked_data

lib/iris/tests/unit/coords/test_Coord.py

@@ -480,7 +480,8 @@ def test_2d_discontiguous_along_x(self):
         contiguous, diffs = coord._discontiguity_in_bounds()
         diffs_along_x, diffs_along_y = diffs
         self.assertFalse(contiguous)
-        self.assertArrayEqual(diffs_along_x, np.array([2, 2, 2]).reshape(3, 1))
+        self.assertArrayEqual(diffs_along_x,
+                              np.array([True, True, True]).reshape(3, 1))
         self.assertTrue(not diffs_along_y.any())
 
     def test_2d_discontiguous_along_y(self):
@@ -490,16 +491,16 @@ def test_2d_discontiguous_along_y(self):
         diffs_along_x, diffs_along_y = diffs
         self.assertFalse(contiguous)
         self.assertTrue(not diffs_along_x.any())
-        self.assertArrayEqual(diffs_along_y, np.array([[2, 2, 2]]))
+        self.assertArrayEqual(diffs_along_y, np.array([[True, True, True]]))
 
     def test_2d_discontiguous_along_x_and_y(self):
         coord = AuxCoord(np.array([[1, 5], [3, 5]]),
                          bounds=np.array([[[0, 2, 2, 0], [4, 6, 6, 4]],
                                           [[2, 4, 4, 2], [4, 6, 6, 4]]]))
         contiguous, diffs = coord._discontiguity_in_bounds()
         diffs_along_x, diffs_along_y = diffs
-        exp_x_diffs = np.array([2, 0]).reshape(2, 1)
-        exp_y_diffs = np.array([2, 0]).reshape(1, 2)
+        exp_x_diffs = np.array([True, False]).reshape(2, 1)
+        exp_y_diffs = np.array([True, False]).reshape(1, 2)
         self.assertFalse(contiguous)
         self.assertArrayEqual(diffs_along_x, exp_x_diffs)
         self.assertArrayEqual(diffs_along_y, exp_y_diffs)
@@ -508,10 +509,11 @@ def test_2d_contiguous_along_x_atol(self):
         coord = AuxCoord(self.points_3by3[:, ::2],
                          bounds=self.lon_bounds_3by3[:, ::2, :])
         # Set a high atol that allows small discontiguities.
-        contiguous, diffs = coord._discontiguity_in_bounds(atol=2)
+        contiguous, diffs = coord._discontiguity_in_bounds(atol=5)
         diffs_along_x, diffs_along_y = diffs
         self.assertTrue(contiguous)
-        self.assertArrayEqual(diffs_along_x, np.array([2, 2, 2]).reshape(3, 1))
+        self.assertArrayEqual(diffs_along_x,
+                              np.array([False, False, False]).reshape(3, 1))
         self.assertTrue(not diffs_along_y.any())
 
     def test_2d_one_cell(self):
@@ -571,7 +573,7 @@ def test_2d_discontiguous_mod_360(self):
         contiguous, diffs = coord._discontiguity_in_bounds()
         diffs_along_x, diffs_along_y = diffs
         self.assertFalse(contiguous)
-        self.assertArrayEqual(diffs_along_x, np.array([[170], [170]]))
+        self.assertArrayEqual(diffs_along_x, np.array([[True], [True]]))
         self.assertTrue(not diffs_along_y.any())
 
     def test_2d_contiguous_mod_360_not_longitude(self):
@@ -599,7 +601,7 @@ def test_2d_discontiguous_mod_360_not_longitude(self):
         contiguous, diffs = coord._discontiguity_in_bounds()
         diffs_along_x, diffs_along_y = diffs
         self.assertFalse(contiguous)
-        self.assertArrayEqual(diffs_along_x, np.array([[100], [100]]))
+        self.assertArrayEqual(diffs_along_x, np.array([[True], [True]]))
         self.assertTrue(not diffs_along_y.any())
 
 

lib/iris/tests/unit/plot/test__check_bounds_contiguity_and_mask.py

@@ -33,6 +33,9 @@
 from iris.tests.stock import (sample_2d_latlons,
                               make_bounds_discontiguous_at_point)
 
+from iris.plot import _check_bounds_contiguity_and_mask
+
+
 if tests.MPL_AVAILABLE:
     import iris.plot as iplt
 
@@ -43,20 +46,21 @@ def test_1d_not_checked(self):
         # Test a 1D coordinate, which is not checked as atol is not set.
         coord = DimCoord([1, 3, 5], bounds=[[0, 2], [2, 4], [5, 6]])
         data = np.array([278, 300, 282])
-        iplt._check_bounds_contiguity_and_mask(coord, data)
+        # Make sure contiguity checking doesn't throw an error
+        _check_bounds_contiguity_and_mask(coord, data)
 
     def test_1d_contiguous(self):
-        # Test a 1D coordinate which is contiguous.
+        # Test that a 1D coordinate which is contiguous does not fail.
         coord = DimCoord([1, 3, 5], bounds=[[0, 2], [2, 4], [4, 6]])
         data = np.array([278, 300, 282])
-        iplt._check_bounds_contiguity_and_mask(coord, data, atol=1e-3)
+        _check_bounds_contiguity_and_mask(coord, data, atol=1e-3)
 
     def test_1d_discontigous_masked(self):
         # Test a 1D coordinate which is discontiguous but masked at
         # discontiguities.
         coord = DimCoord([1, 3, 5], bounds=[[0, 2], [2, 4], [5, 6]])
         data = ma.array(np.array([278, 300, 282]), mask=[0, 1, 0])
-        iplt._check_bounds_contiguity_and_mask(coord, data, atol=1e-3)
+        _check_bounds_contiguity_and_mask(coord, data, atol=1e-3)
 
     def test_1d_discontigous_unmasked(self):
         # Test a 1D coordinate which is discontiguous and unmasked at
@@ -66,13 +70,13 @@ def test_1d_discontigous_unmasked(self):
         msg = 'coordinate are not contiguous and data is not masked where ' \
               'the discontiguity occurs'
         with self.assertRaisesRegexp(ValueError, msg):
-            iplt._check_bounds_contiguity_and_mask(coord, data, atol=1e-3)
+            _check_bounds_contiguity_and_mask(coord, data, atol=1e-3)
 
     def test_2d_contiguous(self):
-        # Test a 2D coordinate which is contiguous.
+        # Test that a 2D coordinate which is contiguous does not throw
+        # an error.
         cube = sample_2d_latlons()
-        iplt._check_bounds_contiguity_and_mask(cube.coord('longitude'),
-                                               cube.data)
+        _check_bounds_contiguity_and_mask(cube.coord('longitude'), cube.data)
 
     def test_2d_contiguous_atol(self):
         # Check the atol is passed correctly.
@@ -82,18 +86,16 @@ def test_2d_contiguous_atol(self):
             # Discontiguity returns two objects that are unpacked in
             # `_check_bounds_contiguity_and_mask`.
             discontiguity_check.return_value = [True, None]
-            iplt._check_bounds_contiguity_and_mask(cube.coord('longitude'),
-                                                   cube.data,
-                                                   atol=1e-3)
+            _check_bounds_contiguity_and_mask(cube.coord('longitude'),
+                                              cube.data, atol=1e-3)
         discontiguity_check.assert_called_with(atol=1e-3)
 
     def test_2d_discontigous_masked(self):
-        # Test a 2D coordinate which is discontiguous but masked at
-        # discontiguities.
+        # Test that a 2D coordinate which is discontiguous but masked at
+        # discontiguities doesn't error.
         cube = sample_2d_latlons()
         make_bounds_discontiguous_at_point(cube, 3, 4)
-        iplt._check_bounds_contiguity_and_mask(cube.coord('longitude'),
-                                               cube.data)
+        _check_bounds_contiguity_and_mask(cube.coord('longitude'), cube.data)
 
     def test_2d_discontigous_unmasked(self):
         # Test a 2D coordinate which is discontiguous and unmasked at
@@ -103,8 +105,8 @@ def test_2d_discontigous_unmasked(self):
         msg = 'coordinate are not contiguous'
         cube.data[3, 4] = ma.nomask
         with self.assertRaisesRegexp(ValueError, msg):
-            iplt._check_bounds_contiguity_and_mask(cube.coord('longitude'),
-                                                   cube.data)
+            _check_bounds_contiguity_and_mask(cube.coord('longitude'),
+                                              cube.data)
 
 
 if __name__ == "__main__":