Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 14 additions & 18 deletions esmf_regrid/experimental/unstructured_scheme.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
import copy

import iris
from iris.analysis._interpolation import get_xy_coords
from iris.analysis._interpolation import get_xy_dim_coords
import numpy as np

# from numpy import ma
Expand All @@ -21,15 +21,6 @@ def _bounds_cf_to_simple_1d(cf_bounds):
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:
# mesh = cube.mesh
# mesh_dim = cube.mesh_dim(mesh)
# return mesh, mesh_dim
pass


def _mesh_to_MeshInfo(mesh):
# Returns a MeshInfo object describing the mesh of the cube.
assert mesh.topology_dimension == 2
Expand Down Expand Up @@ -117,10 +108,14 @@ def _regrid_unstructured_to_rectilinear__prepare(src_mesh_cube, target_grid_cube

# TODO: Record appropriate dimensions (i.e. which dimension the mesh belongs to)

grid_x, grid_y = get_xy_coords(target_grid_cube)
grid_x, grid_y = get_xy_dim_coords(target_grid_cube)
mesh = src_mesh_cube.mesh
# TODO: Improve the checking of mesh validity. Check the mesh location and
# raise appropriate error messages.
assert mesh is not None
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this to be properly handled / exception raised in the future? If so maybe a TODO added?

# From src_mesh_cube, fetch the mesh, and the dimension on the cube which that
# mesh belongs to (mesh_dim).
mesh, mesh_dim = _get_mesh_and_dim(src_mesh_cube)
# mesh belongs to.
mesh_dim = src_mesh_cube.mesh_dim()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Previous pseudocode version of this (L28) took mesh as a parameter. Is it not necessary to do this anymore? Is there a possibility that a cube contains more than one mesh?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The previous pseudocode version was written before this functionality was implemented in iris. We've discussed the idea of a cube containing more than one mesh and we don't anticipate this possibility. In fact I don't believe it's possible to define such a thing under the current UGRID conventions.


meshinfo = _mesh_to_MeshInfo(mesh)
gridinfo = _cube_to_GridInfo(target_grid_cube)
Expand Down Expand Up @@ -183,16 +178,17 @@ def __init__(self, src_mesh_cube, target_grid_cube, mdtol=1):
)

# Store regrid info.
_, _, self.grid_x, self.grid_y, self.regridder = partial_regrid_info
_, self.grid_x, self.grid_y, self.regridder = partial_regrid_info

def __call__(self, cube):
"""TODO: write docstring."""
mesh = cube.mesh
# TODO: Ensure cube has the same mesh as that of the recorded mesh.
# For the time being, we simply check that the mesh exists.
assert mesh is not None
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

TODO here if we plan on adding more error handling

mesh_dim = cube.mesh_dim()

# mesh is probably an iris Mesh object, though it could also be a MeshCoord
mesh, mesh_dim = _get_mesh_and_dim(cube)

regrid_info = (mesh, mesh_dim, self.grid_x, self.grid_y, self.regridder)
regrid_info = (mesh_dim, self.grid_x, self.grid_y, self.regridder)

return _regrid_unstructured_to_rectilinear__perform(
cube, regrid_info, self.mdtol
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
"""Unit tests for :mod:`esmf_regrid.experimental.unstructured_scheme`."""
Original file line number Diff line number Diff line change
@@ -0,0 +1,58 @@
"""Unit tests for :func:`esmf_regrid.experimental.unstructured_scheme.MeshToGridESMFRegridder`."""

from iris.coords import AuxCoord, DimCoord
import numpy as np

from esmf_regrid.experimental.unstructured_scheme import (
MeshToGridESMFRegridder,
)
from esmf_regrid.tests.unit.experimental.unstructured_scheme.test__cube_to_GridInfo import (
_grid_cube,
)
from esmf_regrid.tests.unit.experimental.unstructured_scheme.test__regrid_unstructured_to_rectilinear__prepare import (
_flat_mesh_cube,
)


def test_flat_cubes():
"""
Basic test for :func:`esmf_regrid.experimental.unstructured_scheme.MeshToGridESMFRegridder`.

Tests with flat cubes as input (a 1D mesh cube and a 2D grid cube).
"""
src = _flat_mesh_cube()

n_lons = 6
n_lats = 5
lon_bounds = (-180, 180)
lat_bounds = (-90, 90)
tgt = _grid_cube(n_lons, n_lats, lon_bounds, lat_bounds, circular=True)
# Ensure data in the target grid is different to the expected data.
# i.e. target grid data is all zero, expected data is all one
tgt.data[:] = 0

def _add_metadata(cube):
result = cube.copy()
result.units = "K"
result.attributes = {"a": 1}
result.standard_name = "air_temperature"
scalar_height = AuxCoord([5], units="m", standard_name="height")
scalar_time = DimCoord([10], units="s", standard_name="time")
result.add_aux_coord(scalar_height)
result.add_aux_coord(scalar_time)
return result

src = _add_metadata(src)
src.data[:] = 1 # Ensure all data in the source is one.
regridder = MeshToGridESMFRegridder(src, tgt)
result = regridder(src)

expected_data = np.ones([n_lats, n_lons])
expected_cube = _add_metadata(tgt)

# Lenient check for data.
assert np.allclose(expected_data, result.data)

# Check metadata and scalar coords.
expected_cube.data = result.data
assert expected_cube == result
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
"""Unit tests for :func:`esmf_regrid.experimental.unstructured_scheme._regrid_unstructured_to_rectilinear__prepare`."""

from iris.coords import AuxCoord
from iris.cube import Cube
import numpy as np

from esmf_regrid.esmf_regridder import GridInfo
from esmf_regrid.experimental.unstructured_regrid import MeshInfo
from esmf_regrid.experimental.unstructured_scheme import (
_regrid_unstructured_to_rectilinear__prepare,
)
from esmf_regrid.tests.unit.experimental.unstructured_scheme.test__cube_to_GridInfo import (
_grid_cube,
)
from esmf_regrid.tests.unit.experimental.unstructured_scheme.test__mesh_to_MeshInfo import (
_example_mesh,
)


def _full_mesh():
mesh = _example_mesh()

# In order to add a mesh to a cube, face locations must be added.
# These are not used in calculations and are here given a value of zero.
mesh_length = mesh.connectivity(contains_face=True).shape[0]
dummy_face_lon = AuxCoord(np.zeros(mesh_length), standard_name="longitude")
dummy_face_lat = AuxCoord(np.zeros(mesh_length), standard_name="latitude")
mesh.add_coords(face_x=dummy_face_lon, face_y=dummy_face_lat)
mesh.long_name = "example mesh"
return mesh


def _flat_mesh_cube():
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you document what the object that this returns looks like? i.e. data is all ones?

"""
Return a 1D cube with a mesh attached.

Returned cube has no metadata except for the mesh and two MeshCoords.
Returned cube has data consisting of an array of ones.
"""
mesh = _full_mesh()
mesh_length = mesh.connectivity(contains_face=True).shape[0]

cube = Cube(np.ones([mesh_length]))
mesh_coord_x, mesh_coord_y = mesh.to_MeshCoords("face")
cube.add_aux_coord(mesh_coord_x, 0)
cube.add_aux_coord(mesh_coord_y, 0)
return cube


def test_flat_cubes():
"""
Basic test for :func:`esmf_regrid.experimental.unstructured_scheme._regrid_unstructured_to_rectilinear__prepare`.

Tests with flat cubes as input (a 1D mesh cube and a 2D grid cube).
"""
src = _flat_mesh_cube()

n_lons = 6
n_lats = 5
lon_bounds = (-180, 180)
lat_bounds = (-90, 90)
tgt = _grid_cube(n_lons, n_lats, lon_bounds, lat_bounds, circular=True)
regrid_info = _regrid_unstructured_to_rectilinear__prepare(src, tgt)
mesh_dim, grid_x, grid_y, regridder = regrid_info

assert mesh_dim == 0
assert grid_x == tgt.coord("longitude")
assert grid_y == tgt.coord("latitude")
assert type(regridder.tgt) == GridInfo
assert type(regridder.src) == MeshInfo
Original file line number Diff line number Diff line change
@@ -0,0 +1,57 @@
"""Unit tests for :func:`esmf_regrid.experimental.unstructured_scheme.regrid_unstructured_to_rectilinear`."""

from iris.coords import AuxCoord, DimCoord
import numpy as np

from esmf_regrid.experimental.unstructured_scheme import (
regrid_unstructured_to_rectilinear,
)
from esmf_regrid.tests.unit.experimental.unstructured_scheme.test__cube_to_GridInfo import (
_grid_cube,
)
from esmf_regrid.tests.unit.experimental.unstructured_scheme.test__regrid_unstructured_to_rectilinear__prepare import (
_flat_mesh_cube,
)


def test_flat_cubes():
"""
Basic test for :func:`esmf_regrid.experimental.unstructured_scheme.regrid_unstructured_to_rectilinear`.

Tests with flat cubes as input (a 1D mesh cube and a 2D grid cube).
"""
src = _flat_mesh_cube()

n_lons = 6
n_lats = 5
lon_bounds = (-180, 180)
lat_bounds = (-90, 90)
tgt = _grid_cube(n_lons, n_lats, lon_bounds, lat_bounds, circular=True)
# Ensure data in the target grid is different to the expected data.
# i.e. target grid data is all zero, expected data is all one
tgt.data[:] = 0

def _add_metadata(cube):
result = cube.copy()
result.units = "K"
result.attributes = {"a": 1}
result.standard_name = "air_temperature"
scalar_height = AuxCoord([5], units="m", standard_name="height")
scalar_time = DimCoord([10], units="s", standard_name="time")
result.add_aux_coord(scalar_height)
result.add_aux_coord(scalar_time)
return result

src = _add_metadata(src)
src.data[:] = 1 # Ensure all data in the source is one.
result = regrid_unstructured_to_rectilinear(src, tgt)

expected_data = np.ones([n_lats, n_lons])
expected_cube = _add_metadata(tgt)

# Lenient check for data.
assert np.allclose(expected_data, result.data)

# Check metadata and scalar coords.
expected_cube.data = result.data
assert expected_cube == result