diff --git a/docs/src/whatsnew/latest.rst b/docs/src/whatsnew/latest.rst index 5e8dfef252..7e15674946 100644 --- a/docs/src/whatsnew/latest.rst +++ b/docs/src/whatsnew/latest.rst @@ -72,6 +72,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 index 4245fe55ef..26d2394aad 100644 --- a/lib/iris/coords.py +++ b/lib/iris/coords.py @@ -1936,11 +1936,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() @@ -1949,7 +1950,9 @@ 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: @@ -1957,31 +1960,55 @@ 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 index ff96ecc35e..b7fda5fa40 100644 --- 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 index 72a48437ec..69b6b70c96 100644 --- 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 index e939416e7d..9e043c71bd 100644 --- 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 index d27aa0c722..0bb9a5e365 100644 --- a/lib/iris/util.py +++ b/lib/iris/util.py @@ -1705,14 +1705,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: