diff --git a/CHANGELOG.md b/CHANGELOG.md index be23b5ba..fd53b5c5 100644 --- a/CHANGELOG.md +++ b/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) diff --git a/docs/src/userguide/examples.rst b/docs/src/userguide/examples.rst index 06564656..e66a5a14 100644 --- a/docs/src/userguide/examples.rst +++ b/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) diff --git a/esmf_regrid/_esmf_sdo.py b/esmf_regrid/_esmf_sdo.py index ee657908..2da3128e 100644 --- a/esmf_regrid/_esmf_sdo.py +++ b/esmf_regrid/_esmf_sdo.py @@ -162,7 +162,7 @@ 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 " @@ -170,7 +170,7 @@ def __init__( ) 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,26 +219,31 @@ 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: @@ -246,8 +251,12 @@ def _as_esmf_info(self): # 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) diff --git a/esmf_regrid/experimental/unstructured_scheme.py b/esmf_regrid/experimental/unstructured_scheme.py index 07922235..172bfd4e 100644 --- a/esmf_regrid/experimental/unstructured_scheme.py +++ b/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 diff --git a/esmf_regrid/schemes.py b/esmf_regrid/schemes.py index 7452c4bb..69ba70e1 100644 --- a/esmf_regrid/schemes.py +++ b/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 diff --git a/esmf_regrid/tests/unit/schemes/test_ESMFBilinearRegridder.py b/esmf_regrid/tests/unit/schemes/test_ESMFBilinearRegridder.py index 519eeac6..da60f416 100644 --- a/esmf_regrid/tests/unit/schemes/test_ESMFBilinearRegridder.py +++ b/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 diff --git a/esmf_regrid/tests/unit/schemes/test_ESMFNearestRegridder.py b/esmf_regrid/tests/unit/schemes/test_ESMFNearestRegridder.py index 46540ddf..6c325e4d 100644 --- a/esmf_regrid/tests/unit/schemes/test_ESMFNearestRegridder.py +++ b/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 diff --git a/esmf_regrid/tests/unit/schemes/test__cube_to_GridInfo.py b/esmf_regrid/tests/unit/schemes/test__cube_to_GridInfo.py index fa74d829..386a562c 100644 --- a/esmf_regrid/tests/unit/schemes/test__cube_to_GridInfo.py +++ b/esmf_regrid/tests/unit/schemes/test__cube_to_GridInfo.py @@ -8,6 +8,7 @@ import pytest import scipy.sparse +from esmf_regrid import Constants from esmf_regrid.esmf_regridder import Regridder from esmf_regrid.schemes import _contiguous_masked, _cube_to_GridInfo @@ -200,6 +201,56 @@ def test_curvilinear_grid(): assert np.allclose(expected_weights.todense(), rg_circular.weight_matrix.todense()) +@pytest.mark.parametrize("curvilinear", [True, False]) +def test_center(curvilinear): + """Test the beheaviour when center=True.""" + n_lons = 6 + n_lats = 5 + lon_bounds = (-180, 180) + lat_bounds = (-90, 90) + + if curvilinear: + cube = _curvilinear_cube( + n_lons, + n_lats, + lon_bounds, + lat_bounds, + coord_system=GeogCS(EARTH_RADIUS), + ) + else: + cube = _grid_cube( + n_lons, + n_lats, + lon_bounds, + lat_bounds, + coord_system=GeogCS(EARTH_RADIUS), + ) + cube.coord("latitude").bounds = None + cube.coord("longitude").bounds = None + gridinfo = _cube_to_GridInfo(cube, center=True) + # Ensure conversion to ESMF works without error + _ = gridinfo.make_esmf_field() + + # The following test ensures there are no overlapping cells. + # This catches geometric/topological abnormalities that would arise from, + # for example: switching lat/lon values, using euclidean coords vs spherical. + rg = Regridder(gridinfo, gridinfo, method=Constants.Method.BILINEAR) + expected_weights = scipy.sparse.identity(n_lats * n_lons) + assert np.array_equal(expected_weights.todense(), rg.weight_matrix.todense()) + assert gridinfo.crs == GeogCS(EARTH_RADIUS).as_cartopy_crs() + + # While curvilinear coords do not have the "circular" attribute, the code + # allows "circular" to be True when setting the core regridder directly. + # This describes an ESMF object which is topologically different, but ought + # to be geometrically equivalent to the non-circular case. + circular_gridinfo = _cube_to_GridInfo(cube, center=True) + circular_gridinfo.circular = True + rg_circular = Regridder( + circular_gridinfo, gridinfo, method=Constants.Method.BILINEAR + ) + assert np.allclose(expected_weights.todense(), rg_circular.weight_matrix.todense()) + + @pytest.mark.parametrize("masked", (True, False), ids=("masked", "unmasked")) def test__contiguous_bounds(masked): """Test generation of contiguous bounds."""