From ebee97757122d257aab1ab35a0f412d1de0016e2 Mon Sep 17 00:00:00 2001 From: Hamish Steptoe Date: Thu, 22 Aug 2024 15:05:12 +0100 Subject: [PATCH 1/3] 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 2/3] 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 3/3] [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: