Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Bugfix shape masking #6129

Draft
wants to merge 4 commits into
base: main
Choose a base branch
from
Draft

Bugfix shape masking #6129

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
ebee977
Updates to allow varying shape crs
hsteptoe Aug 22, 2024
e32f911
Ruff fixes
hsteptoe Aug 22, 2024
76eae38
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Aug 22, 2024
bcdca30
Merge branch 'SciTools:main' into bugfix-shape-masking
hsteptoe Nov 22, 2024
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
86 changes: 56 additions & 30 deletions lib/iris/_shapefiles.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,14 +17,16 @@
import shapely.geometry as sgeom
import shapely.ops

import iris.analysis.cartography
from iris.warnings import IrisDefaultingWarning, IrisUserWarning


def create_shapefile_mask(
geometry,
cube,
minimum_weight=0.0,
):
geometry: shapely.Geometry,
cube: iris.cube.Cube,
minimum_weight: float = 0.0,
geometry_crs: cartopy.crs = None,
) -> np.array:
"""Make a mask for a cube from a shape.

Get the mask of the intersection between the
Expand All @@ -43,13 +45,29 @@ def create_shapefile_mask(
to be unmasked.
Requires geometry to be a Polygon or MultiPolygon
Defaults to 0.0 (eg only test intersection).
geometry_crs : :class:`cartopy.crs`, optional
A :class:`~iris.coord_systems` object describing
the coord_system of the shapefile. Defaults to None,
in which case the geometry_crs is assumed to be the
same as the `cube`.

Returns
-------
:class:`np.array`
An array of the shape of the x & y coordinates of the cube, with points
to mask equal to True.

Notes
-----
For best masking results, both the cube _and_ masking geometry should have a
coordinate reference system (CRS) defined. Masking results will be most reliable
when the cube and masking geometry have the same CRS.

If the cube has no coord_system, the default GeogCS is used where
the coordinate units are degrees. For any other coordinate units,
the cube **must** have a coord_system defined.

If a CRS is not provided for the the masking geometry, the CRS of the cube is assumed.
"""
from iris.cube import Cube, CubeList

Expand Down Expand Up @@ -116,7 +134,9 @@ def create_shapefile_mask(
return mask_template


def _transform_coord_system(geometry, cube, geometry_system=None):
def _transform_coord_system(
geometry: shapely.Geometry, cube: iris.cube.Cube, geometry_crs: cartopy.crs = None
) -> shapely.Geometry:
"""Project the shape onto another coordinate system.

Parameters
Expand All @@ -125,36 +145,52 @@ def _transform_coord_system(geometry, cube, geometry_system=None):
cube : :class:`iris.cube.Cube`
:class:`~iris.cube.Cube` with the coord_system to be projected to and
a x coordinate.
geometry_system : :class:`iris.coord_systems`, optional
A :class:`~iris.coord_systems` object describing
geometry_crs : :class:`cartopy.crs`, optional
A :class:`cartopy.crs` object describing
the coord_system of the shapefile. Defaults to None,
which is treated as GeogCS.
in which case the geometry_crs is assumed to be the
same as the `cube`.

Returns
-------
:class:`shapely.Geometry`
A transformed copy of the provided :class:`shapely.Geometry`.

"""
y_name, x_name = _cube_primary_xy_coord_names(cube)
import iris.analysis.cartography
_, x_name = _cube_primary_xy_coord_names(cube)

DEFAULT_CS = iris.coord_systems.GeogCS(
iris.analysis.cartography.DEFAULT_SPHERICAL_EARTH_RADIUS
)
target_system = cube.coord_system()
target_system = cube.coord_system().as_cartopy_projection()
if not target_system:
# If no cube coord_system do our best to guess...
if (
cube.coord(axis="x").units == "degrees"
and cube.coord(axis="y").units == "degrees"
):
# If units of degrees assume GeogCS
target_system = iris.coord_systems.GeogCS(
iris.analysis.cartography.DEFAULT_SPHERICAL_EARTH_RADIUS
)
warnings.warn(
"Cube has no coord_system; using default GeogCS lat/lon.",
category=IrisDefaultingWarning,
)
else:
# For any other units, don't guess and raise an error
target_system = iris.coord_systems.GeogCS(
iris.analysis.cartography.DEFAULT_SPHERICAL_EARTH_RADIUS
)
raise ValueError("Cube has no coord_system; cannot guess coord_system!")

if not geometry_crs:
# If no geometry_crs assume it has the cube coord_system
geometry_crs = target_system
warnings.warn(
"Cube has no coord_system; using default GeogCS lat/lon",
"No geometry coordinate reference system supplied; using cube coord_system instead.",
category=IrisDefaultingWarning,
)
target_system = DEFAULT_CS
if geometry_system is None:
geometry_system = DEFAULT_CS
target_proj = target_system.as_cartopy_projection()
source_proj = geometry_system.as_cartopy_projection()

trans_geometry = target_proj.project_geometry(geometry, source_proj)

# A GeogCS in iris can be either -180 to 180 or 0 to 360. If cube is 0-360, shift geom to match
if (
isinstance(target_system, iris.coord_systems.GeogCS)
Expand All @@ -165,16 +201,6 @@ def _transform_coord_system(geometry, cube, geometry_system=None):
trans_geometry = trans_geometry.difference(prime_meridian_line.buffer(0.00001))
trans_geometry = shapely.transform(trans_geometry, _trans_func)

if (not isinstance(target_system, iris.coord_systems.GeogCS)) and cube.coord(
x_name
).points[-1] > 180:
# this may lead to incorrect masking or not depending on projection type so warn user
warnings.warn(
"""Cube has x-coordinates over 180E and a non-standard projection type.\n
This may lead to incorrect masking. \n
If the result is not as expected, you might want to transform the x coordinate points of your cube to -180-180 """,
category=IrisUserWarning,
)
return trans_geometry


Expand Down
30 changes: 25 additions & 5 deletions lib/iris/util.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2132,9 +2132,17 @@ def _strip_metadata_from_dims(cube, dims):
return reduced_cube


def mask_cube_from_shapefile(cube, shape, minimum_weight=0.0, in_place=False):
def mask_cube_from_shapefile(
cube: iris.cube.Cube,
shape: shapely.Geometry,
shape_crs: cartopy.crs = None,
minimum_weight: float = 0.0,
in_place: bool = False,
):
"""Take a shape object and masks all points not touching it in a cube.

This function allows masking a cube with any shape object,
(e.g. Natural Earth Shapefiles via cartopy).
Finds the overlap between the `shape` and the `cube` in 2D xy space and
masks out any cells with less % overlap with shape than set.
Default behaviour is to count any overlap between shape and cell as valid
Expand All @@ -2143,10 +2151,12 @@ def mask_cube_from_shapefile(cube, shape, minimum_weight=0.0, in_place=False):
----------
cube : :class:`~iris.cube.Cube` object
The `Cube` object to masked. Must be singular, rather than a `CubeList`.
shape : Shapely.Geometry object
shape : shapely.Geometry object
A single `shape` of the area to remain unmasked on the `cube`.
If it a line object of some kind then minimum_weight will be ignored,
because you cannot compare the area of a 1D line and 2D Cell.
shape_crs : cartopy.crs.CRS, default=None
The coordinate reference system of the shape object.
minimum_weight : float , default=0.0
A number between 0-1 describing what % of a cube cell area must
the shape overlap to include it.
Expand All @@ -2166,8 +2176,6 @@ def mask_cube_from_shapefile(cube, shape, minimum_weight=0.0, in_place=False):

Notes
-----
This function allows masking a cube with any cartopy projection by a shape object,
most commonly from Natural Earth Shapefiles via cartopy.
To mask a cube from a shapefile, both must first be on the same coordinate system.
Shapefiles are mostly on a lat/lon grid with a projection very similar to GeogCS
The shapefile is projected to the coord system of the cube using cartopy, then each cell
Expand All @@ -2184,8 +2192,20 @@ def mask_cube_from_shapefile(cube, shape, minimum_weight=0.0, in_place=False):
>>> shape = shapely.geometry.box(-100,30, -80,40) # box between 30N-40N 100W-80W
>>> masked_cube = mask_cube_from_shapefile(cube, shape)

Warning
-------
For best masking results, both the cube _and_ masking geometry should have a
coordinate reference system (CRS) defined. Masking results will be most reliable
when the cube and masking geometry have the same CRS.

If the cube has no coord_system, the default GeogCS is used where
the coordinate units are degrees. For any other coordinate units,
the cube **must** have a coord_system defined.

If a CRS is not provided for the the masking geometry, the CRS of the cube is assumed.

"""
shapefile_mask = create_shapefile_mask(shape, cube, minimum_weight)
shapefile_mask = create_shapefile_mask(shape, cube, minimum_weight, shape_crs)
masked_cube = mask_cube(cube, shapefile_mask, in_place=in_place)
if not in_place:
return masked_cube
Loading