diff --git a/lib/iris/_shapefiles.py b/lib/iris/_shapefiles.py index 74b24b6627..5bc2bb9d12 100644 --- a/lib/iris/_shapefiles.py +++ b/lib/iris/_shapefiles.py @@ -3,46 +3,74 @@ # 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 +from functools import wraps from itertools import product +import sys import warnings +from affine import Affine +import cartopy.crs as ccrs import numpy as np +from pyproj import CRS, Transformer +import rasterio.features as rfeatures import shapely -import shapely.errors import shapely.geometry as sgeom import shapely.ops -from iris.warnings import IrisDefaultingWarning, IrisUserWarning +import iris +from iris.exceptions import IrisError +from iris.warnings import IrisUserWarning + +# from iris.util import mutually_exclusive_keywords + +if "iris.cube" in sys.modules: + import iris.cube +if "iris.analysis.cartography" in sys.modules: + import iris.analysis.cartography +# @mutually_exclusive_keywords("minimum_weight", "all_touched") def create_shapefile_mask( - geometry, - cube, - minimum_weight=0.0, -): + geometry: shapely.Geometry, + geometry_crs: ccrs.CRS | CRS, + cube: iris.cube.Cube, + minimum_weight: float = 0.0, + **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 ---------- geometry : :class:`shapely.Geometry` + 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 :class:`iris.cube.Cube`. 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). + minimum_weight : float, optional + The minimum weight of the geometry to be included in the mask. + If the weight is less than this value, the geometry will not be + included in the mask. Defaults to 0.0. + + Other Parameters + ---------------- + all_touched : bool, optional + If True, all pixels touched by the geometry will be included in the mask. + If False, only pixels fully covered by the geometry will be included in the mask. + Defaults to True. + invert : bool, optional + If True, the mask will be inverted, so that pixels not covered by the geometry + will be included in the mask. Defaults to False. Returns ------- @@ -50,194 +78,377 @@ 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 + ------ + TypeError + If the cube is not a :class:`~iris.cube.Cube`. + IrisError + If the cube does not have a coordinate reference system defined. + TypeError + If the geometry is not a valid shapely geometry. + ValueError + If the :class:`~iris.cube.Cube` has a semi-structured model grid. + ValueError + If the minimum_weight is not between 0.0 and 1.0. + ValueError + If the minimum_weight is greater than 0.0 and all_touched is True. + + Warns + ----- + IrisUserWarning + If the geometry CRS does not match the cube CRS, and the geometry is transformed + to the cube CRS using pyproj. + + Notes + ----- + 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 :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 reference systems. Consult the `pyproj documentation`_ for + more information. + + If a CRS is not provided for the the masking geometry, the CRS of the :class:`iris.cube.Cube` is assumed. + + Note that `minimum_weight` and `all_touched` are mutually exclusive options: an error will be raised if + a `minimum_weight` > 0 *and* `all_touched` is set to `True`. This is because + `all_touched=True` is equivalent to `minimum_weight=0`. + + Warning + ------- + 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` + + .. _`pyproj documentation`: https://pyproj4.github.io/pyproj/stable/api/transformer.html#pyproj-transformer """ - from iris.cube import Cube, CubeList + # Check validity of geometry CRS + is_geometry_valid(geometry, geometry_crs) - try: - 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): + # 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 - \ 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 cube coordinate system + if not cube.coord_system(): + err_msg = "Cube does not have a coordinate references system defined. For reliable results we recommend you add a coordinate system to your cube." + raise IrisError(err_msg) + + # Check minimum_weight is within range + if (minimum_weight < 0.0) or (minimum_weight > 1.0): + msg = "Minimum weight must be between 0.0 and 1.0" + raise ValueError(msg) + + # Check compatibility of function arguments + # all_touched and minimum_weight are mutually exclusive + if (minimum_weight > 0) and (kwargs.get("all_touched") is True): + msg = "Cannot use minimum_weight > 0.0 with all_touched=True." + raise ValueError(msg) + + # Get cube coordinates + x_coord = cube.coord(axis="x", dim_coords=True) + y_coord = cube.coord(axis="y", dim_coords=True) + # Check if cube lons units are in degrees, and if so do they exist in [0, 360] or [-180, 180] + if (x_coord.units.origin == "degrees") and (x_coord.points.max() > 180): + # Convert to [-180, 180] domain + cube = cube.intersection(iris.coords.CoordExtent(x_coord.name(), -180, 180)) + + if geometry_crs is None: + # If no geometry CRS is provided, assume it is the same as the cube CRS + geometry_crs = cube.coord_system().as_cartopy_projection() + + # Check for CRS equality and transform if necessary + 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=IrisUserWarning) + geometry = _transform_geometry( + geometry=geometry, + geometry_crs=geometry_crs, + cube_crs=cube_crs, ) - # prepare 2D cube - 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 = len(x_coord.points) + h = len(y_coord.points) + # Mask by weight if minimum_weight > 0.0 + if minimum_weight > 0: + mask_template = _get_weighted_mask( + geometry=geometry, + cube=cube, + minimum_weight=minimum_weight, + ) + else: + if (minimum_weight == 0) and (kwargs.get("all_touched") is None): + # For speed, if minimum_weight is 0, then + # we can use the geometry_mask function directly + # This is equivalent to all_touched=True + all_touched = True + elif (minimum_weight == 0) and (kwargs.get("all_touched") is not None): + if not isinstance(kwargs.get("all_touched"), bool): + msg = "`all_touched` kwarg must be True or False." + raise TypeError(msg) + all_touched = kwargs.get("all_touched") + # Define raster transform based on cube + # This maps the geometry domain onto the cube domain + tr = _make_raster_cube_transform(cube) + # Generate mask from geometry + mask_template = rfeatures.geometry_mask( + geometries=shapely.get_parts(geometry), + out_shape=(h, w), + transform=tr, + all_touched=all_touched, + ) + + # If cube was on circular domain, then the transformed + # mask template needs shifting to match the cube domain + if x_coord.circular: + mask_template = np.roll(mask_template, w // 2, axis=1) + + if kwargs.get("invert"): + if not isinstance(kwargs.get("invert"), bool): + msg = "`invert` kwarg must be True or False." + raise TypeError(msg) + # Invert the mask + mask_template = np.logical_not(mask_template) + return mask_template -def _transform_coord_system(geometry, cube, geometry_system=None): - """Project the shape onto another coordinate system. +def is_geometry_valid( + geometry: shapely.Geometry, + geometry_crs: ccrs.CRS | CRS, +) -> 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 antimeridian, + based on the assumption that the shape will + cross the antimeridian if the difference between + the shape bound longitudes is greater than 180. + 4) The geometry does not cross the poles. 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_system : :class:`iris.coord_systems`, optional - A :class:`~iris.coord_systems` object describing - the coord_system of the shapefile. Defaults to None, - which is treated as GeogCS. + 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`. - + 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 antimeridian. + 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 if 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. """ - y_name, x_name = _cube_primary_xy_coord_names(cube) - import iris.analysis.cartography + WGS84_crs = CRS.from_epsg(4326) + + # Check crs is valid type + if not isinstance(geometry_crs, (ccrs.CRS | CRS)): + msg = f"Geometry CRS must be a cartopy.crs or pyproj.CRS object, not {type(geometry_crs)}." + raise TypeError(msg) - DEFAULT_CS = iris.coord_systems.GeogCS( - iris.analysis.cartography.DEFAULT_SPHERICAL_EARTH_RADIUS + # Check geometry is valid shapely geometry + if not shapely.is_valid_input(geometry): + msg = f"Shape geometry is not a valid shape (not well formed)." + raise TypeError(msg) + + # Check that the geometry is within the bounds of the coordinate system + # 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 ) - target_system = cube.coord_system() - if not target_system: - warnings.warn( - "Cube has no coord_system; using default GeogCS lat/lon", - 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) - 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) - - 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 + if not geometry_crs.equals(WGS84_crs): + # 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." + raise ValueError(msg) + + # Check if shape crosses the 180th meridian (or equivalent) + # Exception for MultiPoint geometries where sequential points + # may be separated by more than 180 degrees + if not isinstance(geometry, sgeom.MultiPoint): + if bool(abs(geometry.bounds[2] - geometry.bounds[0]) > 180.0): + msg = "Geometry crossing the antimeridian is not supported." + raise ValueError(msg) + + # Check if the geometry crosses the poles + npole = sgeom.Point(0, 90) + spole = sgeom.Point(0, -90) + if geometry.intersects(npole) or geometry.intersects(spole): + msg = "Geometry crossing the poles is not supported." + raise ValueError(msg) + + return + + +def _transform_geometry( + geometry: shapely.Geometry, geometry_crs: ccrs.CRS | CRS, cube_crs: ccrs.CRS +) -> shapely.Geometry: + """Transform a geometry to the cube CRS using pyproj. + Parameters + ---------- + geometry : :class:`shapely.Geometry` + The geometry to transform. + geometry_crs : :class:`cartopy.crs`, :class:`pyproj.CRS` + The coordinate reference system of the geometry. + cube_crs : :class:`cartopy.crs`, :class:`pyproj.CRS` + The coordinate reference system of the cube. -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 + Returns + ------- + :class:`shapely.Geometry` + The transformed geometry. + """ + # Set-up transform via pyproj + t = Transformer.from_crs( + crs_from=geometry_crs, crs_to=cube_crs, always_xy=True + ).transform + # Transform geometry + return shapely.ops.transform(t, geometry) -def _cube_primary_xy_coord_names(cube): - """Return the primary latitude and longitude coordinate names, or long names, from a cube. +def _get_weighted_mask( + cube: iris.cube.Cube, + geometry: shapely.Geometry, + minimum_weight: float, +) -> np.array: + """Get a mask based on the geometry and minimum weight. Parameters ---------- cube : :class:`iris.cube.Cube` + The cube to mask. + geometry : :class:`shapely.Geometry` + The geometry to use for masking. + minimum_weight : float + The minimum weight of the geometry to be included in the mask. + If the weight is less than this value, the geometry will not be + included in the mask. Returns ------- - tuple of str - The names of the primary latitude and longitude coordinates. - + :class:`np.array` + An array of the shape of the x & y coordinates of the cube, with points + to mask equal to True. """ - latc = ( - cube.coords(axis="y", dim_coords=True)[0] - if cube.coords(axis="y", dim_coords=True) - else -1 - ) - lonc = ( - cube.coords(axis="x", dim_coords=True)[0] - if cube.coords(axis="x", dim_coords=True) - else -1 - ) + if not cube.coord(axis="x", dim_coords=True).has_bounds(): + cube.coord(axis="x", dim_coords=True).guess_bounds() + if not cube.coord(axis="y", dim_coords=True).has_bounds(): + cube.coord(axis="y", dim_coords=True).guess_bounds() + # Get the shape of the cube + w = len(cube.coord(axis="x", dim_coords=True).points) + h = len(cube.coord(axis="y", dim_coords=True).points) + # Get the bounds of the cube + # x_bounds = _get_mod_rebased_coord_bounds(cube.coord(x_name)) + # y_bounds = _get_mod_rebased_coord_bounds(cube.coord(y_name)) + x_bounds = cube.coord(axis="x", dim_coords=True).bounds + y_bounds = cube.coord(axis="y", dim_coords=True).bounds + # Generate Sort-Tile-Recursive (STR) packed R-tree of bounding boxes + # https://shapely.readthedocs.io/en/stable/strtree.html + grid_boxes = [ + sgeom.box(x[0], y[0], x[1], y[1]) for y, x in product(y_bounds, x_bounds) + ] + grid_tree = shapely.STRtree(grid_boxes) + # Find grid boxes indexes that intersect with the geometry + idxs = grid_tree.query(geometry, predicate="intersects") + # Get grid box indexes that intersect with the minimum weight criteria + mask_idxs_bool = [ + grid_boxes[idx].intersection(geometry).area / grid_boxes[idx].area + >= minimum_weight + for idx in idxs + ] + mask_idxs = idxs[mask_idxs_bool] + mask_xy = [list(product(range(h), range(w)))[i] for i in mask_idxs] + # Create mask from grid box indexes - if -1 in (latc, lonc): - msg = "Error retrieving 1d xy coordinates in cube: {!r}" - raise ValueError(msg.format(cube)) + weighted_mask_template = np.ones((h, w), dtype=bool) + # Set mask = True for grid box indexes identified above + for xy in mask_xy: + weighted_mask_template[xy] = False + return weighted_mask_template - latitude = latc.name() - longitude = lonc.name() - return latitude, longitude +def _make_raster_cube_transform(cube: iris.cube.Cube) -> Affine: + """Create a rasterio transform for the cube. -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. + Raises + ------ + CoordinateNotRegularError + If the cube dimension coordinates are not regular, + such that :func:`iris.util.regular_step` returns an error. Returns ------- - :class:`np.array` - A 1d Numpy array of [start,end] pairs for bounds of the coord. - + :class:`affine.Affine` + An affine transform object that maps the geometry domain onto the cube domain. """ - 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 + x_points = cube.coord(axis="x", dim_coords=True).points + y_points = cube.coord(axis="y", dim_coords=True).points + dx = iris.util.regular_step(cube.coord(axis="x", dim_coords=True)) + dy = iris.util.regular_step(cube.coord(axis="y", dim_coords=True)) + # Create a rasterio transform based on the cube + # This maps the geometry domain onto the cube domain + trans = Affine.translation(x_points[0] - dx / 2, y_points[0] - dy / 2) + scale = Affine.scale(dx, dy) + return trans * scale diff --git a/lib/iris/tests/integration/test_mask_cube_from_shapefile.py b/lib/iris/tests/integration/test_mask_cube_from_shapefile.py index 52fd02615d..30d3e8f848 100644 --- a/lib/iris/tests/integration/test_mask_cube_from_shapefile.py +++ b/lib/iris/tests/integration/test_mask_cube_from_shapefile.py @@ -6,104 +6,195 @@ import math +import cartopy.crs as ccrs import cartopy.io.shapereader as shpreader import numpy as np +from pyproj import CRS +import pytest +from pytest import approx +from shapely.geometry import LineString, MultiLineString, MultiPoint, Point import iris +from iris.coord_systems import GeogCS import iris.tests as tests from iris.util import mask_cube_from_shapefile +wgs84 = CRS.from_epsg(4326) -@tests.skip_data -class TestCubeMasking(tests.IrisTest): - """integration tests of mask_cube_from_shapefile - using different projections in iris_test_data - - values are the KGO calculated using ASCEND. - """ - - def setUp(self): - ne_countries = shpreader.natural_earth( - resolution="10m", category="cultural", name="admin_0_countries" - ) - self.reader = shpreader.Reader(ne_countries) - - def test_global_proj_russia(self): - path = tests.get_data_path( - ["NetCDF", "global", "xyt", "SMALL_hires_wind_u_for_ipcc4.nc"] - ) - test_global = iris.load_cube(path) - ne_russia = [ - country.geometry - for country in self.reader.records() - if "Russia" in country.attributes["NAME_LONG"] - ][0] - masked_test = mask_cube_from_shapefile(test_global, ne_russia) - print(np.sum(masked_test.data)) - assert math.isclose(np.sum(masked_test.data), 76845.37, rel_tol=0.001), ( - "Global data with Russia mask failed test" - ) - - def test_rotated_pole_proj_germany(self): - path = tests.get_data_path( - ["NetCDF", "rotated", "xy", "rotPole_landAreaFraction.nc"] - ) - test_rotated = iris.load_cube(path) - ne_germany = [ - country.geometry - for country in self.reader.records() - if "Germany" in country.attributes["NAME_LONG"] - ][0] - masked_test = mask_cube_from_shapefile(test_rotated, ne_germany) - assert math.isclose(np.sum(masked_test.data), 179.46872, rel_tol=0.001), ( - "rotated europe data with German mask failed test" - ) - - def test_transverse_mercator_proj_uk(self): - path = tests.get_data_path( - ["NetCDF", "transverse_mercator", "tmean_1910_1910.nc"] - ) - test_transverse = iris.load_cube(path) - ne_uk = [ - country.geometry - for country in self.reader.records() - if "United Kingdom" in country.attributes["NAME_LONG"] - ][0] - masked_test = mask_cube_from_shapefile(test_transverse, ne_uk) - assert math.isclose(np.sum(masked_test.data), 90740.25, rel_tol=0.001), ( - "transverse mercator UK data with UK mask failed test" - ) - - def test_rotated_pole_proj_germany_weighted_area(self): - path = tests.get_data_path( - ["NetCDF", "rotated", "xy", "rotPole_landAreaFraction.nc"] - ) - test_rotated = iris.load_cube(path) - ne_germany = [ - country.geometry - for country in self.reader.records() - if "Germany" in country.attributes["NAME_LONG"] - ][0] - masked_test = mask_cube_from_shapefile( - test_rotated, ne_germany, minimum_weight=0.9 - ) - assert math.isclose(np.sum(masked_test.data), 125.60199, rel_tol=0.001), ( - "rotated europe data with 0.9 weight germany mask failed test" - ) - - def test_4d_global_proj_brazil(self): - path = tests.get_data_path(["NetCDF", "global", "xyz_t", "GEMS_CO2_Apr2006.nc"]) - test_4d_brazil = iris.load_cube(path, "Carbon Dioxide") - ne_brazil = [ - country.geometry - for country in self.reader.records() - if "Brazil" in country.attributes["NAME_LONG"] - ][0] - masked_test = mask_cube_from_shapefile( - test_4d_brazil, - ne_brazil, - ) - print(np.sum(masked_test.data)) - # breakpoint() - assert math.isclose(np.sum(masked_test.data), 18616921.2, rel_tol=0.001), ( - "4d data with brazil mask failed test" - ) +ne_countries = shpreader.natural_earth( + resolution="10m", category="cultural", name="admin_0_countries" +) +reader = shpreader.Reader(ne_countries) + + +@pytest.mark.parametrize( + ("minimum_weight", "all_touched", "invert", "expected_sum"), + [ + (0.0, None, None, 10522684.77), # Minimum weight == 0 + (0.0, None, False, 10522684.77), # Minimum weight == 0 + (0.0, True, False, 10522684.77), # All touched == True + (0.5, None, False, 8965584.47), # Minimum weight == 0.5 + (1.0, None, False, 7504361.29), # Minimum weight == 1 + (0.0, False, False, 8953582.05), # All touched == False + (0.0, True, True, 605637718.12), # All touched == True, Invert == True + ], +) +def test_global_proj_china(minimum_weight, all_touched, invert, expected_sum): + """Test masking with a shape for China with various parameter combinations.""" + path = tests.get_data_path(["NetCDF", "global", "xyt", "SMALL_total_column_co2.nc"]) + test_global = iris.load_cube(path) # Crop to avoid edge effects + test_global.coord("latitude").coord_system = GeogCS(6371229) + test_global.coord("longitude").coord_system = GeogCS(6371229) + ne_china = [ + country.geometry + for country in reader.records() + if "China" in country.attributes["NAME_LONG"] + ][0] + masked_test = mask_cube_from_shapefile( + test_global, + ne_china, + shape_crs=wgs84, + minimum_weight=minimum_weight, + all_touched=all_touched, + invert=invert, + ) + assert masked_test.ndim == 3 + assert approx(np.sum(masked_test.data), rel=0.001) == expected_sum + + +def test_global_proj_russia(): + """Test masking with a shape that crosses the antimeridian.""" + path = tests.get_data_path(["NetCDF", "global", "xyt", "SMALL_total_column_co2.nc"]) + test_global = iris.load_cube(path) + test_global.coord("latitude").coord_system = GeogCS(6371229) + test_global.coord("longitude").coord_system = GeogCS(6371229) + ne_russia = [ + country.geometry + for country in reader.records() + if "Russia" in country.attributes["NAME_LONG"] + ][0] + + with pytest.raises( + ValueError, match="Geometry crossing the antimeridian is not supported." + ): + mask_cube_from_shapefile(test_global, ne_russia, shape_crs=wgs84) + + +def test_rotated_pole_proj_uk(): + """Test masking a rotated pole projection cube for the UK with lat/lon shape.""" + path = tests.get_data_path( + ["NetCDF", "rotated", "xy", "rotPole_landAreaFraction.nc"] + ) + test_rotated = iris.load_cube(path) + ne_uk = [ + country.geometry + for country in reader.records() + if "United Kingdom" in country.attributes["NAME_LONG"] + ][0] + masked_test = mask_cube_from_shapefile(test_rotated, ne_uk, shape_crs=wgs84) + assert masked_test.ndim == 2 + assert approx(np.sum(masked_test.data), rel=0.001) == 102.77 + + +def test_transverse_mercator_proj_uk(): + """Test masking a transverse mercator projection cube for the UK with lat/lon shape.""" + path = tests.get_data_path(["NetCDF", "transverse_mercator", "tmean_1910_1910.nc"]) + test_transverse = iris.load_cube(path) + ne_uk = [ + country.geometry + for country in reader.records() + if "United Kingdom" in country.attributes["NAME_LONG"] + ][0] + masked_test = mask_cube_from_shapefile(test_transverse, ne_uk, shape_crs=wgs84) + assert masked_test.ndim == 3 + assert approx(np.sum(masked_test.data), rel=0.001) == 90740.25 + + +def test_rotated_pole_proj_germany_weighted_area(): + """Test masking a rotated pole projection cube for Germany with weighted area.""" + path = tests.get_data_path( + ["NetCDF", "rotated", "xy", "rotPole_landAreaFraction.nc"] + ) + test_rotated = iris.load_cube(path) + ne_germany = [ + country.geometry + for country in reader.records() + if "Germany" in country.attributes["NAME_LONG"] + ][0] + masked_test = mask_cube_from_shapefile( + test_rotated, ne_germany, shape_crs=wgs84, minimum_weight=0.9 + ) + assert masked_test.ndim == 2 + assert approx(np.sum(masked_test.data), rel=0.001) == 125.60199 + + +def test_4d_global_proj_brazil(): + """Test masking a 4D global projection cube for Brazil with lat/lon shape.""" + path = tests.get_data_path(["NetCDF", "global", "xyz_t", "GEMS_CO2_Apr2006.nc"]) + test_4d_brazil = iris.load_cube(path, "Carbon Dioxide") + test_4d_brazil.coord("latitude").coord_system = GeogCS(6371229) + test_4d_brazil.coord("longitude").coord_system = GeogCS(6371229) + ne_brazil = [ + country.geometry + for country in reader.records() + if "Brazil" in country.attributes["NAME_LONG"] + ][0] + masked_test = mask_cube_from_shapefile( + test_4d_brazil, ne_brazil, shape_crs=wgs84, all_touched=True + ) + assert masked_test.ndim == 4 + assert approx(np.sum(masked_test.data), rel=0.001) == 18616921.2 + + +@pytest.mark.parametrize( + ("shape", "expected_value"), + [ + (Point(-3.475446894622651, 50.72770791320487), 12061.74), # (x,y) + ( + LineString( + [ + (-5.712431305030631, 50.06590599588483), + (-3.0704940433528947, 58.644091639685456), + ] + ), + 120530.41, + ), # (x,y) to (x,y) + ( + MultiPoint( + [ + (-5.712431305030631, 50.06590599588483), + (-3.0704940433528947, 58.644091639685456), + ] + ), + 24097.47, + ), + ( + MultiLineString( + [ + [ + (-5.206405826948041, 49.95891620303525), + (-3.376975634580173, 58.67197421392852), + ], + [ + (-6.2276386132877475, 56.71561805509071), + (1.7626540441873777, 52.48118683241357), + ], + ] + ), + 253248.44, + ), + ], +) +def test_global_proj_uk_shapes(shape, expected_value): + """Test masking with a variety of shape types.""" + path = tests.get_data_path(["NetCDF", "global", "xyt", "SMALL_total_column_co2.nc"]) + test_global = iris.load_cube(path) # Crop to avoid edge effects + test_global.coord("latitude").coord_system = GeogCS(6371229) + test_global.coord("longitude").coord_system = GeogCS(6371229) + masked_test = mask_cube_from_shapefile( + test_global, + shape, + shape_crs=wgs84, + ) + assert masked_test.ndim == 3 + assert approx(np.sum(masked_test.data), rel=0.001) == expected_value 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_create_shapefile_mask.py b/lib/iris/tests/unit/_shapefiles/test_create_shapefile_mask.py new file mode 100644 index 0000000000..96706d7d12 --- /dev/null +++ b/lib/iris/tests/unit/_shapefiles/test_create_shapefile_mask.py @@ -0,0 +1,225 @@ +# 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.create_shapefile_mask`.""" + +# import iris tests first so that some things can be initialised before +# importing anything else +import iris.tests as tests # isort:skip + +import numpy as np +from pyproj import CRS +import pytest +from shapely.geometry import Point, Polygon + +from iris._shapefiles import create_shapefile_mask +from iris.coord_systems import GeogCS +from iris.coords import DimCoord +from iris.cube import Cube, CubeList +from iris.exceptions import IrisError +from iris.warnings import IrisUserWarning + + +@pytest.fixture(scope="session") +def wgs84_crs(): + return CRS.from_epsg(4326) + + +@pytest.fixture +def square_polygon(): + # Create a a 3x3 square polygon from (3,3) to (6,6) using shapely + return Polygon([(3, 3), (6, 3), (6, 6), (3, 6)]) + + +@pytest.fixture +def circle_polygon(): + # Create a a circular polygon centred on (5,5) with radius (2,) using shapely + return Point(5, 5).buffer(2) + + +@pytest.fixture +def mock_cube(): + """Create a mock 10x10 Iris cube for testing.""" + x_points = np.linspace(0, 9, 10) + y_points = np.linspace(0, 9, 10) + x_coord = DimCoord( + x_points, + standard_name="longitude", + units="degrees", + coord_system=GeogCS(6371229), + ) + y_coord = DimCoord( + y_points, + standard_name="latitude", + units="degrees", + coord_system=GeogCS(6371229), + ) + data = np.ones((len(y_points), len(x_points))) + cube = Cube(data, dim_coords_and_dims=[(y_coord, 0), (x_coord, 1)]) + return cube + + +def test_basic_create_shapefile_mask(square_polygon, wgs84_crs, mock_cube): + """Test the create_shapefile_mask function.""" + # Create a mask using the square polygon + mask = create_shapefile_mask(square_polygon, wgs84_crs, mock_cube) + + # Check that the mask is a boolean array with the same shape as the cube data + assert mask.shape == mock_cube.data.shape + assert mask.dtype == np.bool_ + + # Check that the masked area corresponds to the square polygon + expected_mask = np.ones_like(mock_cube.data, dtype=bool) + expected_mask[3:6, 3:6] = False # The square polygon covers this area + assert np.array_equal(mask, expected_mask) + + +def test_invert_create_shapefile_mask(square_polygon, wgs84_crs, mock_cube): + """Test the create_shapefile_mask function.""" + # Create a mask using the square polygon + mask = create_shapefile_mask(square_polygon, wgs84_crs, mock_cube, invert=True) + + # Check that the mask is a boolean array with the same shape as the cube data + assert mask.shape == mock_cube.data.shape + assert mask.dtype == np.bool_ + + # Check that the masked area corresponds to the square polygon + expected_mask = np.zeros_like(mock_cube.data, dtype=bool) + expected_mask[3:6, 3:6] = True # The square polygon covers this area + assert np.array_equal(mask, expected_mask) + + +def test_all_touched_true_create_shapefile_mask(circle_polygon, wgs84_crs, mock_cube): + """Test the create_shapefile_mask function.""" + # Create a mask using the square polygon + mask = create_shapefile_mask(circle_polygon, wgs84_crs, mock_cube, all_touched=True) + + # Check that the mask is a boolean array with the same shape as the cube data + assert mask.shape == mock_cube.data.shape + assert mask.dtype == np.bool_ + + # Check that the masked area corresponds to the circular polygon + expected_mask = np.array( + [ + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 0, 0, 0, 1, 1, 1], + [1, 1, 1, 0, 0, 0, 0, 0, 1, 1], + [1, 1, 1, 0, 0, 0, 0, 0, 1, 1], + [1, 1, 1, 0, 0, 0, 0, 0, 1, 1], + [1, 1, 1, 1, 0, 0, 0, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + ], + dtype=bool, + ) + assert np.array_equal(mask, expected_mask) + + +def test_all_touched_false_create_shapefile_mask(circle_polygon, wgs84_crs, mock_cube): + """Test the create_shapefile_mask function.""" + # Create a mask using the square polygon + mask = create_shapefile_mask( + circle_polygon, wgs84_crs, mock_cube, all_touched=False + ) + + # Check that the mask is a boolean array with the same shape as the cube data + assert mask.shape == mock_cube.data.shape + assert mask.dtype == np.bool_ + + # Check that the masked area corresponds to the circular polygon + expected_mask = np.array( + [ + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 0, 1, 1, 1, 1], + [1, 1, 1, 1, 0, 0, 0, 1, 1, 1], + [1, 1, 1, 0, 0, 0, 0, 1, 1, 1], + [1, 1, 1, 1, 0, 0, 0, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + ], + dtype=bool, + ) + assert np.array_equal(mask, expected_mask) + + +def test_create_shapefile_mask_(square_polygon, wgs84_crs, mock_cube): + """Test the create_shapefile_mask function.""" + # Create a mask using the square polygon + mask = create_shapefile_mask(square_polygon, wgs84_crs, mock_cube, invert=True) + + # Check that the mask is a boolean array with the same shape as the cube data + assert mask.shape == mock_cube.data.shape + assert mask.dtype == np.bool_ + + # Check that the masked area corresponds to the square polygon + expected_mask = np.zeros_like(mock_cube.data, dtype=bool) + expected_mask[3:6, 3:6] = True # The square polygon covers this area + assert np.array_equal(mask, expected_mask) + + +class TestCreateShapefileMaskErrors: + def test_invalid_polygon_type(self, wgs84_crs, mock_cube): + # Pass an invalid geometry type (e.g., a string) + with pytest.raises(TypeError): + create_shapefile_mask("not_a_polygon", wgs84_crs, mock_cube) + + def test_invalid_crs_type(self, square_polygon, mock_cube): + # Pass an invalid CRS type (e.g., a string) + with pytest.raises(TypeError): + create_shapefile_mask(square_polygon, "not_a_crs", mock_cube) + + def test_invalid_cube_type(self, square_polygon, wgs84_crs): + # Pass an invalid cube type (e.g., a string or CubeList) + with pytest.raises(TypeError): + create_shapefile_mask(square_polygon, wgs84_crs, "not_a_cube") + with pytest.raises(TypeError): + create_shapefile_mask(square_polygon, wgs84_crs, CubeList()) + + def test_invalid_cube_crs(self, square_polygon, wgs84_crs): + # Pass a cube without a coordinate system + cube = Cube(np.ones((10, 10)), dim_coords_and_dims=[]) + with pytest.raises(IrisError): + create_shapefile_mask(square_polygon, wgs84_crs, cube) + + def test_invalid_minimum_weight(self, square_polygon, wgs84_crs): + # Pass invalid minimum_weight values + with pytest.raises(TypeError): + create_shapefile_mask( + square_polygon, wgs84_crs, mock_cube, minimum_weight="not_a_number" + ) + with pytest.raises(TypeError): + create_shapefile_mask( + square_polygon, wgs84_crs, mock_cube, minimum_weight=-1 + ) + with pytest.raises(TypeError): + create_shapefile_mask( + square_polygon, wgs84_crs, mock_cube, minimum_weight=2 + ) + + def test_invalid_args(self, square_polygon, wgs84_crs, mock_cube): + # Pass invalid minimum_weight values + with pytest.raises(ValueError): + create_shapefile_mask( + square_polygon, + wgs84_crs, + mock_cube, + minimum_weight=0.5, + all_touched=True, + ) + + def test_incompatible_crs_warning(self, square_polygon, mock_cube): + # Pass a CRS that does not match the cube's CRS + crs = CRS.from_epsg(3857) # Web Mercator, different from WGS84 + warn_message = "Geometry CRS does not match cube CRS. Iris will attempt to transform the geometry onto the cube CRS..." + with pytest.warns(IrisUserWarning, match=warn_message): + create_shapefile_mask(square_polygon, crs, mock_cube) + + +# Note: `minimum_weight` keyword argument is tested under its' own unit test +# `test_mask_cube_from_shapefile.py` in the `lib/iris/tests/unit/_shapefiles/` directory. diff --git a/lib/iris/tests/unit/_shapefiles/test_get_weighted_mask.py b/lib/iris/tests/unit/_shapefiles/test_get_weighted_mask.py new file mode 100644 index 0000000000..ed1b6b1dbb --- /dev/null +++ b/lib/iris/tests/unit/_shapefiles/test_get_weighted_mask.py @@ -0,0 +1,113 @@ +# 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._get_weighted_mask`.""" + +# import iris tests first so that some things can be initialised before +# importing anything else +import iris.tests as tests # isort:skip + +import numpy as np +import pytest +from shapely.geometry import box + +from iris._shapefiles import _get_weighted_mask +from iris.coord_systems import GeogCS +from iris.coords import DimCoord +from iris.cube import Cube + + +@pytest.fixture +def square_polygon(): + # Create a roughly 3x3 square polygon + return box(2.4, 2.4, 6.4, 6.4) + + +@pytest.fixture +def mock_cube(): + """Create a mock 9x9 Iris cube for testing.""" + x_points = np.linspace(1, 9, 9) - 0.5 # Specify cube cell midpoints + y_points = np.linspace(1, 9, 9) - 0.5 + x_coord = DimCoord( + x_points, + standard_name="longitude", + units="degrees", + coord_system=GeogCS(6371229), + ) + y_coord = DimCoord( + y_points, + standard_name="latitude", + units="degrees", + coord_system=GeogCS(6371229), + ) + data = np.ones((len(y_points), len(x_points))) + cube = Cube(data, dim_coords_and_dims=[(y_coord, 0), (x_coord, 1)]) + return cube + + +@pytest.mark.parametrize( + "minimum_weight, expected_mask", + [ + ( + 0.0, + np.array( + [ + [1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 0, 0, 0, 0, 0, 1, 1], + [1, 1, 0, 0, 0, 0, 0, 1, 1], + [1, 1, 0, 0, 0, 0, 0, 1, 1], + [1, 1, 0, 0, 0, 0, 0, 1, 1], + [1, 1, 0, 0, 0, 0, 0, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1], + ], + dtype=bool, + ), + ), + ( + 0.5, + np.array( + [ + [1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 0, 0, 0, 1, 1, 1], + [1, 1, 0, 0, 0, 0, 1, 1, 1], + [1, 1, 0, 0, 0, 0, 1, 1, 1], + [1, 1, 0, 0, 0, 0, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1], + ], + dtype=bool, + ), + ), + ( + 1.0, + np.array( + [ + [1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 0, 0, 0, 1, 1, 1], + [1, 1, 1, 0, 0, 0, 1, 1, 1], + [1, 1, 1, 0, 0, 0, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1], + ], + dtype=bool, + ), + ), + ], +) +def test_basic_create_shapefile_mask( + mock_cube, square_polygon, minimum_weight, expected_mask +): + """Test the create_shapefile_mask function with different minimum weights.""" + # Create a mask using the square polygon + mask = _get_weighted_mask(mock_cube, square_polygon, minimum_weight=minimum_weight) + + # Check that the masked area corresponds to the square polygon + assert np.array_equal(mask, expected_mask) 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..b9ec1c44dc --- /dev/null +++ b/lib/iris/tests/unit/_shapefiles/test_is_geometry_valid.py @@ -0,0 +1,197 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the BSD license. +# See LICENSE in the root of the repository for full licensing details. +"""Unit tests for :func:`iris._shapefiles.is_geometry_valid`.""" + +# import iris tests first so that some things can be initialised before +# importing anything else +import iris.tests as tests # isort:skip + +from pyproj import CRS +import pytest +from shapely.geometry import ( + LineString, + MultiLineString, + MultiPoint, + MultiPolygon, + Point, + box, +) + +from iris._shapefiles import is_geometry_valid + + +# Shareable shape fixtures used in: +# - util/test_mask_cube_from_shapefile.py +# - _shapefiles/test_is_geometry_valid.py +@pytest.fixture(scope="session") +def wgs84_crs(): + return CRS.from_epsg(4326) + + +@pytest.fixture(scope="session") +def osgb_crs(): + return CRS.from_epsg(27700) + + +@pytest.fixture(scope="session") +def basic_polygon_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_multipolygon_geometry(): + # Define the coordinates of a basic rectangle + min_lon = 0 + min_lat = 0 + max_lon = 8 + max_lat = 8 + + # Create the rectangular geometry + return MultiPolygon( + [ + box(min_lon, min_lat, max_lon, max_lat), + box(min_lon + 10, min_lat + 10, max_lon + 10, max_lat + 10), + ] + ) + + +@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_multiline_geometry(): + # Define the coordinates of a basic line + return MultiLineString([[(0, 0), (10, 10)], [(20, 20), (30, 30)]]) + + +@pytest.fixture(scope="session") +def basic_point_collection(): + # Define the coordinates of a basic collection of points + # as (lon, lat) tuples, assuming a WGS84 projection. + points = MultiPoint( + [ + (0, 0), + (10, 10), + (-10, -10), + (-3.476204, 50.727059), + (174.761067, -36.846211), + (-77.032801, 38.892717), + ] + ) + + return points + + +@pytest.fixture(scope="session") +def canada_geometry(): + # Define the coordinates of a rectangle that covers Canada + return box(-143.5, 42.6, -37.8, 84.0) + + +@pytest.fixture(scope="session") +def bering_sea_geometry(): + # Define the coordinates of a rectangle that covers the Bering Sea + return box(148.42, 49.1, -138.74, 73.12) + + +@pytest.fixture(scope="session") +def uk_geometry(): + # Define the coordinates of a rectangle that covers the UK + return box(-10, 49, 2, 61) + + +@pytest.fixture(scope="session") +def invalid_geometry_poles(): + # Define the coordinates of a rectangle that crosses the poles + return box(-10, -90, 10, 90) + + +@pytest.fixture(scope="session") +def invalid_geometry_bounds(): + # Define the coordinates of a rectangle that is outside the bounds of the coordinate system + return box(-200, -100, 200, 100) + + +@pytest.fixture(scope="session") +def not_a_valid_geometry(): + # Return an invalid geometry type + # This is not a valid geometry, e.g., a string + return "This is not a valid geometry" + + +# Test validity of different geometries +@pytest.mark.parametrize( + "test_input", + [ + "basic_polygon_geometry", + "basic_multipolygon_geometry", + "basic_point_geometry", + "basic_point_collection", + "basic_line_geometry", + "basic_multiline_geometry", + "canada_geometry", + ], +) +def test_valid_geometry(test_input, request, wgs84_crs): + # Assert that all valid geometries are return None + assert is_geometry_valid(request.getfixturevalue(test_input), wgs84_crs) is None + + +# Fixtures retrieved from conftest.py +# N.B. error message comparison is done with regex so +# any parentheses in the error message must be escaped (\) +@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, + r"Geometry \[\] is not valid for the given coordinate system EPSG:4326. \nCheck that your coordinates are correctly specified.", + ), + ( + "not_a_valid_geometry", + TypeError, + r"Shape geometry is not a valid shape \(not well formed\).", + ), + ], +) +def test_invalid_geometry(test_input, errortype, errormessage, request, wgs84_crs): + # Assert that all invalid geometries raise the expected error + with pytest.raises(errortype, match=errormessage): + is_geometry_valid(request.getfixturevalue(test_input), wgs84_crs) + + +def test_invalid_crs(basic_polygon_geometry): + # Assert that an invalid crs raise the expected error + msg = "Geometry CRS must be a cartopy.crs or pyproj.CRS object, not ." + # Using a string as an invalid CRS type + with pytest.raises(TypeError, match=msg): + is_geometry_valid(basic_polygon_geometry, "invalid_crs") diff --git a/lib/iris/tests/unit/_shapefiles/test_make_raster_cube_transform.py b/lib/iris/tests/unit/_shapefiles/test_make_raster_cube_transform.py new file mode 100644 index 0000000000..b682708d38 --- /dev/null +++ b/lib/iris/tests/unit/_shapefiles/test_make_raster_cube_transform.py @@ -0,0 +1,72 @@ +# 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._make_raster_cube_transform`.""" + +# import iris tests first so that some things can be initialised before +# importing anything else +import iris.tests as tests # isort:skip + +from affine import Affine +import numpy as np +import pytest + +import iris +from iris._shapefiles import _make_raster_cube_transform +from iris.coords import DimCoord +from iris.cube import Cube +from iris.exceptions import CoordinateNotRegularError +from iris.util import regular_step + + +@pytest.fixture +def mock_cube(): + """Create a mock Iris cube for testing.""" + x_points = np.array([0.0, 1.0, 2.0, 3.0]) + y_points = np.array([0.0, 1.0, 2.0, 3.0]) + x_coord = DimCoord(x_points, standard_name="longitude", units="degrees") + y_coord = DimCoord(y_points, standard_name="latitude", units="degrees") + data = np.zeros((len(y_points), len(x_points))) + cube = Cube(data, dim_coords_and_dims=[(y_coord, 0), (x_coord, 1)]) + return cube + + +@pytest.fixture +def mock_nonregular_cube(): + """Create a mock Iris cube for testing.""" + x_points = np.array( + [0.0, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 9.5, 10.0] + ) + y_points = np.array( + [0.0, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 9.5, 10.0] + ) + x_coord = DimCoord(x_points, standard_name="longitude", units="degrees") + y_coord = DimCoord(y_points, standard_name="latitude", units="degrees") + data = np.zeros((len(y_points), len(x_points))) + cube = Cube(data, dim_coords_and_dims=[(y_coord, 0), (x_coord, 1)]) + return cube + + +def test_make_raster_cube_transform(mock_cube): + """Test the `_make_raster_cube__transform` function.""" + x_name = "longitude" + y_name = "latitude" + + # Call the function + transform = _make_raster_cube_transform(mock_cube) + + # Validate the result + dx = regular_step(mock_cube.coord(x_name)) + dy = regular_step(mock_cube.coord(y_name)) + expected_transform = Affine.translation(-dx / 2, -dy / 2) * Affine.scale(dx, dy) + + assert isinstance(transform, Affine) + assert transform == expected_transform + + +def test_invalid_cube(mock_nonregular_cube): + # Assert that all invalid geometries raise the expected error + errormessage = "Coord longitude is not regular" + with pytest.raises(CoordinateNotRegularError, match=errormessage): + _make_raster_cube_transform(mock_nonregular_cube) diff --git a/lib/iris/tests/unit/_shapefiles/test_transform_geometry.py b/lib/iris/tests/unit/_shapefiles/test_transform_geometry.py new file mode 100644 index 0000000000..e50269da9b --- /dev/null +++ b/lib/iris/tests/unit/_shapefiles/test_transform_geometry.py @@ -0,0 +1,100 @@ +from pyproj import CRS +from pyproj import exceptions as pyproj_exceptions +import pytest +import shapely +from shapely.geometry import LineString, Point, Polygon + +import iris +from iris._shapefiles import _transform_geometry +from iris.tests import _shared_utils, stock + +wgs84 = CRS.from_epsg(4326) # WGS84 coordinate system +osgb = CRS.from_epsg(27700) # OSGB coordinate system + + +@pytest.mark.parametrize( + "input_geometry, input_geometry_crs, input_cube_crs, output_expected_geometry", + [ + ( # Basic geometry in WGS84, no transformation needed + shapely.geometry.box(-10, 50, 2, 60), + wgs84, + stock.simple_pp().coord_system()._crs, + shapely.geometry.box(-10, 50, 2, 60), + ), + ( # Basic geometry in WGS84, transformed to OSGB + shapely.geometry.box(-10, 50, 2, 60), + wgs84, + iris.load_cube( + _shared_utils.get_data_path( + ("NetCDF", "transverse_mercator", "tmean_1910_1910.nc") + ) + ) + .coord_system() + .as_cartopy_projection(), + Polygon( + [ + (686600.5247600826, 18834.835866007765), + (622998.2965261643, 1130592.5248690124), + (-45450.063023168186, 1150844.9676151862), + (-172954.59474739246, 41898.60193228102), + (686600.5247600826, 18834.835866007765), + ] + ), + ), + ( # Basic geometry in WGS84, no transformation needed + LineString([(-10, 50), (2, 60)]), + wgs84, + stock.simple_pp().coord_system()._crs, + LineString([(-10, 50), (2, 60)]), + ), + ( # Basic geometry in WGS84, no transformation needed + Point((-10, 50)), + wgs84, + stock.simple_pp().coord_system()._crs, + Point((-10, 50)), + ), + ], +) +def test_transform_geometry( + input_geometry, + input_geometry_crs, + input_cube_crs, + output_expected_geometry, +): + # Assert that all invalid geometries raise the expected error + out_geometry = _transform_geometry( + input_geometry, input_geometry_crs, input_cube_crs + ) + assert isinstance(out_geometry, shapely.geometry.base.BaseGeometry) + assert output_expected_geometry == out_geometry + + +# Assert that an invalid inputs raise the expected errors +@pytest.mark.parametrize( + "input_geometry, input_geometry_crs, input_cube_crs, expected_error", + [ + ( # Basic geometry in WGS84, no transformation needed + "bad_input_geometry", + wgs84, + stock.simple_pp().coord_system()._crs, + AttributeError, + ), + ( # Basic geometry in WGS84, no transformation needed + shapely.geometry.box(-10, 50, 2, 60), + "bad_input_crs", + stock.simple_pp().coord_system()._crs, + pyproj_exceptions.CRSError, + ), + ( # Basic geometry in WGS84, no transformation needed + shapely.geometry.box(-10, 50, 2, 60), + wgs84, + "bad_input_cube_crs", + pyproj_exceptions.CRSError, + ), + ], +) +def test_transform_geometry_invalid_input( + input_geometry, input_geometry_crs, input_cube_crs, expected_error +): + with pytest.raises(expected_error): + _transform_geometry(input_geometry, input_geometry_crs, input_cube_crs) 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 0bd2afda21..b9520cd67a 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,117 +5,51 @@ """Unit tests for :func:`iris.util.mask_cube_from_shapefile`.""" import numpy as np +from pyproj import CRS import pytest -import shapely +from shapely.geometry import box -from iris.coord_systems import RotatedGeogCS +from iris.coord_systems import GeogCS from iris.coords import DimCoord -import iris.cube +from iris.cube import Cube from iris.util import mask_cube_from_shapefile -from iris.warnings import IrisUserWarning -class TestBasicCubeMasking: - """Unit tests for mask_cube_from_shapefile function.""" - - @pytest.fixture(autouse=True) - def _setup(self): - basic_data = np.array([[1, 2, 3], [4, 8, 12]]) - self.basic_cube = iris.cube.Cube(basic_data) - coord = DimCoord( - np.array([0, 1.0]), - standard_name="projection_y_coordinate", - bounds=[[0, 0.5], [0.5, 1]], - units="1", - ) - self.basic_cube.add_dim_coord(coord, 0) - coord = DimCoord( - np.array([0, 1.0, 1.5]), - standard_name="projection_x_coordinate", - bounds=[[0, 0.5], [0.5, 1], [1, 1.5]], - units="1", - ) - self.basic_cube.add_dim_coord(coord, 1) - - def test_basic_cube_intersect(self): - shape = shapely.geometry.box(0.6, 0.6, 0.9, 0.9) - masked_cube = mask_cube_from_shapefile(self.basic_cube, shape) - assert np.sum(masked_cube.data) == 8, ( - f"basic cube masking failed test - expected 8 got {np.sum(masked_cube.data)}" - ) - - def test_basic_cube_intersect_in_place(self): - shape = shapely.geometry.box(0.6, 0.6, 0.9, 0.9) - cube = self.basic_cube.copy() - mask_cube_from_shapefile(cube, shape, in_place=True) - assert np.sum(cube.data) == 8, ( - f"basic cube masking failed test - expected 8 got {np.sum(cube.data)}" - ) - - def test_basic_cube_intersect_low_weight(self): - shape = shapely.geometry.box(0.1, 0.6, 1, 1) - masked_cube = mask_cube_from_shapefile( - self.basic_cube, shape, minimum_weight=0.2 - ) - assert np.sum(masked_cube.data) == 12, ( - f"basic cube masking weighting failed test - expected 12 got {np.sum(masked_cube.data)}" - ) - - def test_basic_cube_intersect_high_weight(self): - shape = shapely.geometry.box(0.1, 0.6, 1, 1) - masked_cube = mask_cube_from_shapefile( - self.basic_cube, shape, minimum_weight=0.7 - ) - assert np.sum(masked_cube.data) == 8, ( - f"basic cube masking weighting failed test- expected 8 got {np.sum(masked_cube.data)}" - ) - - def test_cube_list_error(self): - cubelist = iris.cube.CubeList([self.basic_cube]) - shape = shapely.geometry.box(1, 1, 2, 2) - with pytest.raises(TypeError, match="CubeList object rather than Cube"): - mask_cube_from_shapefile(cubelist, shape) - - def test_non_cube_error(self): - fake = None - shape = shapely.geometry.box(1, 1, 2, 2) - with pytest.raises(TypeError, match="Received non-Cube object"): - mask_cube_from_shapefile(fake, shape) - - def test_line_shape_warning(self): - shape = shapely.geometry.LineString([(0, 0.75), (2, 0.75)]) - with pytest.warns(IrisUserWarning, match="invalid type"): - masked_cube = mask_cube_from_shapefile( - self.basic_cube, shape, minimum_weight=0.1 - ) - assert np.sum(masked_cube.data) == 24, ( - f"basic cube masking against line failed test - expected 24 got {np.sum(masked_cube.data)}" - ) - - def test_cube_coord_mismatch_warning(self): - shape = shapely.geometry.box(0.6, 0.6, 0.9, 0.9) - cube = self.basic_cube - cube.coord("projection_x_coordinate").points = [180, 360, 540] - cube.coord("projection_x_coordinate").coord_system = RotatedGeogCS(30, 30) - with pytest.warns(IrisUserWarning, match="masking"): - mask_cube_from_shapefile( - cube, - shape, - ) - - def test_missing_xy_coord(self): - shape = shapely.geometry.box(0.6, 0.6, 0.9, 0.9) - cube = self.basic_cube - cube.remove_coord("projection_x_coordinate") - with pytest.raises(ValueError, match="1d xy coordinates"): - mask_cube_from_shapefile(cube, shape) - - def test_shape_not_shape(self): - shape = [5, 6, 7, 8] # random array - with pytest.raises(TypeError, match="valid Shapely"): - mask_cube_from_shapefile(self.basic_cube, shape) - - def test_shape_invalid(self): - shape = shapely.box(0, 1, 1, 1) - with pytest.raises(TypeError, match="valid Shapely"): - mask_cube_from_shapefile(self.basic_cube, shape) +@pytest.fixture +def mock_cube(): + """Create a mock 10x10 Iris cube for testing.""" + x_points = np.linspace(0, 9, 10) + y_points = np.linspace(0, 9, 10) + x_coord = DimCoord( + x_points, + standard_name="longitude", + units="degrees", + coord_system=GeogCS(6371229), + ) + y_coord = DimCoord( + y_points, + standard_name="latitude", + units="degrees", + coord_system=GeogCS(6371229), + ) + data = np.ones((len(y_points), len(x_points))) + cube = Cube(data, dim_coords_and_dims=[(y_coord, 0), (x_coord, 1)]) + return cube + + +def test_mask_cube_from_shapefile_inplace( + mock_cube, +): + shape = box(0, 0, 10, 10) + masked_cube = mask_cube_from_shapefile( + mock_cube, shape, shape_crs=CRS.from_epsg(4326), in_place=True + ) + assert masked_cube is None + + +def test_mask_cube_from_shapefile_not_inplace(mock_cube): + shape = box(0, 0, 10, 10) + masked_cube = mask_cube_from_shapefile( + mock_cube, shape, shape_crs=CRS.from_epsg(4326), in_place=False + ) + assert masked_cube is not None diff --git a/lib/iris/util.py b/lib/iris/util.py index 1795c8d87b..5723d44857 100644 --- a/lib/iris/util.py +++ b/lib/iris/util.py @@ -10,6 +10,7 @@ from collections.abc import Hashable, Iterable from copy import deepcopy import functools +from functools import wraps import inspect import os import os.path @@ -18,10 +19,13 @@ from typing import TYPE_CHECKING, List, Literal from warnings import warn +import cartopy import cf_units from dask import array as da import numpy as np import numpy.ma as ma +import pyproj +import shapely from iris._deprecation import warn_deprecated from iris._lazy_data import is_lazy_data, is_lazy_masked_data @@ -2189,27 +2193,94 @@ def _strip_metadata_from_dims(cube, dims): return reduced_cube -def mask_cube_from_shapefile(cube, shape, minimum_weight=0.0, in_place=False): - """Take a shape object and masks all points not touching it in a cube. +# def mutually_exclusive_keywords(keyword, *keywords): +# """Ensure that exactly one of the given keywords is specified. - 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 +# Used to ensure that only one of a set of mutually exclusive keyword arguments +# is provided to a function. + +# Parameters +# ---------- +# *keywords : str +# Names of mutually exclusive keyword arguments. + +# Returns +# ------- +# decorator +# A decorator that enforces the mutually exclusive constraint. + +# Raises +# ------ +# iris.exceptions.IrisError +# If zero or more than one of the specified keywords are present in kwargs. +# """ +# keywords = (keyword,) + keywords + +# def wrapper(func): +# @wraps(func) +# def inner(*args, **kwargs): +# if sum(k in keywords for k in kwargs) != 1: +# raise iris.exceptions.IrisError( +# "You must specify exactly one of {}".format(", ".join(keywords)) +# ) +# return func(*args, **kwargs) + +# return inner + +# return wrapper + + +# @mutually_exclusive_keywords("minimum_weight", "all_touched") +def mask_cube_from_shapefile( + cube: iris.cube.Cube, + shape: shapely.Geometry, + shape_crs: cartopy.crs | pyproj.CRS, + in_place: bool = False, + minimum_weight: float = 0.0, + **kwargs, +): + """Mask all points in a cube that do not intersect a shapefile object. + + Mask a :class:`~iris.cube.Cube` with any shapefile object, (e.g. Natural Earth Shapefiles via cartopy). + Finds the overlap between the `shapefile` shape and the :class:`~iris.cube.Cube` and + masks out any cells that __do not__ intersect the shape. + + Shapes can be polygons, lines or points. + + By default, all cells touched by geometries are kept (equivalent to `minimum_weight=0`). This behaviour + can be changed by increasing the `minimum_weight` keyword argument or setting `all_touched=False`, + then only the only cells whose center is within the polygon or that are selected by Bresenham’s line algorithm + (for line type shapes) are kept. For points, the `minimum_weight` is ignored, and the cell that intersects the point + is kept. Parameters ---------- cube : :class:`~iris.cube.Cube` object The `Cube` object to masked. Must be singular, rather than a `CubeList`. - shape : Shapely.Geometry object + shape : shapely.Geometry object A single `shape` of the area to remain unmasked on the `cube`. If it a line object of some kind then minimum_weight will be ignored, because you cannot compare the area of a 1D line and 2D Cell. - minimum_weight : float , default=0.0 - A number between 0-1 describing what % of a cube cell area must - the shape overlap to include it. + shape_crs : cartopy.crs.CRS, default=None + The coordinate reference system of the shape object. in_place : bool, default=False Whether to mask the `cube` in-place or return a newly masked `cube`. Defaults to False. + minimum_weight : float, default=0.0 + A number between 0-1 describing what % of a cube cell area must the shape overlap to be masked. + Only applied to polygon shapes. If the shape is a line or point then this is ignored. + + Other Parameters + ---------------- + all_touched : bool, default=None + If True, all cells touched by the shape are kept. If False, only cells whose + center is within the polygon or that are selected by Bresenham’s line algorithm + (for line type shape) are kept. + invert : bool, default=False + If True, the mask is inverted, meaning that cells that intersect the shape are masked out + and cells that do not intersect the shape are kept. If False, the mask is applied normally, + meaning that cells that intersect the shape are kept and cells that do not intersect the shape + are masked out. Returns ------- @@ -2221,28 +2292,65 @@ def mask_cube_from_shapefile(cube, shape, minimum_weight=0.0, in_place=False): :func:`~iris.util.mask_cube` Mask any cells in the cube’s data array. - 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 - is compared to the shapefile to determine overlap and populate a true/false array - This array is then used to mask the cube using the `iris.util.mask_cube` function - This uses numpy arithmetic logic for broadcasting, so you may encounter unexpected - results if your cube has other dimensions the same length as the x/y dimensions - Examples -------- >>> import shapely + >>> from pyproj import CRS >>> from iris.util import mask_cube_from_shapefile + + Extract a rectangular region covering the UK from a stereographic projection cube: + + >>> cube = iris.load_cube(iris.sample_data_path("toa_brightness_stereographic.nc")) + >>> shape = shapely.geometry.box(-10, 50, 2, 60) # box around the UK + >>> wgs84 = CRS.from_epsg(4326) # WGS84 coordinate system + >>> masked_cube = mask_cube_from_shapefile(cube, shape, wgs84) + + Extract a trajectory by using a line shapefile: + + >>> from shapely import LineString + >>> line = LineString([(-45, 40), (-28, 53), (-2, 55), (19, 45)]) + >>> masked_cube = mask_cube_from_shapefile(cube, line, wgs84) + + Standard shapely manipulations can be applied. For example, to extract a trajectory + with a 1 degree buffer around it: + + >>> buffer = line.buffer(1) + >>> masked_cube = mask_cube_from_shapefile(cube, buffer, wgs84) + + You can load more complex shapes from other libraries. For example, to extract the + Canadian provience of Ontario from a cube: + + >>> import cartopy.io.shapereader as shpreader + >>> admin1 = shpreader.natural_earth(resolution='110m', + category='cultural', + name='admin_1_states_provinces_lakes') + >>> admin1shp = shpreader.Reader(admin1).geometries() + >>> cube = iris.load_cube(iris.sample_data_path("E1_north_america.nc")) >>> shape = shapely.geometry.box(-100,30, -80,40) # box between 30N-40N 100W-80W - >>> masked_cube = mask_cube_from_shapefile(cube, shape) + >>> wgs84 = CRS.from_epsg(4326) + >>> masked_cube = mask_cube_from_shapefile(cube, shape, wgs84) + + Warning + ------- + For best masking results, both the cube _and_ masking geometry should have a + coordinate reference system (CRS) defined. Masking results will be most reliable + when the cube and masking geometry have the same CRS. + + If the cube has no coord_system, the default GeogCS is used where + the coordinate units are degrees. For any other coordinate units, + the cube **must** have a coord_system defined. + + If a CRS is not provided for the the masking geometry, the CRS of the cube is assumed. """ - shapefile_mask = create_shapefile_mask(shape, cube, minimum_weight) + shapefile_mask = create_shapefile_mask( + geometry=shape, + geometry_crs=shape_crs, + cube=cube, + minimum_weight=minimum_weight, + **kwargs, + ) masked_cube = mask_cube(cube, shapefile_mask, in_place=in_place) if not in_place: return masked_cube diff --git a/pyproject.toml b/pyproject.toml index ae99aaa554..72657a3536 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -52,6 +52,9 @@ license-files = ["LICENSE"] name = "scitools-iris" requires-python = ">=3.11" +[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 d14b9976bc..f5141337d0 100644 --- a/requirements/py312.yml +++ b/requirements/py312.yml @@ -34,6 +34,8 @@ dependencies: - pandas - pip - python-stratify + - rasterio + - affine # Test dependencies. - asv_runner diff --git a/requirements/py313.yml b/requirements/py313.yml index 1b6249b97b..61f7126bc7 100644 --- a/requirements/py313.yml +++ b/requirements/py313.yml @@ -34,6 +34,8 @@ dependencies: - pandas - pip - python-stratify + - raterio + - affine # Test dependencies. - asv_runner