Fixes to _discontiguity_in_bounds (attempt 2) (#4975)

stephenworsley · web-flow · commit 6b4ba73cd306 · 2023-04-12T16:06:36.000+01:00
diff --git a/docs/src/whatsnew/latest.rst b/docs/src/whatsnew/latest.rst
@@ -75,6 +75,14 @@ This document explains the changes made to Iris for this release
    :class:`~iris.coords.CellMethod` are printed to be more CF compliant.
    (:pull:`5224`)
 
+#. `@stephenworsley`_ fixed the way discontiguities were discovered for 2D coords.
+   Previously, the only bounds being compared were the bottom right bound in one
+   cell with the bottom left bound in the cell to its right, and the top left bound
+   in a cell with the bottom left bound in the cell above it. Now all bounds are
+   compared with all adjacent bounds from neighbouring cells. This affects
+   :meth:`~iris.coords.Coord.is_contiguous` and :func:`iris.util.find_discontiguities`
+   where additional discontiguities may be detected which previously were not.
+
 
 💣 Incompatible Changes
 =======================
diff --git a/lib/iris/coords.py b/lib/iris/coords.py
@@ -1932,11 +1932,12 @@ def _discontiguity_in_bounds(self, rtol=1e-5, atol=1e-8):
         * contiguous: (boolean)
             True if there are no discontiguities.
         * diffs: (array or tuple of arrays)
-            The diffs along the bounds of the coordinate. If self is a 2D
-            coord of shape (Y, X), a tuple of arrays is returned, where the
-            first is an array of differences along the x-axis, of the shape
-            (Y, X-1) and the second is an array of differences along the
-            y-axis, of the shape (Y-1, X).
+            A boolean array or tuple of boolean arrays which are true where
+            there are discontiguities between neighbouring bounds. If self is
+            a 2D coord of shape (Y, X), a pair of arrays is returned, where
+            the first is an array of differences along the x-axis, of the
+            shape (Y, X-1) and the second is an array of differences along
+            the y-axis, of the shape (Y-1, X).
 
         """
         self._sanity_check_bounds()
@@ -1945,39 +1946,65 @@ def _discontiguity_in_bounds(self, rtol=1e-5, atol=1e-8):
             contiguous = np.allclose(
                 self.bounds[1:, 0], self.bounds[:-1, 1], rtol=rtol, atol=atol
             )
-            diffs = np.abs(self.bounds[:-1, 1] - self.bounds[1:, 0])
+            diffs = ~np.isclose(
+                self.bounds[1:, 0], self.bounds[:-1, 1], rtol=rtol, atol=atol
+            )
 
         elif self.ndim == 2:
 
             def mod360_adjust(compare_axis):
                 bounds = self.bounds.copy()
 
                 if compare_axis == "x":
-                    upper_bounds = bounds[:, :-1, 1]
-                    lower_bounds = bounds[:, 1:, 0]
+                    # Extract the pairs of upper bounds and lower bounds which
+                    # connect along the "x" axis. These connect along indices
+                    # as shown by the following diagram:
+                    #
+                    # 3---2 + 3---2
+                    # |   |   |   |
+                    # 0---1 + 0---1
+                    upper_bounds = np.stack(
+                        (bounds[:, :-1, 1], bounds[:, :-1, 2])
+                    )
+                    lower_bounds = np.stack(
+                        (bounds[:, 1:, 0], bounds[:, 1:, 3])
+                    )
                 elif compare_axis == "y":
-                    upper_bounds = bounds[:-1, :, 3]
-                    lower_bounds = bounds[1:, :, 0]
+                    # Extract the pairs of upper bounds and lower bounds which
+                    # connect along the "y" axis. These connect along indices
+                    # as shown by the following diagram:
+                    #
+                    # 3---2
+                    # |   |
+                    # 0---1
+                    # +   +
+                    # 3---2
+                    # |   |
+                    # 0---1
+                    upper_bounds = np.stack(
+                        (bounds[:-1, :, 3], bounds[:-1, :, 2])
+                    )
+                    lower_bounds = np.stack(
+                        (bounds[1:, :, 0], bounds[1:, :, 1])
+                    )
 
                 if self.name() in ["longitude", "grid_longitude"]:
                     # If longitude, adjust for longitude wrapping
                     diffs = upper_bounds - lower_bounds
-                    index = diffs > 180
+                    index = np.abs(diffs) > 180
                     if index.any():
                         sign = np.sign(diffs)
                         modification = (index.astype(int) * 360) * sign
                         upper_bounds -= modification
 
-                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
+                diffs_along_bounds = ~np.isclose(
+                    upper_bounds, lower_bounds, rtol=rtol, atol=atol
                 )
-
-                points_close_enough = diffs_along_axis <= (
-                    atol + rtol * cell_size
+                diffs_along_axis = np.logical_or(
+                    diffs_along_bounds[0], diffs_along_bounds[1]
                 )
-                contiguous_along_axis = np.all(points_close_enough)
+
+                contiguous_along_axis = ~np.any(diffs_along_axis)
                 return diffs_along_axis, contiguous_along_axis
 
             diffs_along_x, match_cell_x1 = mod360_adjust(compare_axis="x")
diff --git a/lib/iris/tests/stock/_stock_2d_latlons.py b/lib/iris/tests/stock/_stock_2d_latlons.py
@@ -296,7 +296,9 @@ def sample_cube(xargs, yargs):
     return cube
 
 
-def make_bounds_discontiguous_at_point(cube, at_iy, at_ix, in_y=False):
+def make_bounds_discontiguous_at_point(
+    cube, at_iy, at_ix, in_y=False, upper=True
+):
     """
     Meddle with the XY grid bounds of a 2D cube to make the grid discontiguous.
 
@@ -325,16 +327,22 @@ def adjust_coord(coord):
         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
+            new_bds_b = 0.5 * (bds_bl + bds_br)
+            new_bds_t = 0.5 * (bds_tl + bds_tr)
+            if upper:
+                bds_br, bds_tr = new_bds_b, new_bds_t
+            else:
+                bds_bl, bds_tl = new_bds_b, new_bds_t
         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
+            new_bds_l = 0.5 * (bds_bl + bds_tl)
+            new_bds_r = 0.5 * (bds_br + bds_tr)
+            if upper:
+                bds_tl, bds_tr = new_bds_l, new_bds_r
+            else:
+                bds_bl, bds_br = new_bds_l, new_bds_r
 
         # Write in the new bounds (all 4 corners).
         bds[at_iy, at_ix] = [bds_bl, bds_br, bds_tr, bds_tl]
@@ -355,7 +363,16 @@ def adjust_coord(coord):
             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
+    masked_data = ma.masked_array(cube.data)
+
+    # Mask all points which would be found discontiguous.
+    # Note that find_discontiguities finds all instances where a cell is
+    # discontiguous with a neighbouring cell to its *right* or *above*
+    # that cell.
+    masked_data[at_iy, at_ix] = ma.masked
+    if in_y or not upper:
+        masked_data[at_iy, at_ix - 1] = ma.masked
+    if not in_y or not upper:
+        masked_data[at_iy - 1, at_ix] = ma.masked
 
     cube.data = masked_data
diff --git a/lib/iris/tests/unit/coords/test_Coord.py b/lib/iris/tests/unit/coords/test_Coord.py
@@ -708,7 +708,7 @@ def test_1d_discontiguous(self):
         coord = DimCoord([10, 20, 40], bounds=[[5, 15], [15, 25], [35, 45]])
         contiguous, diffs = coord._discontiguity_in_bounds()
         self.assertFalse(contiguous)
-        self.assertArrayEqual(diffs, np.array([0, 10]))
+        self.assertArrayEqual(diffs, np.array([False, True]))
 
     def test_1d_one_cell(self):
         # Test a 1D coord with a single cell.
diff --git a/lib/iris/tests/unit/util/test_find_discontiguities.py b/lib/iris/tests/unit/util/test_find_discontiguities.py
@@ -29,26 +29,55 @@ def setUp(self):
         # Set up a 2d lat-lon cube with 2d coordinates that have been
         # transformed so they are not in a regular lat-lon grid.
         # Then generate a discontiguity at a single lat-lon point.
-        self.testcube_discontig = full2d_global()
-        make_bounds_discontiguous_at_point(self.testcube_discontig, 3, 3)
-        # Repeat that for a discontiguity in the grid 'Y' direction.
-        self.testcube_discontig_along_y = full2d_global()
+        # Discontiguities will be caused at the rightmost bounds.
+        self.testcube_discontig_right = full2d_global()
+        make_bounds_discontiguous_at_point(self.testcube_discontig_right, 3, 3)
+
+        # Repeat for a discontiguity on the leftmost bounds.
+        self.testcube_discontig_left = full2d_global()
+        make_bounds_discontiguous_at_point(
+            self.testcube_discontig_left, 2, 4, upper=False
+        )
+        # Repeat for a discontiguity on the topmost bounds.
+        self.testcube_discontig_top = full2d_global()
         make_bounds_discontiguous_at_point(
-            self.testcube_discontig_along_y, 2, 4, in_y=True
+            self.testcube_discontig_top, 2, 4, in_y=True
         )
 
-    def test_find_discontiguities(self):
+        # Repeat for a discontiguity on the botommost bounds.
+        self.testcube_discontig_along_bottom = full2d_global()
+        make_bounds_discontiguous_at_point(
+            self.testcube_discontig_along_bottom, 2, 4, in_y=True, upper=False
+        )
+
+    def test_find_discontiguities_right(self):
+        # Check that the mask we generate when making the discontiguity
+        # matches that generated by find_discontiguities
+        cube = self.testcube_discontig_right
+        expected = cube.data.mask
+        returned = find_discontiguities(cube)
+        self.assertTrue(np.all(expected == returned))
+
+    def test_find_discontiguities_left(self):
+        # Check that the mask we generate when making the discontiguity
+        # matches that generated by find_discontiguities
+        cube = self.testcube_discontig_left
+        expected = cube.data.mask
+        returned = find_discontiguities(cube)
+        self.assertTrue(np.all(expected == returned))
+
+    def test_find_discontiguities_top(self):
         # Check that the mask we generate when making the discontiguity
         # matches that generated by find_discontiguities
-        cube = self.testcube_discontig
+        cube = self.testcube_discontig_top
         expected = cube.data.mask
         returned = find_discontiguities(cube)
         self.assertTrue(np.all(expected == returned))
 
-    def test_find_discontiguities_in_y(self):
+    def test_find_discontiguities_bottom(self):
         # Check that the mask we generate when making the discontiguity
         # matches that generated by find_discontiguities
-        cube = self.testcube_discontig_along_y
+        cube = self.testcube_discontig_along_bottom
         expected = cube.data.mask
         returned = find_discontiguities(cube)
         self.assertTrue(np.all(expected == returned))
@@ -61,7 +90,7 @@ def test_find_discontiguities_1d_coord(self):
             find_discontiguities(cube)
 
     def test_find_discontiguities_with_atol(self):
-        cube = self.testcube_discontig
+        cube = self.testcube_discontig_right
         # Choose a very large absolute tolerance which will result in fine
         # discontiguities being disregarded
         atol = 100
@@ -72,7 +101,7 @@ def test_find_discontiguities_with_atol(self):
         self.assertTrue(np.all(expected == returned))
 
     def test_find_discontiguities_with_rtol(self):
-        cube = self.testcube_discontig
+        cube = self.testcube_discontig_right
         # Choose a very large relative tolerance which will result in fine
         # discontiguities being disregarded
         rtol = 1000
diff --git a/lib/iris/util.py b/lib/iris/util.py
@@ -1820,14 +1820,15 @@ def _meshgrid(*xi, **kwargs):
 
 def find_discontiguities(cube, rel_tol=1e-5, abs_tol=1e-8):
     """
-    Searches coord for discontiguities in the bounds array, returned as a
-    boolean array (True where discontiguities are present).
+    Searches the 'x' and 'y' coord on the cube for discontiguities in the
+    bounds array, returned as a boolean array (True for all cells which are
+    discontiguous with the cell immediately above them or to their right).
 
     Args:
 
     * cube (`iris.cube.Cube`):
         The cube to be checked for discontinuities in its 'x' and 'y'
-        coordinates.
+        coordinates. These coordinates must be 2D.
 
     Kwargs:
 
