Bugfix shape masking

lib/iris/tests/unit/_shapefiles/__init__.py

@@ -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."""

lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py

@@ -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)

lib/iris/tests/unit/conftest.py

@@ -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

lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py

@@ -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."""

lib/iris/util.py

@@ -2129,9 +2129,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
@@ -2140,10 +2148,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.
@@ -2163,8 +2173,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
@@ -2181,8 +2189,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

pyproject.toml

@@ -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"

requirements/py312.yml

@@ -34,6 +34,7 @@ dependencies:
   - pandas
   - pip
   - python-stratify
+  - rasterio
 
 # Test dependencies.
   - asv_runner