From ebee97757122d257aab1ab35a0f412d1de0016e2 Mon Sep 17 00:00:00 2001 From: Hamish Steptoe Date: Thu, 22 Aug 2024 15:05:12 +0100 Subject: [PATCH 01/15] Updates to allow varying shape crs --- lib/iris/_shapefiles.py | 88 +++++++++++++++++++++++++++-------------- lib/iris/util.py | 31 ++++++++++++--- 2 files changed, 84 insertions(+), 35 deletions(-) diff --git a/lib/iris/_shapefiles.py b/lib/iris/_shapefiles.py index 74b24b6627..45d3c1a0fb 100644 --- a/lib/iris/_shapefiles.py +++ b/lib/iris/_shapefiles.py @@ -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 @@ -43,6 +45,11 @@ 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 ------- @@ -50,6 +57,17 @@ def create_shapefile_mask( 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 @@ -116,7 +134,11 @@ 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 @@ -125,10 +147,11 @@ 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 ------- @@ -136,25 +159,40 @@ def _transform_coord_system(geometry, cube, geometry_system=None): 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) @@ -165,16 +203,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 diff --git a/lib/iris/util.py b/lib/iris/util.py index 4924ca68d2..332c354652 100644 --- a/lib/iris/util.py +++ b/lib/iris/util.py @@ -2122,9 +2122,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 @@ -2133,10 +2141,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. @@ -2156,8 +2166,7 @@ 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 @@ -2174,8 +2183,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 From e32f911b1e08f94e8d2afc6f38a5806ccf2d0d81 Mon Sep 17 00:00:00 2001 From: Hamish Steptoe Date: Thu, 22 Aug 2024 15:06:34 +0100 Subject: [PATCH 02/15] Ruff fixes --- lib/iris/_shapefiles.py | 4 ++-- lib/iris/util.py | 1 - 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/lib/iris/_shapefiles.py b/lib/iris/_shapefiles.py index 45d3c1a0fb..ab07e4e727 100644 --- a/lib/iris/_shapefiles.py +++ b/lib/iris/_shapefiles.py @@ -135,8 +135,8 @@ def create_shapefile_mask( def _transform_coord_system( - geometry: shapely.Geometry, - cube: iris.cube.Cube, + geometry: shapely.Geometry, + cube: iris.cube.Cube, geometry_crs: cartopy.crs = None ) -> shapely.Geometry: """Project the shape onto another coordinate system. diff --git a/lib/iris/util.py b/lib/iris/util.py index 332c354652..6d1c0bcf40 100644 --- a/lib/iris/util.py +++ b/lib/iris/util.py @@ -2166,7 +2166,6 @@ def mask_cube_from_shapefile( Notes ----- - 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 From 76eae38439c4b9a8a5d1e63a9ea4cdbe737248e5 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 22 Aug 2024 14:36:40 +0000 Subject: [PATCH 03/15] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- lib/iris/_shapefiles.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/lib/iris/_shapefiles.py b/lib/iris/_shapefiles.py index ab07e4e727..735686cc7c 100644 --- a/lib/iris/_shapefiles.py +++ b/lib/iris/_shapefiles.py @@ -135,9 +135,7 @@ def create_shapefile_mask( def _transform_coord_system( - geometry: shapely.Geometry, - cube: iris.cube.Cube, - geometry_crs: cartopy.crs = None + geometry: shapely.Geometry, cube: iris.cube.Cube, geometry_crs: cartopy.crs = None ) -> shapely.Geometry: """Project the shape onto another coordinate system. @@ -159,7 +157,7 @@ def _transform_coord_system( A transformed copy of the provided :class:`shapely.Geometry`. """ - _ , x_name = _cube_primary_xy_coord_names(cube) + _, x_name = _cube_primary_xy_coord_names(cube) target_system = cube.coord_system().as_cartopy_projection() if not target_system: From ebaa1ee8d7a095ca9b7d1ac0cb19a3a2ed129524 Mon Sep 17 00:00:00 2001 From: Hamish Steptoe Date: Fri, 6 Dec 2024 14:06:25 +0000 Subject: [PATCH 04/15] Adding rasterio as optional dependency --- pyproject.toml | 3 +++ requirements/py312.yml | 1 + 2 files changed, 4 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 74e514ad20..e1715f04fa 100644 --- a/pyproject.toml +++ b/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" diff --git a/requirements/py312.yml b/requirements/py312.yml index 74277e417f..3882bba3d4 100644 --- a/requirements/py312.yml +++ b/requirements/py312.yml @@ -34,6 +34,7 @@ dependencies: - pandas - pip - python-stratify + - rasterio # Test dependencies. - asv_runner From 103e3d7bee970d25ca916782a4b943964620405c Mon Sep 17 00:00:00 2001 From: Hamish Steptoe Date: Thu, 19 Dec 2024 16:54:13 +0000 Subject: [PATCH 05/15] Adding rasterio functionality --- lib/iris/_shapefiles.py | 245 ++++++++++++++++------------------------ 1 file changed, 100 insertions(+), 145 deletions(-) diff --git a/lib/iris/_shapefiles.py b/lib/iris/_shapefiles.py index 735686cc7c..2fd6fcc602 100644 --- a/lib/iris/_shapefiles.py +++ b/lib/iris/_shapefiles.py @@ -7,31 +7,48 @@ # the Met Office by Chris Kent, Emilie Vanvyve, David Bentley, Joana Mendes # many thanks to them. Converted to iris by Alex Chamberlain-Clay +from __future__ import annotations +import importlib from itertools import product import warnings +import sys import numpy as np import shapely import shapely.errors import shapely.geometry as sgeom import shapely.ops +from pyproj import CRS -import iris.analysis.cartography +import rasterio.features as rfeatures +import rasterio.transform as rtransform +import rasterio.warp as rwarp + +import iris from iris.warnings import IrisDefaultingWarning, IrisUserWarning +import pdb + +if "iris.cube" in sys.modules: + import iris.cube +if "iris.analysis.cartography" in sys.modules: + import iris.analysis.cartography + def create_shapefile_mask( geometry: shapely.Geometry, + geometry_crs: cartopy.crs | CRS, cube: iris.cube.Cube, - minimum_weight: float = 0.0, - geometry_crs: cartopy.crs = None, + **kwargs, ) -> np.array: """Make a mask for a cube from a shape. Get the mask of the intersection between the given shapely geometry and cube with x/y DimCoords. - Can take a minimum weight and evaluate area overlaps instead + Can take a minimum weight and evaluate area overlaps instead. + + Transforming is performed by GDAL warp. Parameters ---------- @@ -50,6 +67,11 @@ def create_shapefile_mask( the coord_system of the shapefile. Defaults to None, in which case the geometry_crs is assumed to be the same as the `cube`. + **kwargs + Additional keyword arguments to pass to `rasterio.features.geometry_mask`. + Valid keyword arguments are: + all_touched : bool, optional + invert : bool, optional Returns ------- @@ -57,6 +79,13 @@ def create_shapefile_mask( An array of the shape of the x & y coordinates of the cube, with points to mask equal to True. + Raises + ------ + ValueError + If the geometry is not valid for the given coordinate system. + This most likely occurs when the geometry coordinates are not within the bounds of the + geometry coordinates reference system. + Notes ----- For best masking results, both the cube _and_ masking geometry should have a @@ -69,150 +98,100 @@ def create_shapefile_mask( If a CRS is not provided for the the masking geometry, the CRS of the cube is assumed. """ - from iris.cube import Cube, CubeList - - try: + if not shapely.is_valid(geometry): msg = "Geometry is not a valid Shapely object" - if not shapely.is_valid(geometry): - raise TypeError(msg) - except Exception: raise TypeError(msg) - if not isinstance(cube, Cube): - if isinstance(cube, CubeList): + + if not isinstance(cube, iris.cube.Cube): + if isinstance(cube, iris.cube.CubeList): msg = "Received CubeList object rather than Cube - \ to mask a CubeList iterate over each Cube" raise TypeError(msg) else: msg = "Received non-Cube object where a Cube is expected" raise TypeError(msg) - if minimum_weight > 0.0 and isinstance( - geometry, - ( - sgeom.Point, - sgeom.LineString, - sgeom.LinearRing, - sgeom.MultiPoint, - sgeom.MultiLineString, - ), - ): - minimum_weight = 0.0 - warnings.warn( - """Shape is of invalid type for minimum weight masking, - must use a Polygon rather than Line shape.\n - Masking based off intersection instead. """, - category=IrisDefaultingWarning, + + # Check validity of geometry + geom_valid = _is_geometry_valid(geometry, geometry_crs) + if not geom_valid.all(): + raise ValueError( + f"Geometry {shapely.get_parts(geometry)[~geom_valid]} is not valid for the given coordinate system {geometry_crs.to_string()}. Check that your coordinates are correctly specified." ) - # prepare 2D cube + # Define raster transform based on cube + dst_crs = cube.coord_system().as_cartopy_projection() y_name, x_name = _cube_primary_xy_coord_names(cube) - cube_2d = cube.slices([y_name, x_name]).next() - for coord in cube_2d.dim_coords: - if not coord.has_bounds(): - coord.guess_bounds() - trans_geo = _transform_coord_system(geometry, cube_2d) - - y_coord, x_coord = [cube_2d.coord(n) for n in (y_name, x_name)] - x_bounds = _get_mod_rebased_coord_bounds(x_coord) - y_bounds = _get_mod_rebased_coord_bounds(y_coord) - # prepare array for dark - box_template = [ - sgeom.box(x[0], y[0], x[1], y[1]) for x, y in product(x_bounds, y_bounds) - ] - # shapely can do lazy evaluation of intersections if it's given a list of grid box shapes - # delayed lets us do it in parallel - intersect_template = shapely.intersects(trans_geo, box_template) - # we want areas not under shapefile to be True (to mask) - intersect_template = np.invert(intersect_template) - # now calc area overlaps if doing weights and adjust mask - if minimum_weight > 0.0: - intersections = np.array(box_template)[~intersect_template] - intersect_template[~intersect_template] = [ - trans_geo.intersection(box).area / box.area <= minimum_weight - for box in intersections - ] - mask_template = np.reshape(intersect_template, cube_2d.shape[::-1]).T + w = cube.coord(x_name).points.size + h = cube.coord(y_name).points.size + + tr, w, h = rwarp.calculate_default_transform( + src_crs=geometry_crs, + dst_crs=dst_crs, + width=w, + height=h, + dst_width=w, + dst_height=h, + src_geoloc_array=( + np.meshgrid( + cube.coord(x_name).points, cube.coord(y_name).points, indexing="xy" + ) + ), + ) + + # Check for CRS equality and transform if necessary + if not geometry_crs.equals(dst_crs): + trans_geometry = rwarp.transform_geom( + geom=geometry, src_crs=geometry_crs, dst_crs=dst_crs + ) + geometry = sgeom.shape(trans_geometry) + + mask_template = rfeatures.geometry_mask( + geometries=shapely.get_parts(geometry), out_shape=(h, w), transform=tr + ) + return mask_template -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. +def _is_geometry_valid( + geometry: shapely.Geometry, + geometry_crs: cartopy.crs, +) -> np.array: + """Check if a shape geometry is valid given its coordinate system. + + Achieved by asserting that the geometry, falls within bounds equivalent to + lon = [-180, 180] and lat = [-90, 90]. Parameters ---------- - geometry : :class:`shapely.Geometry` - cube : :class:`iris.cube.Cube` - :class:`~iris.cube.Cube` with the coord_system to be projected to and - a x coordinate. - geometry_crs : :class:`cartopy.crs`, optional - A :class:`cartopy.crs` 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`. + geometry : :class:`shapely.geometry.base.BaseGeometry` + The geometry to check. + geometry_crs : :class:`cartopy.crs` + The coordinate reference system of the geometry. Returns ------- - :class:`shapely.Geometry` - A transformed copy of the provided :class:`shapely.Geometry`. + np.array of bool + True if the geometry is valid, False otherwise for each part of geometry. """ - _, x_name = _cube_primary_xy_coord_names(cube) - - 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( - "No geometry coordinate reference system supplied; using cube coord_system instead.", - category=IrisDefaultingWarning, - ) - - trans_geometry = target_proj.project_geometry(geometry, source_proj) + WGS84_crs = CRS.from_epsg(4326) - # 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) - and cube.coord(x_name).points[-1] > 180 - ): - # chop geom at 0 degree line very finely then transform - prime_meridian_line = shapely.LineString([(0, 90), (0, -90)]) - trans_geometry = trans_geometry.difference(prime_meridian_line.buffer(0.00001)) - trans_geometry = shapely.transform(trans_geometry, _trans_func) - - return trans_geometry + lon_lat_bounds = shapely.geometry.Polygon.from_bounds( + xmin=-180.0, ymin=-90.0, xmax=180.0, ymax=90.0 + ) + # If the geometry is not in WGS84, transform the validation bounds + # to the geometry's CRS + if not geometry_crs.equals(WGS84_crs): + lon_lat_bounds = rwarp.transform_geom( + src_crs=WGS84_crs, dst_crs=geometry_crs, geom=lon_lat_bounds + ) + lon_lat_bounds = sgeom.shape(lon_lat_bounds) -def _trans_func(geometry): - """Pocket function for transforming the x coord of a geometry from -180 to 180 to 0-360.""" - for point in geometry: - if point[0] < 0: - point[0] = 360 - np.abs(point[0]) - return geometry + return lon_lat_bounds.contains(shapely.get_parts(geometry)) -def _cube_primary_xy_coord_names(cube): +def _cube_primary_xy_coord_names(cube: iris.cube.Cube) -> tuple[str, str]: """Return the primary latitude and longitude coordinate names, or long names, from a cube. Parameters @@ -243,27 +222,3 @@ def _cube_primary_xy_coord_names(cube): latitude = latc.name() longitude = lonc.name() return latitude, longitude - - -def _get_mod_rebased_coord_bounds(coord): - """Take in a coord and returns a array of the bounds of that coord rebased to the modulus. - - Parameters - ---------- - coord : :class:`iris.coords.Coord` - An Iris coordinate with a modulus. - - Returns - ------- - :class:`np.array` - A 1d Numpy array of [start,end] pairs for bounds of the coord. - - """ - modulus = coord.units.modulus - # Force realisation (rather than core_bounds) - more efficient for the - # repeated indexing happening downstream. - result = np.array(coord.bounds) - if modulus: - result[result < 0.0] = (np.abs(result[result < 0.0]) % modulus) * -1 - result[np.isclose(result, modulus, 1e-10)] = 0.0 - return result From 3c65dab26fe91e12fc16ca03b3f57bccacabe828 Mon Sep 17 00:00:00 2001 From: Hamish Steptoe Date: Fri, 20 Dec 2024 15:55:38 +0000 Subject: [PATCH 06/15] First working version with rasterio --- lib/iris/_shapefiles.py | 61 +++++++++++++++++++++++++++++++---------- 1 file changed, 47 insertions(+), 14 deletions(-) diff --git a/lib/iris/_shapefiles.py b/lib/iris/_shapefiles.py index 2fd6fcc602..af98b0d7ce 100644 --- a/lib/iris/_shapefiles.py +++ b/lib/iris/_shapefiles.py @@ -118,14 +118,28 @@ def create_shapefile_mask( f"Geometry {shapely.get_parts(geometry)[~geom_valid]} is not valid for the given coordinate system {geometry_crs.to_string()}. Check that your coordinates are correctly specified." ) - # Define raster transform based on cube + # Check for CRS equality and transform if necessary dst_crs = cube.coord_system().as_cartopy_projection() + if not geometry_crs.equals(dst_crs): + trans_geometry = rwarp.transform_geom( + geom=geometry, src_crs=geometry_crs, dst_crs=dst_crs + ) + geometry = sgeom.shape(trans_geometry) + + # Get cube coordinates y_name, x_name = _cube_primary_xy_coord_names(cube) - w = cube.coord(x_name).points.size - h = cube.coord(y_name).points.size + x_points = iris.analysis.cartography.wrap_lons( + cube.coord(x_name).points, base=-180, period=360 + ) + y_points = cube.coord(y_name).points + w = len(x_points) + h = len(y_points) + # Define raster transform based on cube + # This maps the geometry domain onto the cube domain + # using an Affine transformation tr, w, h = rwarp.calculate_default_transform( - src_crs=geometry_crs, + src_crs=dst_crs, dst_crs=dst_crs, width=w, height=h, @@ -133,23 +147,18 @@ def create_shapefile_mask( dst_height=h, src_geoloc_array=( np.meshgrid( - cube.coord(x_name).points, cube.coord(y_name).points, indexing="xy" + x_points, + y_points, + indexing="xy", ) ), ) - - # Check for CRS equality and transform if necessary - if not geometry_crs.equals(dst_crs): - trans_geometry = rwarp.transform_geom( - geom=geometry, src_crs=geometry_crs, dst_crs=dst_crs - ) - geometry = sgeom.shape(trans_geometry) - + # Generate mask from geometry mask_template = rfeatures.geometry_mask( geometries=shapely.get_parts(geometry), out_shape=(h, w), transform=tr ) - return mask_template + return mask_template[::-1,:] #TODO: check why need to reverse i dimension def _is_geometry_valid( @@ -222,3 +231,27 @@ def _cube_primary_xy_coord_names(cube: iris.cube.Cube) -> tuple[str, str]: latitude = latc.name() longitude = lonc.name() return latitude, longitude + + +def _get_mod_rebased_coord_bounds(coord): + """Take in a coord and returns a array of the bounds of that coord rebased to the modulus. + + Parameters + ---------- + coord : :class:`iris.coords.Coord` + An Iris coordinate with a modulus. + + Returns + ------- + :class:`np.array` + A 1d Numpy array of [start,end] pairs for bounds of the coord. + + """ + modulus = coord.units.modulus + # Force realisation (rather than core_bounds) - more efficient for the + # repeated indexing happening downstream. + result = np.array(coord.bounds) + if modulus: + result[result < 0.0] = (np.abs(result[result < 0.0]) % modulus) * -1 + result[np.isclose(result, modulus, 1e-10)] = 0.0 + return result From 25de18dfe0192f046bbb548d9925bfcb5a16d4b9 Mon Sep 17 00:00:00 2001 From: Hamish Steptoe Date: Fri, 3 Jan 2025 11:46:22 +0000 Subject: [PATCH 07/15] Add further shape geometry checks --- lib/iris/_shapefiles.py | 166 +++++++++++++++++++++++++--------------- 1 file changed, 104 insertions(+), 62 deletions(-) diff --git a/lib/iris/_shapefiles.py b/lib/iris/_shapefiles.py index af98b0d7ce..fe14d2a068 100644 --- a/lib/iris/_shapefiles.py +++ b/lib/iris/_shapefiles.py @@ -11,25 +11,22 @@ import importlib from itertools import product -import warnings import sys +import warnings import numpy as np -import shapely -import shapely.errors -import shapely.geometry as sgeom -import shapely.ops from pyproj import CRS - import rasterio.features as rfeatures import rasterio.transform as rtransform import rasterio.warp as rwarp +import shapely +import shapely.errors +import shapely.geometry as sgeom +import shapely.ops import iris from iris.warnings import IrisDefaultingWarning, IrisUserWarning -import pdb - if "iris.cube" in sys.modules: import iris.cube if "iris.analysis.cartography" in sys.modules: @@ -81,10 +78,8 @@ def create_shapefile_mask( Raises ------ - ValueError - If the geometry is not valid for the given coordinate system. - This most likely occurs when the geometry coordinates are not within the bounds of the - geometry coordinates reference system. + TypeError + If the cube is not a :class:`~iris.cube.Cube`. Notes ----- @@ -97,11 +92,16 @@ def create_shapefile_mask( 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. + + Warning + ------- + Using a shapefile that crosses the 180th meridian may lead to unexpected masking behaviour. + For best results, ensure that the shapefile is pre-split at the 180th meridian. """ - if not shapely.is_valid(geometry): - msg = "Geometry is not a valid Shapely object" - raise TypeError(msg) + # Check validity of geometry CRS + is_geometry_valid(geometry, geometry_crs) + # Check cube is a Cube if not isinstance(cube, iris.cube.Cube): if isinstance(cube, iris.cube.CubeList): msg = "Received CubeList object rather than Cube - \ @@ -111,13 +111,6 @@ def create_shapefile_mask( msg = "Received non-Cube object where a Cube is expected" raise TypeError(msg) - # Check validity of geometry - geom_valid = _is_geometry_valid(geometry, geometry_crs) - if not geom_valid.all(): - raise ValueError( - f"Geometry {shapely.get_parts(geometry)[~geom_valid]} is not valid for the given coordinate system {geometry_crs.to_string()}. Check that your coordinates are correctly specified." - ) - # Check for CRS equality and transform if necessary dst_crs = cube.coord_system().as_cartopy_projection() if not geometry_crs.equals(dst_crs): @@ -128,9 +121,10 @@ def create_shapefile_mask( # Get cube coordinates y_name, x_name = _cube_primary_xy_coord_names(cube) - x_points = iris.analysis.cartography.wrap_lons( - cube.coord(x_name).points, base=-180, period=360 - ) + # Check if cube lons exist in [0, 360] or [-180, 180] + if cube.coord(x_name).circular: + cube = cube.intersection(longitude=(-180, 180)) + x_points = cube.coord(x_name).points y_points = cube.coord(y_name).points w = len(x_points) h = len(y_points) @@ -158,17 +152,24 @@ def create_shapefile_mask( geometries=shapely.get_parts(geometry), out_shape=(h, w), transform=tr ) - return mask_template[::-1,:] #TODO: check why need to reverse i dimension + return mask_template -def _is_geometry_valid( +def is_geometry_valid( geometry: shapely.Geometry, geometry_crs: cartopy.crs, -) -> np.array: - """Check if a shape geometry is valid given its coordinate system. - - Achieved by asserting that the geometry, falls within bounds equivalent to - lon = [-180, 180] and lat = [-90, 90]. +) -> None: + """Check the validity of the shape geometry. + + This function checks that: + 1) The geometry is a valid shapely geometry. + 2) The geometry falls within bounds equivalent to + lon = [-180, 180] and lat = [-90, 90]. + 3) The geometry does not cross the 180th meridian, + based on the assumption that the shape will + cross the 180th meridian if the difference between + the shape bound longitudes is greater than 180. + 4) The geometry does not cross the poles. Parameters ---------- @@ -179,25 +180,90 @@ def _is_geometry_valid( Returns ------- - np.array of bool - True if the geometry is valid, False otherwise for each part of geometry. + None if the geometry is valid. + Raises + ------ + TypeError + If the geometry is not a valid shapely geometry. + ValueError + If the geometry is not valid for the given coordinate system. + This most likely occurs when the geometry coordinates are not within the bounds of the + geometry coordinates reference system. + ValueError + If the geometry crosses the 180th meridian. + ValueError + If the geometry crosses the poles. + + Examples + -------- + >>> from shapely.geometry import box + >>> from pyproj import CRS + >>> from iris._shapefiles import is_geometry_valid + + Create a valid geometry covering Canada, and check + its validity for the WGS84 coordinate system: + + >>> canada = box(-143.5,42.6,-37.8,84.0) + >>> wgs84 = CRS.from_epsg(4326) + >>> is_geometry_valid(canada, wgs84) + + The function returns silently as the geometry is valid. + + Now create an invalid geometry covering the Bering Sea, + and check its validity for the WGS84 coordinate system. + + >>> bering_sea = box(148.42,49.1,-138.74,73.12) + >>> wgs84 = CRS.from_epsg(4326) + >>> is_geometry_valid(bering_sea, wgs84) # doctest: +IGNORE_EXCEPTION_DETAIL + Traceback (most recent call last) + ValueError: Geometry crossing the 180th meridian is not supported. """ WGS84_crs = CRS.from_epsg(4326) + # Check geometry is valid shapely geometry + if not shapely.is_valid(geometry): + msg = "Geometry is not a valid Shapely object" + raise TypeError(msg) + + # Check that the geometry is within the bounds of the coordinate system + # If the geometry is not in WGS84, transform the validation bounds + # to the geometry's CR lon_lat_bounds = shapely.geometry.Polygon.from_bounds( xmin=-180.0, ymin=-90.0, xmax=180.0, ymax=90.0 ) - - # If the geometry is not in WGS84, transform the validation bounds - # to the geometry's CRS if not geometry_crs.equals(WGS84_crs): lon_lat_bounds = rwarp.transform_geom( src_crs=WGS84_crs, dst_crs=geometry_crs, geom=lon_lat_bounds ) lon_lat_bounds = sgeom.shape(lon_lat_bounds) + if not lon_lat_bounds.contains(shapely.get_parts(geometry)).all(): + msg = f"Geometry {shapely.get_parts(geometry)[~geom_valid]} is not valid for the given coordinate system \ + {geometry_crs.to_string()}. Check that your coordinates are correctly specified." + raise ValueError(msg) + + # Check if shape crosses the 180th meridian (or equivalent) + if bool(abs(geometry.bounds[2] - geometry.bounds[0]) > 180.0): + msg = "Geometry crossing the 180th meridian is not supported." + raise ValueError(msg) + + # Check if the geometry crosses the poles + npole = sgeom.Point(0, 90) + spole = sgeom.Point(0, -90) + if not geometry_crs.equals(WGS84_crs): + npole = rwarp.transform_geom( + src_crs=WGS84_crs, dst_crs=geometry_crs, geom=npole + ) + spole = rwarp.transform_geom( + src_crs=WGS84_crs, dst_crs=geometry_crs, geom=spole + ) + npole = sgeom.shape(npole) + spole = sgeom.shape(spole) + if geometry.intersects(npole) or geometry.intersects(spole): + msg = "Geometry crossing the poles is not supported." + raise ValueError(msg) - return lon_lat_bounds.contains(shapely.get_parts(geometry)) + return def _cube_primary_xy_coord_names(cube: iris.cube.Cube) -> tuple[str, str]: @@ -231,27 +297,3 @@ def _cube_primary_xy_coord_names(cube: iris.cube.Cube) -> tuple[str, str]: latitude = latc.name() longitude = lonc.name() return latitude, longitude - - -def _get_mod_rebased_coord_bounds(coord): - """Take in a coord and returns a array of the bounds of that coord rebased to the modulus. - - Parameters - ---------- - coord : :class:`iris.coords.Coord` - An Iris coordinate with a modulus. - - Returns - ------- - :class:`np.array` - A 1d Numpy array of [start,end] pairs for bounds of the coord. - - """ - modulus = coord.units.modulus - # Force realisation (rather than core_bounds) - more efficient for the - # repeated indexing happening downstream. - result = np.array(coord.bounds) - if modulus: - result[result < 0.0] = (np.abs(result[result < 0.0]) % modulus) * -1 - result[np.isclose(result, modulus, 1e-10)] = 0.0 - return result From 10e4413b1d63b8ef5f1c78ab923d5a2cabcb2f6c Mon Sep 17 00:00:00 2001 From: Hamish Steptoe Date: Fri, 3 Jan 2025 14:34:19 +0000 Subject: [PATCH 08/15] Update shapefile tests --- lib/iris/_shapefiles.py | 26 ++++-- lib/iris/tests/unit/_shapefiles/__init__.py | 5 ++ .../_shapefiles/test_is_geometry_valid.py | 56 +++++++++++++ lib/iris/tests/unit/conftest.py | 83 +++++++++++++++++++ .../util/test_mask_cube_from_shapefile.py | 73 ++++++++++++++++ 5 files changed, 234 insertions(+), 9 deletions(-) create mode 100644 lib/iris/tests/unit/_shapefiles/__init__.py create mode 100644 lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py diff --git a/lib/iris/_shapefiles.py b/lib/iris/_shapefiles.py index fe14d2a068..fc24964bfe 100644 --- a/lib/iris/_shapefiles.py +++ b/lib/iris/_shapefiles.py @@ -95,8 +95,16 @@ def create_shapefile_mask( Warning ------- - Using a shapefile that crosses the 180th meridian may lead to unexpected masking behaviour. - For best results, ensure that the shapefile is pre-split at the 180th meridian. + Because shape vectors are inherently Cartesian in nature, they contain no inherent + understanding of the spherical geometry underpinning geographic coordinate systems. + For this reason, shapefiles or shape vectors that cross the antimeridian or poles + are not supported by this function to avoid unexpected masking behaviour. + + Shape geometries can be checked prior to masking using the :func:`is_geometry_valid`. + + See Also + -------- + :func:`is_geometry_valid` """ # Check validity of geometry CRS is_geometry_valid(geometry, geometry_crs) @@ -165,9 +173,9 @@ def is_geometry_valid( 1) The geometry is a valid shapely geometry. 2) The geometry falls within bounds equivalent to lon = [-180, 180] and lat = [-90, 90]. - 3) The geometry does not cross the 180th meridian, + 3) The geometry does not cross the antimeridian, based on the assumption that the shape will - cross the 180th meridian if the difference between + cross the antimeridian if the difference between the shape bound longitudes is greater than 180. 4) The geometry does not cross the poles. @@ -191,7 +199,7 @@ def is_geometry_valid( This most likely occurs when the geometry coordinates are not within the bounds of the geometry coordinates reference system. ValueError - If the geometry crosses the 180th meridian. + If the geometry crosses the antimeridian. ValueError If the geometry crosses the poles. @@ -237,14 +245,14 @@ def is_geometry_valid( src_crs=WGS84_crs, dst_crs=geometry_crs, geom=lon_lat_bounds ) lon_lat_bounds = sgeom.shape(lon_lat_bounds) - if not lon_lat_bounds.contains(shapely.get_parts(geometry)).all(): - msg = f"Geometry {shapely.get_parts(geometry)[~geom_valid]} is not valid for the given coordinate system \ - {geometry_crs.to_string()}. Check that your coordinates are correctly specified." + geom_valid = lon_lat_bounds.contains(shapely.get_parts(geometry)) + if not geom_valid.all(): + msg = f"Geometry {shapely.get_parts(geometry)[~geom_valid]} is not valid for the given coordinate system {geometry_crs.to_string()}. \nCheck that your coordinates are correctly specified." raise ValueError(msg) # Check if shape crosses the 180th meridian (or equivalent) if bool(abs(geometry.bounds[2] - geometry.bounds[0]) > 180.0): - msg = "Geometry crossing the 180th meridian is not supported." + msg = "Geometry crossing the antimeridian is not supported." raise ValueError(msg) # Check if the geometry crosses the poles diff --git a/lib/iris/tests/unit/_shapefiles/__init__.py b/lib/iris/tests/unit/_shapefiles/__init__.py new file mode 100644 index 0000000000..0d51a683b8 --- /dev/null +++ b/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.""" diff --git a/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py b/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py new file mode 100644 index 0000000000..bcfceaccf7 --- /dev/null +++ b/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py @@ -0,0 +1,56 @@ +# 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`.""" + +from pyproj import CRS +import pytest +from shapely.geometry import box + +from iris._shapefiles import is_geometry_valid + + +# Fixtures retrieved from conftest.py +@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(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 [] 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(test_input, wgs84) diff --git a/lib/iris/tests/unit/conftest.py b/lib/iris/tests/unit/conftest.py index 524ca53ce8..ee6e194c4c 100644 --- a/lib/iris/tests/unit/conftest.py +++ b/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 @@ -13,3 +14,85 @@ def sample_coord(): sample_coord = iris.coords.DimCoord(points=(1, 2, 3, 4, 5)) return sample_coord + + +# 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 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 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) \ No newline at end of file diff --git a/lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py b/lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py index bdd5c5fc56..34ab462df0 100644 --- a/lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py +++ b/lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py @@ -5,8 +5,12 @@ """Unit tests for :func:`iris.util.mask_cube_from_shapefile`.""" import numpy as np +from pyproj import CRS import pytest import shapely +from shapely import Point, MultiPoint, LineString +from shapely.geometry import box +from cartopy.feature import ShapelyFeature from iris.coord_systems import RotatedGeogCS from iris.coords import DimCoord @@ -14,7 +18,76 @@ 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 + +@pytest.fixture +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 +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 + rectangle = box(min_lon, min_lat, max_lon, max_lat) + + return rectangle + + +@pytest.fixture +def basic_point_geometry(): + # Define the coordinates of a basic point (lon, lat) + point = Point((-3.476204, 50.727059)) + + return point + +@pytest.fixture +def basic_line_geometry(): + # Define the coordinates of a basic line + line = LineString([(0, 0), (10, 10)]) + + return line + +@pytest.fixture +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 +def basic_line_geometry() class TestBasicCubeMasking: """Unit tests for mask_cube_from_shapefile function.""" From cd5cf991d72b22af0c500e953e5bd82040eef48b Mon Sep 17 00:00:00 2001 From: Hamish Steptoe Date: Fri, 3 Jan 2025 14:42:01 +0000 Subject: [PATCH 09/15] Reorganise test fixtures --- .../_shapefiles/test_is_geometry_valid.py | 82 ++++++++++++++++++- 1 file changed, 81 insertions(+), 1 deletion(-) diff --git a/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py b/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py index bcfceaccf7..db0591a813 100644 --- a/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py +++ b/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py @@ -11,7 +11,87 @@ from iris._shapefiles import is_geometry_valid -# Fixtures retrieved from conftest.py +# 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 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 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) + @pytest.mark.parametrize( "test_input", [ From edc25cc5bf1d358ce345c6bdc08b41d68fc563f5 Mon Sep 17 00:00:00 2001 From: Hamish Steptoe Date: Fri, 3 Jan 2025 14:43:11 +0000 Subject: [PATCH 10/15] Reorganise test fixtures --- lib/iris/tests/unit/conftest.py | 81 --------------------------------- 1 file changed, 81 deletions(-) diff --git a/lib/iris/tests/unit/conftest.py b/lib/iris/tests/unit/conftest.py index ee6e194c4c..acfc074684 100644 --- a/lib/iris/tests/unit/conftest.py +++ b/lib/iris/tests/unit/conftest.py @@ -15,84 +15,3 @@ def sample_coord(): sample_coord = iris.coords.DimCoord(points=(1, 2, 3, 4, 5)) return sample_coord - -# 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 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 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) \ No newline at end of file From 789c98d21aecdd31b6f9b6fdfac5778a655a2237 Mon Sep 17 00:00:00 2001 From: Hamish Steptoe Date: Fri, 3 Jan 2025 15:03:03 +0000 Subject: [PATCH 11/15] Fix fixture fetching --- .../_shapefiles/test_is_geometry_valid.py | 107 +++++++++++------- .../util/test_mask_cube_from_shapefile.py | 58 +--------- 2 files changed, 66 insertions(+), 99 deletions(-) diff --git a/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py b/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py index db0591a813..7cadae9e7e 100644 --- a/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py +++ b/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py @@ -4,6 +4,10 @@ # 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 @@ -18,6 +22,12 @@ 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 @@ -82,55 +92,68 @@ 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) -@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(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 [] 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(test_input, wgs84) + +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 [] 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) diff --git a/lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py b/lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py index 34ab462df0..3d19d8c9a4 100644 --- a/lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py +++ b/lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py @@ -30,64 +30,8 @@ # geometry["crs"] = crs # return geometry -@pytest.fixture -def basic_circular_geometry(): - # Define geometry of a basic circle with center at (0, 0) and radius 10 - center = (0, 0) - radius = 10 +# Fixtures retrieved from _shapefiles/test_is_geometry_valid.py - circle = Point(center).buffer(radius) - - return circle - - -@pytest.fixture -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 - rectangle = box(min_lon, min_lat, max_lon, max_lat) - - return rectangle - - -@pytest.fixture -def basic_point_geometry(): - # Define the coordinates of a basic point (lon, lat) - point = Point((-3.476204, 50.727059)) - - return point - -@pytest.fixture -def basic_line_geometry(): - # Define the coordinates of a basic line - line = LineString([(0, 0), (10, 10)]) - - return line - -@pytest.fixture -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 -def basic_line_geometry() class TestBasicCubeMasking: """Unit tests for mask_cube_from_shapefile function.""" From dc090a5028e9ef403f3e2a00aecc2a356901e1d7 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 3 Jan 2025 15:04:26 +0000 Subject: [PATCH 12/15] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py | 1 + lib/iris/tests/unit/conftest.py | 1 - lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py | 5 +++-- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py b/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py index 7cadae9e7e..bb54de15cf 100644 --- a/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py +++ b/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py @@ -92,6 +92,7 @@ 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 diff --git a/lib/iris/tests/unit/conftest.py b/lib/iris/tests/unit/conftest.py index acfc074684..219b23b1d1 100644 --- a/lib/iris/tests/unit/conftest.py +++ b/lib/iris/tests/unit/conftest.py @@ -14,4 +14,3 @@ def sample_coord(): sample_coord = iris.coords.DimCoord(points=(1, 2, 3, 4, 5)) return sample_coord - diff --git a/lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py b/lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py index 3d19d8c9a4..2cb498a312 100644 --- a/lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py +++ b/lib/iris/tests/unit/util/test_mask_cube_from_shapefile.py @@ -4,13 +4,13 @@ # 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 Point, MultiPoint, LineString +from shapely import LineString, MultiPoint, Point from shapely.geometry import box -from cartopy.feature import ShapelyFeature from iris.coord_systems import RotatedGeogCS from iris.coords import DimCoord @@ -32,6 +32,7 @@ # Fixtures retrieved from _shapefiles/test_is_geometry_valid.py + class TestBasicCubeMasking: """Unit tests for mask_cube_from_shapefile function.""" From 8948ce6ce36f914b6f0660186cb1fe31a3a8e299 Mon Sep 17 00:00:00 2001 From: Hamish Steptoe Date: Thu, 13 Feb 2025 16:39:11 +0000 Subject: [PATCH 13/15] Recording Affine transform trials --- lib/iris/_shapefiles.py | 166 +++++++++++++++++++++++++--------------- 1 file changed, 106 insertions(+), 60 deletions(-) diff --git a/lib/iris/_shapefiles.py b/lib/iris/_shapefiles.py index fc24964bfe..06a444d1b1 100644 --- a/lib/iris/_shapefiles.py +++ b/lib/iris/_shapefiles.py @@ -9,20 +9,22 @@ from __future__ import annotations +import pdb +from typing import overload import importlib from itertools import product import sys import warnings import numpy as np -from pyproj import CRS +from pyproj import CRS, Transformer import rasterio.features as rfeatures import rasterio.transform as rtransform import rasterio.warp as rwarp import shapely import shapely.errors import shapely.geometry as sgeom -import shapely.ops +import shapely.ops as sops import iris from iris.warnings import IrisDefaultingWarning, IrisUserWarning @@ -33,6 +35,14 @@ import iris.analysis.cartography +@overload +def create_shapefile_mask( + geometry: shapely.Geometry, + geometry_crs: cartopy.crs | CRS, + cube: iris.cube.Cube, + all_touched: bool = False, + invert: bool = False, +) -> np.array: ... def create_shapefile_mask( geometry: shapely.Geometry, geometry_crs: cartopy.crs | CRS, @@ -50,21 +60,14 @@ def create_shapefile_mask( Parameters ---------- geometry : :class:`shapely.Geometry` - cube : :class:`iris.cube.Cube` - A :class:`~iris.cube.Cube` which has 1d x and y coordinates. - minimum_weight : float, default 0.0 - A float between 0 and 1 determining what % of a cell - a shape must cover for the cell to remain unmasked. - eg: 0.1 means that at least 10% of the shape overlaps the cell - 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`. - **kwargs + same as the :class:`iris.cube.Cube`. + cube : :class:`iris.cube.Cube` + A :class:`~iris.cube.Cube` which has 1d x and y coordinates. + **kwargs : Additional keyword arguments to pass to `rasterio.features.geometry_mask`. Valid keyword arguments are: all_touched : bool, optional @@ -83,15 +86,15 @@ def create_shapefile_mask( Notes ----- - For best masking results, both the cube _and_ masking geometry should have a + For best masking results, both the :class:`iris.cube.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. + when the :class:`iris.cube.Cube` and masking geometry have the same CRS. - If the cube has no coord_system, the default GeogCS is used where + If the :class:`iris.cube.Cube` has no :class:`~iris.coord_systems`, the default GeogCS is used where the coordinate units are degrees. For any other coordinate units, - the cube **must** have a coord_system defined. + the cube **must** have a :class:`~iris.coord_systems` defined. - If a CRS is not provided for the the masking geometry, the CRS of the cube is assumed. + If a CRS is not provided for the the masking geometry, the CRS of the :class:`iris.cube.Cube` is assumed. Warning ------- @@ -120,18 +123,31 @@ def create_shapefile_mask( raise TypeError(msg) # Check for CRS equality and transform if necessary - dst_crs = cube.coord_system().as_cartopy_projection() - if not geometry_crs.equals(dst_crs): - trans_geometry = rwarp.transform_geom( - geom=geometry, src_crs=geometry_crs, dst_crs=dst_crs - ) - geometry = sgeom.shape(trans_geometry) + cube_crs = cube.coord_system().as_cartopy_projection() + if not geometry_crs.equals(cube_crs): + transform_warning_msg = "Geometry CRS does not match cube CRS. Iris will attempt to transform the geometry onto the cube CRS..." + warnings.warn(transform_warning_msg, category=iris.warnings.IrisUserWarning) + # Set-up transform via pyproj + t = Transformer.from_crs( + crs_from=geometry_crs, crs_to=cube_crs, always_xy=True + ).transform + # Transform geometry + geometry = shapely.ops.transform(t, geometry) + # Recheck geometry validity + if not shapely.is_valid(geometry): + msg = f"Shape geometry is invalid (not well formed): {shapely.is_valid_reason(geometry)}." + raise TypeError(msg) # Get cube coordinates + isLonFix = False y_name, x_name = _cube_primary_xy_coord_names(cube) - # Check if cube lons exist in [0, 360] or [-180, 180] - if cube.coord(x_name).circular: - cube = cube.intersection(longitude=(-180, 180)) + # Check if cube lons units are in degrees, and if so do they exist in [0, 360] or [-180, 180] + if (cube.coord(x_name).units.origin == "degrees") and ( + cube.coord(x_name).points.max() > 180 + ): + isLonFix = True + cube = cube.intersection(iris.coords.CoordExtent(x_name, -180, 180)) + x_points = cube.coord(x_name).points y_points = cube.coord(y_name).points w = len(x_points) @@ -140,26 +156,48 @@ def create_shapefile_mask( # Define raster transform based on cube # This maps the geometry domain onto the cube domain # using an Affine transformation - tr, w, h = rwarp.calculate_default_transform( - src_crs=dst_crs, - dst_crs=dst_crs, - width=w, - height=h, - dst_width=w, - dst_height=h, - src_geoloc_array=( - np.meshgrid( - x_points, - y_points, - indexing="xy", - ) - ), - ) + # tr, w, h = rwarp.calculate_default_transform( + # src_crs=cube_crs, + # dst_crs=cube_crs, + # width=w, + # height=h, + # dst_width=w, + # dst_height=h, + # src_geoloc_array=( + # np.meshgrid( + # x_points, + # y_points, + # indexing="ij", + # ) + # ), + # ) + + # tr = rtransform.from_bounds( + # west=x_points.min(), + # south=y_points.min(), + # east=x_points.max(), + # north=y_points.max(), + # width=h, + # height=-w, + # ) + + tr = _transform_from_latlon(x_points, y_points) + # Generate mask from geometry mask_template = rfeatures.geometry_mask( - geometries=shapely.get_parts(geometry), out_shape=(h, w), transform=tr + geometries=shapely.get_parts(geometry), + out_shape=(h, w), + transform=tr, + **kwargs, ) + # If cube was on [0, 360] domain, then shift mask template + # to match the cube domain + # if isLonFix: + # mask_template = np.roll(mask_template, w // 2, axis=1) + + # pdb.set_trace() + # return mask_template[::-1,:] return mask_template @@ -188,7 +226,7 @@ def is_geometry_valid( Returns ------- - None if the geometry is valid. + None if the geometry is valid.OSTN15_NTv2OSGBtoETRS.gsb Raises ------ @@ -231,20 +269,21 @@ def is_geometry_valid( # Check geometry is valid shapely geometry if not shapely.is_valid(geometry): - msg = "Geometry is not a valid Shapely object" + msg = f"Shape geometry is invalid (not well formed): {shapely.is_valid_reason(geometry)}." raise TypeError(msg) # Check that the geometry is within the bounds of the coordinate system - # If the geometry is not in WGS84, transform the validation bounds - # to the geometry's CR + # If the geometry is not in WGS84, transform the geometry to WGS84 + # This is more reliable than transforming the lon_lat_bounds to the geometry CRS lon_lat_bounds = shapely.geometry.Polygon.from_bounds( xmin=-180.0, ymin=-90.0, xmax=180.0, ymax=90.0 ) if not geometry_crs.equals(WGS84_crs): - lon_lat_bounds = rwarp.transform_geom( - src_crs=WGS84_crs, dst_crs=geometry_crs, geom=lon_lat_bounds - ) - lon_lat_bounds = sgeom.shape(lon_lat_bounds) + # Make pyproj transformer + # Transforms the input geometry to the WGS84 coordinate system + t = Transformer.from_crs(geometry_crs, WGS84_crs, always_xy=True).transform + geometry = shapely.ops.transform(t, geometry) + geom_valid = lon_lat_bounds.contains(shapely.get_parts(geometry)) if not geom_valid.all(): msg = f"Geometry {shapely.get_parts(geometry)[~geom_valid]} is not valid for the given coordinate system {geometry_crs.to_string()}. \nCheck that your coordinates are correctly specified." @@ -258,15 +297,6 @@ def is_geometry_valid( # Check if the geometry crosses the poles npole = sgeom.Point(0, 90) spole = sgeom.Point(0, -90) - if not geometry_crs.equals(WGS84_crs): - npole = rwarp.transform_geom( - src_crs=WGS84_crs, dst_crs=geometry_crs, geom=npole - ) - spole = rwarp.transform_geom( - src_crs=WGS84_crs, dst_crs=geometry_crs, geom=spole - ) - npole = sgeom.shape(npole) - spole = sgeom.shape(spole) if geometry.intersects(npole) or geometry.intersects(spole): msg = "Geometry crossing the poles is not supported." raise ValueError(msg) @@ -305,3 +335,19 @@ def _cube_primary_xy_coord_names(cube: iris.cube.Cube) -> tuple[str, str]: latitude = latc.name() longitude = lonc.name() return latitude, longitude + + +def _transform_from_latlon(lon, lat): + """perform an affine transformation to the latitude/longitude coordinates""" + + from affine import Affine + + lat = np.asarray(lat) + lon = np.asarray(lon) + + d_lon = lon[1] - lon[0] + d_lat = lat[1] - lat[0] + + trans = Affine.translation(lon[0] - d_lon / 2, lat[0] - d_lat / 2) + scale = Affine.scale(d_lon, d_lat) + return trans * scale From a77c6282850fc513011a7cd21a788322d2372619 Mon Sep 17 00:00:00 2001 From: Hamish Steptoe Date: Fri, 14 Feb 2025 16:23:21 +0000 Subject: [PATCH 14/15] Working masking --- lib/iris/_shapefiles.py | 95 +++++++++++++---------------------------- lib/iris/util.py | 76 ++++++++++++++++++++++++++------- 2 files changed, 91 insertions(+), 80 deletions(-) diff --git a/lib/iris/_shapefiles.py b/lib/iris/_shapefiles.py index 06a444d1b1..50db0c3336 100644 --- a/lib/iris/_shapefiles.py +++ b/lib/iris/_shapefiles.py @@ -3,19 +3,15 @@ # 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. -# Much of this code is originally based off the ASCEND library, developed in -# the Met Office by Chris Kent, Emilie Vanvyve, David Bentley, Joana Mendes -# many thanks to them. Converted to iris by Alex Chamberlain-Clay - from __future__ import annotations -import pdb -from typing import overload import importlib from itertools import product import sys +from typing import overload import warnings +from affine import Affine import numpy as np from pyproj import CRS, Transformer import rasterio.features as rfeatures @@ -83,6 +79,8 @@ def create_shapefile_mask( ------ TypeError If the cube is not a :class:`~iris.cube.Cube`. + ValueError + If the :class:`~iris.cube.Cube` has a semi-structured model grid. Notes ----- @@ -90,6 +88,12 @@ def create_shapefile_mask( coordinate reference system (CRS) defined. Masking results will be most reliable when the :class:`iris.cube.Cube` and masking geometry have the same CRS. + Where the :class:`iris.cube.Cube` CRS and the geometry CRS differ, the geometry will be + transformed to the cube CRS using the pyproj library. This is a best-effort + transformation and may not be perfect, especially for complex geometries and + non-standard coordinate referece systems. Consult the `pyproj documentation`_ for + more information. + If the :class:`iris.cube.Cube` has no :class:`~iris.coord_systems`, the default GeogCS is used where the coordinate units are degrees. For any other coordinate units, the cube **must** have a :class:`~iris.coord_systems` defined. @@ -108,6 +112,8 @@ def create_shapefile_mask( See Also -------- :func:`is_geometry_valid` + + .. _`pyproj documentation`: https://pyproj4.github.io/pyproj/stable/api/transformer.html#pyproj-transformer """ # Check validity of geometry CRS is_geometry_valid(geometry, geometry_crs) @@ -122,6 +128,15 @@ def create_shapefile_mask( msg = "Received non-Cube object where a Cube is expected" raise TypeError(msg) + # Get cube coordinates + y_name, x_name = _cube_primary_xy_coord_names(cube) + # Check if cube lons units are in degrees, and if so do they exist in [0, 360] or [-180, 180] + if (cube.coord(x_name).units.origin == "degrees") and ( + cube.coord(x_name).points.max() > 180 + ): + # Convert to [-180, 180] domain + cube = cube.intersection(iris.coords.CoordExtent(x_name, -180, 180)) + # Check for CRS equality and transform if necessary cube_crs = cube.coord_system().as_cartopy_projection() if not geometry_crs.equals(cube_crs): @@ -138,50 +153,18 @@ def create_shapefile_mask( msg = f"Shape geometry is invalid (not well formed): {shapely.is_valid_reason(geometry)}." raise TypeError(msg) - # Get cube coordinates - isLonFix = False - y_name, x_name = _cube_primary_xy_coord_names(cube) - # Check if cube lons units are in degrees, and if so do they exist in [0, 360] or [-180, 180] - if (cube.coord(x_name).units.origin == "degrees") and ( - cube.coord(x_name).points.max() > 180 - ): - isLonFix = True - cube = cube.intersection(iris.coords.CoordExtent(x_name, -180, 180)) - x_points = cube.coord(x_name).points y_points = cube.coord(y_name).points + dx = iris.util.regular_step(cube.coord(x_name)) + dy = iris.util.regular_step(cube.coord(y_name)) w = len(x_points) h = len(y_points) # Define raster transform based on cube # This maps the geometry domain onto the cube domain - # using an Affine transformation - # tr, w, h = rwarp.calculate_default_transform( - # src_crs=cube_crs, - # dst_crs=cube_crs, - # width=w, - # height=h, - # dst_width=w, - # dst_height=h, - # src_geoloc_array=( - # np.meshgrid( - # x_points, - # y_points, - # indexing="ij", - # ) - # ), - # ) - - # tr = rtransform.from_bounds( - # west=x_points.min(), - # south=y_points.min(), - # east=x_points.max(), - # north=y_points.max(), - # width=h, - # height=-w, - # ) - - tr = _transform_from_latlon(x_points, y_points) + trans = Affine.translation(x_points[0] - dx / 2, y_points[0] - dy / 2) + scale = Affine.scale(dx, dy) + tr = trans * scale # Generate mask from geometry mask_template = rfeatures.geometry_mask( @@ -191,13 +174,11 @@ def create_shapefile_mask( **kwargs, ) - # If cube was on [0, 360] domain, then shift mask template - # to match the cube domain - # if isLonFix: - # mask_template = np.roll(mask_template, w // 2, axis=1) + # If cube was on circular domain, then the transformed + # mask template needs shifting to match the cube domain + if cube.coord(x_name).circular: + mask_template = np.roll(mask_template, w // 2, axis=1) - # pdb.set_trace() - # return mask_template[::-1,:] return mask_template @@ -335,19 +316,3 @@ def _cube_primary_xy_coord_names(cube: iris.cube.Cube) -> tuple[str, str]: latitude = latc.name() longitude = lonc.name() return latitude, longitude - - -def _transform_from_latlon(lon, lat): - """perform an affine transformation to the latitude/longitude coordinates""" - - from affine import Affine - - lat = np.asarray(lat) - lon = np.asarray(lon) - - d_lon = lon[1] - lon[0] - d_lat = lat[1] - lat[0] - - trans = Affine.translation(lon[0] - d_lon / 2, lat[0] - d_lat / 2) - scale = Affine.scale(d_lon, d_lat) - return trans * scale diff --git a/lib/iris/util.py b/lib/iris/util.py index ac27fc00f9..115505d0ba 100644 --- a/lib/iris/util.py +++ b/lib/iris/util.py @@ -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 @@ -2128,21 +2128,34 @@ def _strip_metadata_from_dims(cube, dims): return reduced_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 = None, - minimum_weight: float = 0.0, + shape_crs: cartopy.crs | CRS, in_place: bool = False, + **kwargs, ): - """Take a shape object and masks all points not touching it in a cube. + """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. - 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 + 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 ---------- @@ -2154,12 +2167,14 @@ def mask_cube_from_shapefile( 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. 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 ------- @@ -2184,10 +2199,41 @@ def mask_cube_from_shapefile( 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 ------- @@ -2202,7 +2248,7 @@ def mask_cube_from_shapefile( 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, shape_crs) + 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 From 1793e1f25bc514e1e1e6f16154b86c096d00f358 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 14 Feb 2025 16:24:03 +0000 Subject: [PATCH 15/15] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- lib/iris/util.py | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/lib/iris/util.py b/lib/iris/util.py index 115505d0ba..e1b32637b9 100644 --- a/lib/iris/util.py +++ b/lib/iris/util.py @@ -2128,6 +2128,7 @@ def _strip_metadata_from_dims(cube, dims): return reduced_cube + @overload def mask_cube_from_shapefile( cube: iris.cube.Cube, @@ -2136,8 +2137,7 @@ def mask_cube_from_shapefile( in_place: bool = False, all_touched: bool = False, invert: bool = False, -): - ... +): ... def mask_cube_from_shapefile( cube: iris.cube.Cube, shape: shapely.Geometry, @@ -2151,9 +2151,9 @@ def mask_cube_from_shapefile( 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. + Shapes can be polygons, lines or points. - By default, only the only cells whose center is within the polygon or that are selected by + 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. @@ -2208,28 +2208,28 @@ def mask_cube_from_shapefile( >>> 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 + 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: + 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', + 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 >>> wgs84 = CRS.from_epsg(4326) @@ -2248,7 +2248,9 @@ def mask_cube_from_shapefile( If a CRS is not provided for the the masking geometry, the CRS of the cube is assumed. """ - shapefile_mask = create_shapefile_mask(geometry=shape, geometry_crs=shape_crs, cube=cube, **kwargs) + 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