Coordinate system agnostic region extraction #6855
Description
✨ Feature Request
I'd like to be able to extract a horizontal region from a cube by specifying a bounding box in some coordinate system (typically normal latitude and longitude) that does not necessarily match the coordinate system of the cube.
Motivation
I often want to extract a horizontal subregion from a larger region. If the cube is in a regular lat-lon coordinate system I can do this without much thought cube.intersection(longitude=(), latitude=()). However, if the coordinate system is something other than regular lat-lon I then need to think about transforming the bounding box. Intersection also does not work if the coordinate does not have a modulus.
Additional context
This is basically my only use for the Met Office internal ascend library, e.g.,
import shapely.geometry as sgeom
import ascend.shape as ashp
bounds = (...) # coordinates of bounding box
cubes = iris.cube.CubeList(...) # some cubes
box = ashp.create(sgeom.box(*bounds), {}, "Polygon")
box.extract_subcubes(cubes)Would be nice if there was equivalent functionality in Iris, or if intersection could be extended to be coordinate system agnostic (and also relax the requirement for the coordinate to have a modulus). I have a function built on top of intersection which basically gets me there, but it would be better if this could just be a simple call to an Iris method/function.
def extract_region(cube: iris.cube.Cube, area: list | tuple | np.ndarray,
crs: cartopy.crs.CRS | None = None) -> iris.cube.Cube:
"""
Extract a horizontal subregion.
Arguments:
area: an iterable of length 2 or 4, containing the (lon, lat)
for a single point, or (min_x, min_y, max_x, max_y) for
a rectangular region.
Keyword Arguments:
crs: the coordinate reference system of the `area` coordinates
provided. If None, :class:`~cartopy.crs.Geodetic` is assumed.
"""
if crs is None:
crs = iris.coord_systems.GeogCS(pp.EARTH_RADIUS).as_cartopy_crs()
else:
if not isinstance(crs, ccrs.CRS):
msg = f"{crs} is not a valid coordinate reference system"
raise TypeError(msg)
cube_crs = cube.coord_system().as_cartopy_crs()
if isinstance(area, Iterable) and not isinstance(area, str):
if len(area) == 2:
x_t, y_t = cube_crs.transform_point(area[0], area[1], crs)
xt_min = xt_max = x_t
yt_min = yt_max = y_t
elif len(area) == 4:
x_min, y_min, x_max, y_max = area
points = np.array(list(itertools.product((x_min, x_max),
(y_min, y_max))))
points_t = cube_crs.transform_points(crs, points[:, 0],
points[:, 1])[:, :2]
xt_min = points_t[:, 0].min()
xt_max = points_t[:, 0].max()
yt_min = points_t[:, 1].min()
yt_max = points_t[:, 1].max()
else:
msg = f"""area should have length 2 or 4, but instead
has length {len(area)}"""
raise ValueError(msg)
for coord in ("x", "y"):
if not cube.coord(axis=coord).has_bounds():
cube.coord(axis=coord).guess_bounds()
x_extent = iris.coords.CoordExtent(cube.coord(axis="x"),
xt_min, xt_max)
y_extent = iris.coords.CoordExtent(cube.coord(axis="y"),
yt_min, yt_max)
try:
return cube.intersection(x_extent, y_extent)
except:
# if coordinate does not have a modulus
x_name = cube.coord(axis="x").name()
y_name = cube.coord(axis="y").name()
coord_values = {x_name: lambda cell: x_extent.minimum <= cell <= x_extent.maximum,
y_name: lambda cell: y_extent.minimum <= cell <= y_extent.maximum}
region_constraint = iris.Constraint(coord_values=coord_values)
return cube.extract(region_constraint)
else:
msg = f"{area} is of type {type(area)} which is not iterable"
raise TypeError(msg)