Allow boundless coordinates for bilinear and nearest regridders (#276)

* allow boundless coordinates for bilinear and nearest regridders

* add documentation

* fix test

* adapt docs to iris changes

* fix docs

CHANGELOG.md

@@ -51,6 +51,13 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
   Added support for Mesh to Mesh regridding.
   [@HGWright](https://github.com/HGWright)
 
+### Added
+
+- [PR#276](https://github.com/SciTools-incubator/iris-esmf-regrid/pull/276)
+  Allow regridding for grids defined on coordinates without bounds for
+  nearest neighbour and bilinear methods.
+  [@stephenworsley](https://github.com/stephenworsley)
+
 ### Fixed
 
 - [PR#301](https://github.com/SciTools-incubator/iris-esmf-regrid/pull/301)

docs/src/userguide/examples.rst

@@ -31,7 +31,7 @@ Saving and Loading a Regridder
 A regridder can be set up for reuse, this saves time performing the
 computationally expensive initialisation process::
 
-    from esmf_regrid.experimental.unstructured_scheme import ESMFAreaWeighted
+    from esmf_regrid.schemes import ESMFAreaWeighted
 
     # Initialise the regridder with a source mesh and target grid.
     regridder = ESMFAreaWeighted().regridder(source_mesh_cube, target_grid_cube)

esmf_regrid/_esmf_sdo.py

@@ -162,15 +162,15 @@ def __init__(
         self.lons = lons
         self.lats = lats
         londims = len(self.lons.shape)
-        if len(lonbounds.shape) != londims:
+        if lonbounds is not None and len(lonbounds.shape) != londims:
             msg = (
                 f"The dimensionality of longitude bounds "
                 f"({len(lonbounds.shape)}) is incompatible with the "
                 f"dimensionality of the longitude ({londims})."
             )
             raise ValueError(msg)
         latdims = len(self.lats.shape)
-        if len(latbounds.shape) != latdims:
+        if latbounds is not None and len(latbounds.shape) != latdims:
             msg = (
                 f"The dimensionality of latitude bounds "
                 f"({len(latbounds.shape)}) is incompatible with the "
@@ -219,35 +219,44 @@ def _as_esmf_info(self):
         londims = len(self.lons.shape)
 
         if londims == 1:
-            if self.circular:
-                adjustedlonbounds = self._refined_lonbounds[:-1]
-            else:
-                adjustedlonbounds = self._refined_lonbounds
+            if self._refined_lonbounds is not None:
+                if self.circular:
+                    adjustedlonbounds = self._refined_lonbounds[:-1]
+                else:
+                    adjustedlonbounds = self._refined_lonbounds
+                cornerlons, cornerlats = np.meshgrid(
+                    adjustedlonbounds, self._refined_latbounds
+                )
             centerlons, centerlats = np.meshgrid(self.lons, self.lats)
-            cornerlons, cornerlats = np.meshgrid(
-                adjustedlonbounds, self._refined_latbounds
-            )
         elif londims == 2:
             if self.circular:
                 slice = np.s_[:, :-1]
             else:
                 slice = np.s_[:]
-            centerlons = self.lons[slice]
-            centerlats = self.lats[slice]
-            cornerlons = self._refined_lonbounds[slice]
-            cornerlats = self._refined_latbounds[slice]
+            centerlons = self.lons[:]
+            centerlats = self.lats[:]
+            if self._refined_lonbounds is not None:
+                cornerlons = self._refined_lonbounds[slice]
+                cornerlats = self._refined_latbounds[slice]
 
         truecenters = ccrs.Geodetic().transform_points(self.crs, centerlons, centerlats)
-        truecorners = ccrs.Geodetic().transform_points(self.crs, cornerlons, cornerlats)
+        if self._refined_lonbounds is not None:
+            truecorners = ccrs.Geodetic().transform_points(
+                self.crs, cornerlons, cornerlats
+            )
 
         # The following note in xESMF suggests that the arrays passed to ESMPy ought to
         # be fortran ordered:
         # https://xesmf.readthedocs.io/en/latest/internal_api.html#xesmf.backend.warn_f_contiguous
         # It is yet to be determined what effect this has on performance.
         truecenterlons = np.asfortranarray(truecenters[..., 0])
         truecenterlats = np.asfortranarray(truecenters[..., 1])
-        truecornerlons = np.asfortranarray(truecorners[..., 0])
-        truecornerlats = np.asfortranarray(truecorners[..., 1])
+        if self._refined_lonbounds is not None:
+            truecornerlons = np.asfortranarray(truecorners[..., 0])
+            truecornerlats = np.asfortranarray(truecorners[..., 1])
+        else:
+            truecornerlons = None
+            truecornerlats = None
 
         info = (
             shape,
@@ -283,12 +292,6 @@ def _make_esmf_sdo(self):
         else:
             grid = esmpy.Grid(shape, pole_kind=[1, 1])
 
-        grid.add_coords(staggerloc=esmpy.StaggerLoc.CORNER)
-        grid_corner_x = grid.get_coords(0, staggerloc=esmpy.StaggerLoc.CORNER)
-        grid_corner_x[:] = truecornerlons
-        grid_corner_y = grid.get_coords(1, staggerloc=esmpy.StaggerLoc.CORNER)
-        grid_corner_y[:] = truecornerlats
-
         # Grid center points are added here, this is not necessary for
         # conservative area weighted regridding
         if self.center:
@@ -297,6 +300,12 @@ def _make_esmf_sdo(self):
             grid_center_x[:] = truecenterlons
             grid_center_y = grid.get_coords(1, staggerloc=esmpy.StaggerLoc.CENTER)
             grid_center_y[:] = truecenterlats
+        else:
+            grid.add_coords(staggerloc=esmpy.StaggerLoc.CORNER)
+            grid_corner_x = grid.get_coords(0, staggerloc=esmpy.StaggerLoc.CORNER)
+            grid_corner_x[:] = truecornerlons
+            grid_corner_y = grid.get_coords(1, staggerloc=esmpy.StaggerLoc.CORNER)
+            grid_corner_y[:] = truecornerlats
 
         def add_get_item(grid, **kwargs):
             grid.add_item(**kwargs)

esmf_regrid/experimental/unstructured_scheme.py

@@ -294,9 +294,9 @@ def __init__(
         ----------
         src : :class:`iris.cube.Cube`
             The rectilinear :class:`~iris.cube.Cube` cube providing the source grid.
-        tgt : :class:`iris.cube.Cube` or :class:`iris.experimental.ugrid.MeshXY`
+        tgt : :class:`iris.cube.Cube` or :class:`iris.mesh.MeshXY`
             The unstructured :class:`~iris.cube.Cube`or
-            :class:`~iris.experimental.ugrid.MeshXY` providing the target mesh.
+            :class:`~iris.mesh.MeshXY` providing the target mesh.
         mdtol : float, optional
             Tolerance of missing data. The value returned in each element of
             the returned array will be masked if the fraction of masked data

esmf_regrid/schemes.py

@@ -179,31 +179,40 @@ def _cube_to_GridInfo(cube, center=False, resolution=None, mask=None):
     londim, latdim = len(lon.shape), len(lat.shape)
     assert londim == latdim
     assert londim in (1, 2)
-    if londim == 1:
-        # Ensure coords come from a proper grid.
-        assert isinstance(lon, iris.coords.DimCoord)
-        assert isinstance(lat, iris.coords.DimCoord)
-        circular = lon.circular
-        lon_bound_array = lon.contiguous_bounds()
-        lat_bound_array = lat.contiguous_bounds()
-        # TODO: perform checks on lat/lon.
-    elif londim == 2:
-        assert cube.coord_dims(lon) == cube.coord_dims(lat)
-        if not np.any(mask):
-            assert lon.is_contiguous()
-            assert lat.is_contiguous()
+    if not center:
+        if londim == 1:
+            # Ensure coords come from a proper grid.
+            assert isinstance(lon, iris.coords.DimCoord)
+            assert isinstance(lat, iris.coords.DimCoord)
+            circular = lon.circular
             lon_bound_array = lon.contiguous_bounds()
             lat_bound_array = lat.contiguous_bounds()
+            # TODO: perform checks on lat/lon.
+        elif londim == 2:
+            assert cube.coord_dims(lon) == cube.coord_dims(lat)
+            if not np.any(mask):
+                assert lon.is_contiguous()
+                assert lat.is_contiguous()
+                lon_bound_array = lon.contiguous_bounds()
+                lat_bound_array = lat.contiguous_bounds()
+            else:
+                lon_bound_array = _contiguous_masked(lon.bounds, mask)
+                lat_bound_array = _contiguous_masked(lat.bounds, mask)
+            # 2D coords must be AuxCoords, which do not have a circular attribute.
+            circular = False
+        if crs is None:
+            lon_bound_array = lon.units.convert(lon_bound_array, Unit("degrees"))
+            lat_bound_array = lat.units.convert(lat_bound_array, Unit("degrees"))
+    else:
+        lon_bound_array = None
+        lat_bound_array = None
+        if londim == 1:
+            circular = lon.circular
         else:
-            lon_bound_array = _contiguous_masked(lon.bounds, mask)
-            lat_bound_array = _contiguous_masked(lat.bounds, mask)
-        # 2D coords must be AuxCoords, which do not have a circular attribute.
-        circular = False
+            circular = False
     lon_points = lon.points
     lat_points = lat.points
     if crs is None:
-        lon_bound_array = lon.units.convert(lon_bound_array, Unit("degrees"))
-        lat_bound_array = lat.units.convert(lat_bound_array, Unit("degrees"))
         lon_points = lon.units.convert(lon_points, Unit("degrees"))
         lat_points = lon.units.convert(lat_points, Unit("degrees"))
     if resolution is None:
@@ -414,7 +423,7 @@ def _create_cube(data, src_cube, src_dims, tgt_coords, num_tgt_dims):
         The dimensions of the X and Y coordinate within the source Cube.
     tgt_coords : tuple of :class:`iris.coords.Coord`\\ 's
         Either two 1D :class:`iris.coords.DimCoord`\\ 's, two 1D
-        :class:`iris.experimental.ugrid.DimCoord`\\ 's or two 2D
+        :class:`iris.mesh.MeshCoord`\\ 's or two 2D
         :class:`iris.coords.AuxCoord`\\ 's representing the new grid's
         X and Y coordinates.
     num_tgt_dims : int
@@ -1005,9 +1014,13 @@ def regridder(
         ----------
         src_grid : :class:`iris.cube.Cube`
             The :class:`~iris.cube.Cube` defining the source.
-        tgt_grid : :class:`iris.cube.Cube` or :class:`iris.experimental.ugrid.MeshXY`
+            If this cube has a grid defined by latitude/longitude coordinates, those
+            coordinates must have bounds.
+        tgt_grid : :class:`iris.cube.Cube` or :class:`iris.mesh.MeshXY`
             The unstructured :class:`~iris.cube.Cube`or
-            :class:`~iris.experimental.ugrid.MeshXY` defining the target.
+            :class:`~iris.mesh.MeshXY` defining the target.
+            If this cube has a grid defined by latitude/longitude coordinates, those
+            coordinates must have bounds.
         src_resolution, tgt_resolution : int, optional
             If present, represents the amount of latitude slices per source/target cell
             given to ESMF for calculation. If resolution is set, ``src`` and ``tgt``
@@ -1117,9 +1130,9 @@ def regridder(
         ----------
         src_grid : :class:`iris.cube.Cube`
             The :class:`~iris.cube.Cube` defining the source.
-        tgt_grid : :class:`iris.cube.Cube` or :class:`iris.experimental.ugrid.MeshXY`
+        tgt_grid : :class:`iris.cube.Cube` or :class:`iris.mesh.MeshXY`
             The unstructured :class:`~iris.cube.Cube`or
-            :class:`~iris.experimental.ugrid.MeshXY` defining the target.
+            :class:`~iris.mesh.MeshXY` defining the target.
         use_src_mask : :obj:`~numpy.typing.ArrayLike` or bool, optional
             Array describing which elements :mod:`esmpy` will ignore on the src_grid.
             If True, the mask will be derived from src_grid.
@@ -1225,9 +1238,9 @@ def regridder(
         ----------
         src_grid : :class:`iris.cube.Cube`
             The :class:`~iris.cube.Cube` defining the source.
-        tgt_grid : :class:`iris.cube.Cube` or :class:`iris.experimental.ugrid.MeshXY`
+        tgt_grid : :class:`iris.cube.Cube` or :class:`iris.mesh.MeshXY`
             The unstructured :class:`~iris.cube.Cube`or
-            :class:`~iris.experimental.ugrid.MeshXY` defining the target.
+            :class:`~iris.mesh.MeshXY` defining the target.
         use_src_mask : :obj:`~numpy.typing.ArrayLike` or bool, optional
             Array describing which elements :mod:`esmpy` will ignore on the src_grid.
             If True, the mask will be derived from src_grid.
@@ -1288,7 +1301,7 @@ def __init__(
         ----------
         src : :class:`iris.cube.Cube`
             The rectilinear :class:`~iris.cube.Cube` providing the source grid.
-        tgt : :class:`iris.cube.Cube` or :class:`iris.experimental.ugrid.MeshXY`
+        tgt : :class:`iris.cube.Cube` or :class:`iris.mesh.MeshXY`
             The rectilinear :class:`~iris.cube.Cube` providing the target grid.
         method : :class:`Constants.Method`
             The method to be used to calculate weights.
@@ -1465,9 +1478,13 @@ def __init__(
         ----------
         src : :class:`iris.cube.Cube`
             The rectilinear :class:`~iris.cube.Cube` providing the source.
-        tgt : :class:`iris.cube.Cube` or :class:`iris.experimental.ugrid.MeshXY`
+            If this cube has a grid defined by latitude/longitude coordinates, those
+            coordinates must have bounds.
+        tgt : :class:`iris.cube.Cube` or :class:`iris.mesh.MeshXY`
             The unstructured :class:`~iris.cube.Cube`or
-            :class:`~iris.experimental.ugrid.MeshXY` defining the target.
+            :class:`~iris.mesh.MeshXY` defining the target.
+            If this cube has a grid defined by latitude/longitude coordinates, those
+            coordinates must have bounds.
         mdtol : float, default=0
             Tolerance of missing data. The value returned in each element of
             the returned array will be masked if the fraction of masked data
@@ -1543,9 +1560,9 @@ def __init__(
         ----------
         src : :class:`iris.cube.Cube`
             The rectilinear :class:`~iris.cube.Cube` providing the source.
-        tgt : :class:`iris.cube.Cube` or :class:`iris.experimental.ugrid.Mesh`
+        tgt : :class:`iris.cube.Cube` or :class:`iris.mesh.MeshXY`
             The unstructured :class:`~iris.cube.Cube`or
-            :class:`~iris.experimental.ugrid.MeshXY` defining the target.
+            :class:`~iris.mesh.MeshXY` defining the target.
         mdtol : float, default=0
             Tolerance of missing data. The value returned in each element of
             the returned array will be masked if the fraction of masked data
@@ -1602,9 +1619,9 @@ def __init__(
         ----------
         src : :class:`iris.cube.Cube`
             The rectilinear :class:`~iris.cube.Cube` providing the source.
-        tgt : :class:`iris.cube.Cube` or :class:`iris.experimental.ugrid.MeshXY`
+        tgt : :class:`iris.cube.Cube` or :class:`iris.mesh.MeshXY`
             The unstructured :class:`~iris.cube.Cube`or
-            :class:`~iris.experimental.ugrid.MeshXY` defining the target.
+            :class:`~iris.mesh.MeshXY` defining the target.
         precomputed_weights : :class:`scipy.sparse.spmatrix`, optional
             If ``None``, :mod:`esmpy` will be used to
             calculate regridding weights. Otherwise, :mod:`esmpy` will be bypassed

esmf_regrid/tests/unit/schemes/test_ESMFBilinearRegridder.py

@@ -88,7 +88,8 @@ def test_invalid_mdtol():
         _ = ESMFBilinearRegridder(src, tgt, mdtol=-1)
 
 
-def test_curvilinear_equivalence():
+@pytest.mark.parametrize("with_bounds", (False, True))
+def test_curvilinear_equivalence(with_bounds):
     """
     Test initialisation of :class:`esmf_regrid.schemes.ESMFBilinearRegridder`.
 
@@ -105,6 +106,11 @@ def test_curvilinear_equivalence():
     grid_tgt = _grid_cube(n_lons_tgt, n_lats_tgt, lon_bounds, lat_bounds, circular=True)
     curv_src = _curvilinear_cube(n_lons_src, n_lats_src, lon_bounds, lat_bounds)
     curv_tgt = _curvilinear_cube(n_lons_tgt, n_lats_tgt, lon_bounds, lat_bounds)
+    if not with_bounds:
+        curv_src.coord("latitude").bounds = None
+        curv_src.coord("longitude").bounds = None
+        curv_tgt.coord("latitude").bounds = None
+        curv_tgt.coord("longitude").bounds = None
 
     grid_to_grid = ESMFBilinearRegridder(grid_src, grid_tgt)
     grid_to_curv = ESMFBilinearRegridder(grid_src, curv_tgt)
@@ -153,10 +159,10 @@ def test_curvilinear_and_rectilinear():
     tgt.add_dim_coord(grid_lon_tgt, 1)
 
     # Change the AuxCoords to check that the DimCoords have priority.
-    src.coord("latitude").bounds[:] = 0
-    src.coord("longitude").bounds[:] = 0
-    tgt.coord("latitude").bounds[:] = 0
-    tgt.coord("longitude").bounds[:] = 0
+    src.coord("latitude").points[:] = 0
+    src.coord("longitude").points[:] = 0
+    tgt.coord("latitude").points[:] = 0
+    tgt.coord("longitude").points[:] = 0
 
     # Ensure the source data is all ones.
     src.data[:] = 1

esmf_regrid/tests/unit/schemes/test_ESMFNearestRegridder.py

@@ -71,7 +71,8 @@ def test_differing_grids():
         _ = regridder(src_dif_circ)
 
 
-def test_curvilinear_equivalence():
+@pytest.mark.parametrize("with_bounds", (False, True))
+def test_curvilinear_equivalence(with_bounds):
     """
     Test initialisation of :class:`esmf_regrid.schemes.ESMFNearestRegridder`.
 
@@ -88,6 +89,11 @@ def test_curvilinear_equivalence():
     grid_tgt = _grid_cube(n_lons_tgt, n_lats_tgt, lon_bounds, lat_bounds, circular=True)
     curv_src = _curvilinear_cube(n_lons_src, n_lats_src, lon_bounds, lat_bounds)
     curv_tgt = _curvilinear_cube(n_lons_tgt, n_lats_tgt, lon_bounds, lat_bounds)
+    if not with_bounds:
+        curv_src.coord("latitude").bounds = None
+        curv_src.coord("longitude").bounds = None
+        curv_tgt.coord("latitude").bounds = None
+        curv_tgt.coord("longitude").bounds = None
 
     grid_to_grid = ESMFNearestRegridder(grid_src, grid_tgt)
     grid_to_curv = ESMFNearestRegridder(grid_src, curv_tgt)
@@ -136,10 +142,10 @@ def test_curvilinear_and_rectilinear():
     tgt.add_dim_coord(grid_lon_tgt, 1)
 
     # Change the AuxCoords to check that the DimCoords have priority.
-    src.coord("latitude").bounds[:] = 0
-    src.coord("longitude").bounds[:] = 0
-    tgt.coord("latitude").bounds[:] = 0
-    tgt.coord("longitude").bounds[:] = 0
+    src.coord("latitude").points[:] = 0
+    src.coord("longitude").points[:] = 0
+    tgt.coord("latitude").points[:] = 0
+    tgt.coord("longitude").points[:] = 0
 
     # Ensure the source data is all ones.
     src.data[:] = 1