Repository: ESMValGroup/ESMValCore

Issue:
Dataset problem: Nemo grid datasets cannot be regridded using iris-esmf
**Describe the dataset issue**
When trying to regrid ocean variables that are in a dataset that use the Nemo grid (CNRM, CMCC, EC-Earth, IPSL), the regridder in iris-esmf  fails due to the coordinates not being contiguous. In those particular datasets, removing the first grid point (which is always masked in ocean variables) seems to be enough in order to make the the coordinates contiguous and have the regridder work properly.

It would be good to have a fix that removes this point, even though it would mean that the data after fixing would have a different shape than the original data. I already worked on the fix, but it would be good to know if anyone has further comments.

**Assign the right people**
If you are already a member of the ESMValTool GitHub project, please assign Valeriu Predoi (valeriupredoi) and
Klaus Zimmermann (zklaus) to the issue. They will then check the issue raised and propagate it further to the
data model developers.

**Please attach**
  - The recipe that you are trying to run, you can find a copy in the `run` directory in the output directory
  - The `main_log_debug.txt` file, this can also be found in the `run` directory in the output directory
  - Attach a text file with the output of `ncdump -h /path/to/file.nc`

Cheers :beer:


Discussion:
sloosvel: > Would it make sense to add some documentation for this feature?

Good point but I cannot make it happen for today, it's a holiday here and I am having trouble connecting with my laptop. In any case `use_src_mask` can always be called from the recipe.
valeriupredoi: and you're also having lunch :rofl: :pizza: 
bouweandela: Thanks for adding the documentation. Would it make sense to add it to the `esmvalcore.preprocessor.regrid` as well, so it shows up in the API docs?

Patch:
diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst
index 44eb9d0dd..1c9140868 100644
--- a/doc/recipe/preprocessor.rst
+++ b/doc/recipe/preprocessor.rst
@@ -1015,6 +1015,13 @@ scheme available in :doc:`iris-esmf-regrid:index`:
           reference: esmf_regrid.schemes:regrid_rectilinear_to_rectilinear
           mdtol: 0.7
 
+Since version 0.6 of :doc:`iris-esmf-regrid:index`, the regridder is able
+to ignore discontinuities in the source grids if the data defined in the
+discontiguous points is masked.
+The :func:`~esmvalcore.preprocessor.regrid` preprocessor automatically detects
+if discontinuities are present, and configures the regridding scheme in order to take
+into account the mask of the source grid to ignore them.
+
 .. _ensemble statistics:
 
 Ensemble statistics
diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py
index fe94744ff..cb5205559 100644
--- a/esmvalcore/preprocessor/_regrid.py
+++ b/esmvalcore/preprocessor/_regrid.py
@@ -384,6 +384,21 @@ def extract_location(cube, location, scheme):
                          scheme)
 
 
+def _check_grid_discontiguities(cube, scheme):
+    """Check if there are grid discontiguities and set use_src_mask to True."""
+    scheme = dict(scheme)
+    try:
+        discontiguities = iris.util.find_discontiguities(cube)
+    except NotImplementedError:
+        pass
+    else:
+        if discontiguities.any():
+            scheme['use_src_mask'] = True
+            logger.debug('Grid discontinuities were found in the source grid. '
+                         'Setting scheme argument `use_src_mask` to True.')
+    return scheme
+
+
 def extract_point(cube, latitude, longitude, scheme):
     """Extract a point, with interpolation.
 
@@ -561,6 +576,13 @@ def regrid(cube, target_grid, scheme, lat_offset=True, lon_offset=True):
             target: 1x1
             scheme:
               reference: esmf_regrid.schemes:ESMFAreaWeighted
+
+    Since version 0.6 of :doc:`iris-esmf-regrid:index`, the regridder is able
+    to ignore discontinuities in the source grids if the data defined in the
+    discontiguous points is masked.
+    The preprocessor automatically detects if discontinuities are present,
+    and configures the regridding scheme in order to take into account
+    the mask of the source grid to ignore them.
     """
     if is_dataset(target_grid):
         target_grid = target_grid.copy()
@@ -615,6 +637,8 @@ def regrid(cube, target_grid, scheme, lat_offset=True, lon_offset=True):
             scheme['src_cube'] = cube
         if 'grid_cube' in scheme_args:
             scheme['grid_cube'] = target_grid
+        if 'use_src_mask' in scheme_args:
+            scheme = _check_grid_discontiguities(cube, scheme)
 
         loaded_scheme = obj(**scheme)
     else:


Test patch:
diff --git a/tests/unit/preprocessor/_regrid/test_regrid.py b/tests/unit/preprocessor/_regrid/test_regrid.py
index 9a7383551..a964eaf64 100644
--- a/tests/unit/preprocessor/_regrid/test_regrid.py
+++ b/tests/unit/preprocessor/_regrid/test_regrid.py
@@ -1,6 +1,7 @@
 """Unit tests for the :func:`esmvalcore.preprocessor.regrid.regrid`
 function."""
 
+import logging
 import unittest
 from unittest import mock
 
@@ -15,12 +16,14 @@ from esmvalcore.preprocessor import regrid
 from esmvalcore.preprocessor._regrid import (
     _CACHE,
     HORIZONTAL_SCHEMES,
+    _check_grid_discontiguities,
     _horizontal_grid_is_close,
     _rechunk,
 )
 
 
 class Test(tests.Test):
+
     def _check(self, tgt_grid, scheme, spec=False):
         expected_scheme = HORIZONTAL_SCHEMES[scheme]
 
@@ -70,8 +73,8 @@ class Test(tests.Test):
         )
         self.src_cube.ndim = 1
         self.tgt_grid_coord = mock.Mock()
-        self.tgt_grid = mock.Mock(
-            spec=iris.cube.Cube, coord=self.tgt_grid_coord)
+        self.tgt_grid = mock.Mock(spec=iris.cube.Cube,
+                                  coord=self.tgt_grid_coord)
         self.regrid_schemes = [
             'linear', 'linear_extrapolate', 'nearest', 'area_weighted',
             'unstructured_nearest'
@@ -111,8 +114,8 @@ class Test(tests.Test):
             regrid(self.src_cube, self.src_cube, 'wibble')
 
     def test_horizontal_schemes(self):
-        self.assertEqual(
-            set(HORIZONTAL_SCHEMES.keys()), set(self.regrid_schemes))
+        self.assertEqual(set(HORIZONTAL_SCHEMES.keys()),
+                         set(self.regrid_schemes))
 
     def test_regrid__horizontal_schemes(self):
         for scheme in self.regrid_schemes:
@@ -154,12 +157,12 @@ class Test(tests.Test):
     third_party_regridder = mock.Mock()
 
     def test_regrid_generic_third_party(self):
-        regrid(self.src_cube, '1x1',
-               {"reference":
-                "tests.unit.preprocessor._regrid.test_regrid:"
+        regrid(
+            self.src_cube, '1x1', {
+                "reference": "tests.unit.preprocessor._regrid.test_regrid:"
                 "Test.third_party_regridder",
                 "method": "good",
-                })
+            })
         self.third_party_regridder.assert_called_once_with(method="good")
 
 
@@ -266,6 +269,62 @@ def test_regrid_is_skipped_if_grids_are_the_same():
     assert expected_different_cube is not cube
 
 
+def test_no_discontiguities_in_coords():
+    """Test that no mask is used if there are no discontinuities in coords."""
+    cube = _make_cube(lat=LAT_SPEC1, lon=LON_SPEC1)
+    scheme = {}
+    scheme = _check_grid_discontiguities(cube, scheme)
+    assert scheme == {}
+
+
+def test_use_mask_if_discontiguities_in_coords(caplog):
+    """Test use_src_mask is added to the scheme."""
+    lat_bounds = np.array(
+        [[[-43.48076211, -34.01923789, -22.00961894, -31.47114317],
+          [-34.01923789, -10.0, 2.00961894, -22.00961894],
+          [-10.0, -0.53847577, 11.47114317, 2.00961894]],
+         [[-31.47114317, -22.00961894, -10.0, -19.46152423],
+          [-22.00961894, 2.00961894, 14.01923789, -10.0],
+          [2.00961894, 11.47114317, 23.48076211, 14.01923789]]])
+    lat_coord = iris.coords.AuxCoord(
+        [[-40.0, -20.0, 0.0], [-20.0, 0.0, 20.0]],
+        var_name='lat',
+        standard_name='latitude',
+        units='degrees_north',
+        bounds=lat_bounds,
+    )
+    lon_bounds = np.array([[[140.625, 99.375, 99.375, 140.625],
+                            [99.375, 140.625, 140.625, 99.375],
+                            [140.625, 99.375, 99.375, 140.625]],
+                           [[140., 99.375, 99.375, 140.],
+                            [99.375, 140.625, 140.625, 99.375],
+                            [140., 99.375, 99.375, 140.]]])
+    lon_coord = iris.coords.AuxCoord(
+        [[100.0, 140.0, 180.0], [80.0, 100.0, 120.0]],
+        var_name='lon',
+        standard_name='longitude',
+        units='degrees_east',
+        bounds=lon_bounds,
+    )
+    data = np.ma.array(
+        [[-40.0, -20.0, 0.0], [-20.0, 0.0, 20.0]],
+        mask=[[True, False, True], [False, True, False]],
+    )
+    cube = iris.cube.Cube(
+        data,
+        aux_coords_and_dims=[(lat_coord, (0, 1)), (lon_coord, (0, 1))],
+    )
+
+    scheme = {}
+    with caplog.at_level(logging.DEBUG):
+        scheme = _check_grid_discontiguities(cube, scheme)
+    assert scheme == {'use_src_mask': True}
+
+    msg = ('Grid discontinuities were found in the source grid. '
+           'Setting scheme argument `use_src_mask` to True.')
+    assert msg in caplog.text
+
+
 def test_rechunk_on_increased_grid():
     """Test that an increase in grid size rechunks."""
     with dask.config.set({'array.chunk-size': '128 M'}):