Unstructured Scheme - Basic GridInfo Handling

esmf_regrid/experimental/unstructured_scheme.py

@@ -2,15 +2,24 @@
 
 import iris
 from iris.analysis._interpolation import get_xy_coords
+import numpy as np
 
 # from numpy import ma
 
-from esmf_regrid.esmf_regridder import Regridder
+from esmf_regrid.esmf_regridder import GridInfo, Regridder
 
-# from esmf_regrid.esmf_regridder import GridInfo
 # from esmf_regrid.experimental.unstructured_regrid import MeshInfo
 
 
+# Taken from PR #26
+def _bounds_cf_to_simple_1d(cf_bounds):
+    assert (cf_bounds[1:, 0] == cf_bounds[:-1, 1]).all()
+    simple_bounds = np.empty((cf_bounds.shape[0] + 1,), dtype=np.float64)
+    simple_bounds[:-1] = cf_bounds[:, 0]
+    simple_bounds[-1] = cf_bounds[-1, 1]
+    return simple_bounds
+
+
 def _get_mesh_and_dim(cube):
     # Returns the cube's mesh and the dimension that mesh belongs to.
     # Likely to be of the form:
@@ -26,9 +35,26 @@ def _cube_to_MeshInfo(cube):
 
 
 def _cube_to_GridInfo(cube):
+    # This is a simplified version of an equivalent function/method in PR #26.
+    # It is anticipated that this function will be replaced by the one in PR #26.
+    #
     # Returns a GridInfo object describing the horizontal grid of the cube.
     # This may be inherited from code written for the rectilinear regridding scheme.
-    pass
+    lon = cube.coord("longitude")
+    lat = cube.coord("latitude")
+    # Ensure coords come from a proper grid.
+    assert isinstance(lon, iris.coords.DimCoord)
+    assert isinstance(lat, iris.coords.DimCoord)
+    # TODO: accomodate other x/y coords.
+    # TODO: perform checks on lat/lon.
+    #  Checks may cover units, coord systems (e.g. rotated pole), contiguous bounds.
+    return GridInfo(
+        lon.points,
+        lat.points,
+        _bounds_cf_to_simple_1d(lon.bounds),
+        _bounds_cf_to_simple_1d(lat.bounds),
+        circular=lon.circular,
+    )
 
 
 # def _regrid_along_dims(regridder, data, src_dim, mdtol):

esmf_regrid/tests/unit/experimental/unstructured_scheme/test__cube_to_GridInfo.py

@@ -0,0 +1,94 @@
+"""Unit tests for miscellaneous helper functions in `esmf_regrid.experimental.unstructured_scheme`."""
+
+from iris.coords import DimCoord
+from iris.cube import Cube
+import numpy as np
+import scipy.sparse
+
+from esmf_regrid.esmf_regridder import Regridder
+from esmf_regrid.experimental.unstructured_scheme import _cube_to_GridInfo
+
+
+def _generate_points_and_bounds(n, outer_bounds):
+    lower, upper = outer_bounds
+    full_span = np.linspace(lower, upper, n * 2 + 1)
+    points = full_span[1::2]
+    bound_span = full_span[::2]
+    bounds = np.stack([bound_span[:-1], bound_span[1:]], axis=-1)
+    return points, bounds
+
+
+def _grid_cube(n_lons, n_lats, lon_outer_bounds, lat_outer_bounds, circular=False):
+    lon_points, lon_bounds = _generate_points_and_bounds(n_lons, lon_outer_bounds)
+    lon = DimCoord(
+        lon_points, "longitude", units="degrees", bounds=lon_bounds, circular=circular
+    )
+    lat_points, lat_bounds = _generate_points_and_bounds(n_lats, lat_outer_bounds)
+    lat = DimCoord(lat_points, "latitude", units="degrees", bounds=lat_bounds)
+
+    data = np.zeros([n_lats, n_lons])
+    cube = Cube(data)
+    cube.add_dim_coord(lon, 1)
+    cube.add_dim_coord(lat, 0)
+    return cube
+
+
+def test_global_grid():
+    """Test conversion of a global grid."""
+    n_lons = 6
+    n_lats = 5
+    lon_bounds = (-180, 180)
+    lat_bounds = (-90, 90)
+
+    cube = _grid_cube(n_lons, n_lats, lon_bounds, lat_bounds, circular=True)
+    gridinfo = _cube_to_GridInfo(cube)
+    # 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)
+    expected_weights = scipy.sparse.identity(n_lats * n_lons)
+    assert np.array_equal(expected_weights.todense(), rg.weight_matrix.todense())
+
+
+def test_local_grid():
+    """Test conversion of a local grid."""
+    n_lons = 6
+    n_lats = 5
+    lon_bounds = (-20, 20)
+    lat_bounds = (20, 60)
+
+    cube = _grid_cube(n_lons, n_lats, lon_bounds, lat_bounds)
+    gridinfo = _cube_to_GridInfo(cube)
+    # Ensure conversion to ESMF works without error
+    _ = gridinfo.make_esmf_field()
+
+    # The following test ensures there are no overlapping cells.
+    # Note that this test fails when longitude is circular.
+    rg = Regridder(gridinfo, gridinfo)
+    expected_weights = scipy.sparse.identity(n_lats * n_lons)
+    assert np.array_equal(expected_weights.todense(), rg.weight_matrix.todense())
+
+
+def test_grid_with_scalars():
+    """Test conversion of a grid with scalar coords."""
+    n_lons = 1
+    n_lats = 5
+    lon_bounds = (-20, 20)
+    lat_bounds = (20, 60)
+
+    cube = _grid_cube(n_lons, n_lats, lon_bounds, lat_bounds)
+    # Convert longitude to a scalar
+    cube = cube[:, 0]
+    assert len(cube.shape) == 1
+
+    gridinfo = _cube_to_GridInfo(cube)
+    # Ensure conversion to ESMF works without error
+    _ = gridinfo.make_esmf_field()
+
+    # The following test ensures there are no overlapping cells.
+    rg = Regridder(gridinfo, gridinfo)
+    expected_weights = scipy.sparse.identity(n_lats * n_lons)
+    assert np.array_equal(expected_weights.todense(), rg.weight_matrix.todense())