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.

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

Bugfix shape masking #6129

Draft
wants to merge 19 commits into
base: main
Choose a base branch
from
Draft
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
389 changes: 232 additions & 157 deletions lib/iris/_shapefiles.py

Large diffs are not rendered by default.

5 changes: 5 additions & 0 deletions lib/iris/tests/unit/_shapefiles/__init__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
# Copyright Iris contributors
#
# This file is part of Iris and is released under the BSD license.
# See LICENSE in the root of the repository for full licensing details.
"""Unit tests for the :mod:`iris._shapefiles` module."""
160 changes: 160 additions & 0 deletions lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,160 @@
# Copyright Iris contributors
#
# This file is part of Iris and is released under the BSD license.
# See LICENSE in the root of the repository for full licensing details.
"""Unit tests for :func:`iris._shapefiles.is_geometry_valid`."""

# import iris tests first so that some things can be initialised before
# importing anything else
import iris.tests as tests # isort:skip

from pyproj import CRS
import pytest
from shapely.geometry import box

from iris._shapefiles import is_geometry_valid


# Shareable shape fixtures used in:
# - util/test_mask_cube_from_shapefile.py
# - _shapefiles/test_is_geometry_valid.py
@pytest.fixture(scope="session")
def wgs84_crs():
return CRS.from_epsg(4326)


@pytest.fixture(scope="session")
def osgb_crs():
return CRS.from_epsg(27700)


@pytest.fixture(scope="session")
def basic_circular_geometry():
# Define geometry of a basic circle with center at (0, 0) and radius 10
center = (0, 0)
radius = 10

circle = Point(center).buffer(radius)

return circle


@pytest.fixture(scope="session")
def basic_rectangular_geometry():
# Define the coordinates of a basic rectangle
min_lon = -90
min_lat = -45
max_lon = 90
max_lat = 45

# Create the rectangular geometry
return box(min_lon, min_lat, max_lon, max_lat)


@pytest.fixture(scope="session")
def basic_point_geometry():
# Define the coordinates of a basic point (lon, lat)
return Point((-3.476204, 50.727059))


@pytest.fixture(scope="session")
def basic_line_geometry():
# Define the coordinates of a basic line
return LineString([(0, 0), (10, 10)])


@pytest.fixture(scope="session")
def basic_point_collection():
# Define the coordinates of a basic collection of points
# as (lon, lat) tuples, assuming a WGS84 projection.
points = MultiPoint(
[
(0, 0),
(10, 10),
(-10, -10),
(-3.476204, 50.727059),
(174.761067, -36.846211),
(-77.032801, 38.892717),
]
)

return points


@pytest.fixture(scope="session")
def canada_geometry():
# Define the coordinates of a rectangle that covers Canada
return box(-143.5, 42.6, -37.8, 84.0)


@pytest.fixture(scope="session")
def bering_sea_geometry():
# Define the coordinates of a rectangle that covers the Bering Sea
return box(148.42, 49.1, -138.74, 73.12)


@pytest.fixture(scope="session")
def uk_geometry():
# Define the coordinates of a rectangle that covers the UK
return box(-10, 49, 2, 61)


@pytest.fixture(scope="session")
def invalid_geometry_poles():
# Define the coordinates of a rectangle that crosses the poles
return box(-10, -90, 10, 90)


@pytest.fixture(scope="session")
def invalid_geometry_bounds():
# Define the coordinates of a rectangle that is outside the bounds of the coordinate system
return box(-200, -100, 200, 100)


class TestGeometry(tests.IrisTest):
# Test validity of different geometries
@pytest.mark.parametrize(
"test_input",
[
"basic_circular_geometry",
"basic_rectangular_geometry",
"basic_point_geometry",
"basic_line_geometry",
"basic_point_collection",
"canada_geometry",
],
)
def test_valid_geometry(test_input, expected):
# Assert that all valid geometries are return None
assert is_geometry_valid(request.getfixturevalue(test_input), wgs84) is None

# Fixtures retrieved from conftest.py
@pytest.mark.parametrize(
"test_input, errortype, errormessage",
[
(
"bering_sea_geometry",
ValueError,
"Geometry crossing the antimeridian is not supported.",
),
(
"invalid_geometry_poles",
ValueError,
"Geometry crossing the poles is not supported.",
),
(
"invalid_geometry_bounds",
ValueError,
"Geometry [<POLYGON ((200 -100, 200 100, -200 100, -200 -100, 200 -100))>] is not valid for the given coordinate system EPSG:4326. Check that your coordinates are correctly specified.",
),
(
"not a valid geometry",
TypeError,
"Geometry is not a valid Shapely object",
),
],
)
def test_invalid_geometry(test_input, errortype, errormessage):
# Assert that all invalid geometries raise the expected error
with pytest.raises(errortype, match=errormessage):
is_geometry_valid(request.getfixturevalue(test_input), wgs84)
1 change: 1 addition & 0 deletions lib/iris/tests/unit/conftest.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
# See LICENSE in the root of the repository for full licensing details.
"""Unit tests fixture infra-structure."""

from pyproj import CRS
import pytest

import iris
Expand Down
18 changes: 18 additions & 0 deletions lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,16 +4,34 @@
# See LICENSE in the root of the repository for full licensing details.
"""Unit tests for :func:`iris.util.mask_cube_from_shapefile`."""

from cartopy.feature import ShapelyFeature
import numpy as np
from pyproj import CRS
import pytest
import shapely
from shapely import LineString, MultiPoint, Point
from shapely.geometry import box

from iris.coord_systems import RotatedGeogCS
from iris.coords import DimCoord
import iris.cube
from iris.util import mask_cube_from_shapefile
from iris.warnings import IrisUserWarning

# def _apply_WGS84_crs(geometry):
# """Apply a WGS84 CRS to a geometry."""
# crs = CRS.from_epsg(4326)
# return ShapelyFeature(geom

# def _apply_OSGB36_crs(geometry):
# """Apply an OSGB36 CRS to a geometry."""
# crs = CRS.from_epsg(27700)
# geometry = geometry.__geo_interface__
# geometry["crs"] = crs
# return geometry

# Fixtures retrieved from _shapefiles/test_is_geometry_valid.py


class TestBasicCubeMasking:
"""Unit tests for mask_cube_from_shapefile function."""
Expand Down
96 changes: 82 additions & 14 deletions lib/iris/util.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
import os.path
import sys
import tempfile
from typing import Literal
from typing import Literal, overload
from warnings import warn

import cf_units
Expand Down Expand Up @@ -2129,27 +2129,52 @@ def _strip_metadata_from_dims(cube, dims):
return reduced_cube


def mask_cube_from_shapefile(cube, shape, minimum_weight=0.0, in_place=False):
"""Take a shape object and masks all points not touching it in a cube.
@overload
def mask_cube_from_shapefile(
cube: iris.cube.Cube,
shape: shapely.Geometry,
shape_crs: cartopy.crs | CRS,
in_place: bool = False,
all_touched: bool = False,
invert: bool = False,
): ...
def mask_cube_from_shapefile(
cube: iris.cube.Cube,
shape: shapely.Geometry,
shape_crs: cartopy.crs | CRS,
in_place: bool = False,
**kwargs,
):
"""Mask all points in a cube that do not intersect a shapefile object.

Mask a :class:`~iris.cube.Cube` with any shapefile object, (e.g. Natural Earth Shapefiles via cartopy).
Finds the overlap between the `shapefile` shape and the :class:`~iris.cube.Cube` and
masks out any cells that __do not__ intersect the shape.

Shapes can be polygons, lines or points.

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
By default, only the only cells whose center is within the polygon or that are selected by
Bresenham’s line algorithm (for line type shape) are kept. This behaviour can be changed by
using the `all_touched=True` keyword argument. Then all cells touched by geometries are kept.

Parameters
----------
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.
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.
shape_crs : cartopy.crs.CRS, default=None
The coordinate reference system of the shape object.
in_place : bool, default=False
Whether to mask the `cube` in-place or return a newly masked `cube`.
Defaults to False.
**kwargs
Additional keyword arguments to pass to `rasterio.features.geometry_mask`.
Valid keyword arguments are:
all_touched : bool, optional
invert : bool, optional

Returns
-------
Expand All @@ -2163,8 +2188,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 @@ -2176,13 +2199,58 @@ def mask_cube_from_shapefile(cube, shape, minimum_weight=0.0, in_place=False):
Examples
--------
>>> import shapely
>>> from pyproj import CRS
>>> from iris.util import mask_cube_from_shapefile

Extract a rectangular region covering the UK from a stereographic projection cube:

>>> cube = iris.load_cube(iris.sample_data_path("toa_brightness_stereographic.nc"))
>>> shape = shapely.geometry.box(-10, 50, 2, 60) # box around the UK
>>> wgs84 = CRS.from_epsg(4326) # WGS84 coordinate system
>>> masked_cube = mask_cube_from_shapefile(cube, shape, wgs84)

Extract a trajectory by using a line shapefile:

>>> from shapely import LineString
>>> line = LineString([(-45, 40), (-28, 53), (-2, 55), (19, 45)])
>>> masked_cube = mask_cube_from_shapefile(cube, line, wgs84)

Standard shapely manipulations can be applied. For example, to extract a trajectory
with a 1 degree buffer around it:

>>> buffer = line.buffer(1)
>>> masked_cube = mask_cube_from_shapefile(cube, buffer, wgs84)

You can load more complex shapes from other libraries. For example, to extract the
Canadian provience of Ontario from a cube:

>>> import cartopy.io.shapereader as shpreader
>>> admin1 = shpreader.natural_earth(resolution='110m',
category='cultural',
name='admin_1_states_provinces_lakes')
>>> admin1shp = shpreader.Reader(admin1).geometries()

>>> cube = iris.load_cube(iris.sample_data_path("E1_north_america.nc"))
>>> shape = shapely.geometry.box(-100,30, -80,40) # box between 30N-40N 100W-80W
>>> masked_cube = mask_cube_from_shapefile(cube, shape)
>>> wgs84 = CRS.from_epsg(4326)
>>> masked_cube = mask_cube_from_shapefile(cube, shape, wgs84)

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(
geometry=shape, geometry_crs=shape_crs, cube=cube, **kwargs
)
masked_cube = mask_cube(cube, shapefile_mask, in_place=in_place)
if not in_place:
return masked_cube
Expand Down
3 changes: 3 additions & 0 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,9 @@ license = {text = "BSD-3-Clause"}
name = "scitools-iris"
requires-python = ">=3.10"

[project.optional-dependencies]
RASTERIO = ["rasterio"]

[project.urls]
Code = "https://github.com/SciTools/iris"
Discussions = "https://github.com/SciTools/iris/discussions"
Expand Down
1 change: 1 addition & 0 deletions requirements/py312.yml
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ dependencies:
- pandas
- pip
- python-stratify
- rasterio

# Test dependencies.
- asv_runner
Expand Down
Loading