diff --git a/lib/iris/_shapefiles.py b/lib/iris/_shapefiles.py index 75b69921d9..6d79ac96a4 100644 --- a/lib/iris/_shapefiles.py +++ b/lib/iris/_shapefiles.py @@ -13,6 +13,8 @@ import shapely.ops import shapely.prepared as prepped +from iris.exceptions import IrisUserWarning + # ------------------------------------------------------------------------------ # GLOBAL VARIABLES # ------------------------------------------------------------------------------ @@ -37,9 +39,14 @@ def create_shapefile_mask( Returns: A numpy array of the shape of the x & y coordinates of the cube, with points to mask equal to True """ - # verification + from iris.cube import Cube, CubeList + try: + if geometry.is_valid is False: + raise TypeError("Geometry is not a valid Shapely object") + except Exception: + raise TypeError("Geometry is not a valid Shapely object") if not isinstance(cube, Cube): if isinstance(cube, CubeList): msg = "Received CubeList object rather than Cube - to mask a CubeList iterate over each Cube" @@ -61,7 +68,8 @@ def create_shapefile_mask( ): weights = False 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. " + "Shape is of invalid type for minimum weight masking, must use a Polygon rather than Line shape.\n Masking based off intersection instead. ", + IrisUserWarning, ) else: weights = True @@ -76,7 +84,6 @@ def create_shapefile_mask( _cube_xy_guessbounds(cube_2d) xmod = cube.coord(x_name).units.modulus ymod = cube.coord(y_name).units.modulus - cube.coord("longitude") mask_template = np.zeros(cube_2d.shape, dtype=np.float64) # perform the masking @@ -122,7 +129,10 @@ def _transform_coord_system(geometry, cube, geometry_system=None): DEFAULT_CS = _get_global_cs() target_system = cube.coord_system() if not target_system: - warnings.warn("Cube has no coord_system; using default GeogCS lat/lon") + warnings.warn( + "Cube has no coord_system; using default GeogCS lat/lon", + IrisUserWarning, + ) target_system = DEFAULT_CS if geometry_system is None: geometry_system = DEFAULT_CS diff --git a/lib/iris/tests/integration/test_shapefiles.py b/lib/iris/tests/integration/test_shapefiles.py new file mode 100644 index 0000000000..a9cd29336f --- /dev/null +++ b/lib/iris/tests/integration/test_shapefiles.py @@ -0,0 +1,79 @@ +import math + +import cartopy.io.shapereader as shpreader +import numpy as np + +import iris +import iris.tests as tests +from iris.util import apply_shapefile + + +@tests.skip_data +class TestCubeMasking(tests.IrisTest): + ne_countries = shpreader.natural_earth( + resolution="10m", category="cultural", name="admin_0_countries" + ) + reader = shpreader.Reader(ne_countries) + + def testGlobal(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 = apply_shapefile(ne_russia, test_global) + print(np.sum(masked_test.data)) + assert math.isclose( + np.sum(masked_test.data), 76845.37, rel_tol=0.00001 + ), "Global data with Russia mask failed test" + + def testRotated(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 = apply_shapefile(ne_germany, test_rotated) + assert math.isclose( + np.sum(masked_test.data), 179.46872, rel_tol=0.00001 + ), "rotated europe data with German mask failed test" + + def testTransverseMercator(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 = apply_shapefile(ne_uk, test_transverse) + assert math.isclose( + np.sum(masked_test.data), 90740.25, rel_tol=0.00001 + ), "transverse mercator UK data with UK mask failed test" + + def testRotatedWeighted(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 = apply_shapefile( + ne_germany, test_rotated, minimum_weight=0.9 + ) + assert math.isclose( + np.sum(masked_test.data), 125.60199, rel_tol=0.00001 + ), "rotated europe data with 0.9 weight germany mask failed test" diff --git a/lib/iris/tests/test_shapefiles.py b/lib/iris/tests/test_shapefiles.py index a9cd29336f..460488a298 100644 --- a/lib/iris/tests/test_shapefiles.py +++ b/lib/iris/tests/test_shapefiles.py @@ -1,79 +1,82 @@ -import math - -import cartopy.io.shapereader as shpreader import numpy as np +import pytest +import shapely -import iris +from iris.coords import DimCoord +import iris.cube +from iris.exceptions import IrisUserWarning import iris.tests as tests from iris.util import apply_shapefile -@tests.skip_data -class TestCubeMasking(tests.IrisTest): - ne_countries = shpreader.natural_earth( - resolution="10m", category="cultural", name="admin_0_countries" +class TestBasicCubeMasking(tests.IrisTest): + basic_data = np.array([[1, 2], [3, 4]]) + basic_cube = iris.cube.Cube(basic_data) + coord = DimCoord( + np.array([0, 1]), + standard_name="projection_y_coordinate", + bounds=[[0, 0.5], [0.5, 1]], + units="1", + ) + basic_cube.add_dim_coord(coord, 0) + coord = DimCoord( + np.array([0, 1]), + standard_name="projection_x_coordinate", + bounds=[[0, 0.5], [0.5, 1]], + units="1", ) - reader = shpreader.Reader(ne_countries) + basic_cube.add_dim_coord(coord, 1) - def testGlobal(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 = apply_shapefile(ne_russia, test_global) - print(np.sum(masked_test.data)) - assert math.isclose( - np.sum(masked_test.data), 76845.37, rel_tol=0.00001 - ), "Global data with Russia mask failed test" + def testBasicCubeIntersect(self): + shape = shapely.geometry.box(0.6, 0.6, 1, 1) + masked_cube = apply_shapefile(shape, self.basic_cube) + assert ( + np.sum(masked_cube.data) == 4 + ), f"basic cube masking failed test - expected 4 got {np.sum(masked_cube.data)}" - def testRotated(self): - path = tests.get_data_path( - ["NetCDF", "rotated", "xy", "rotPole_landAreaFraction.nc"] + def testBasicCubeIntersectLowWeight(self): + shape = shapely.geometry.box(0.5, 0.5, 1, 1) + masked_cube = apply_shapefile( + shape, self.basic_cube, minimum_weight=0.2 ) - 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 = apply_shapefile(ne_germany, test_rotated) - assert math.isclose( - np.sum(masked_test.data), 179.46872, rel_tol=0.00001 - ), "rotated europe data with German mask failed test" + assert ( + np.sum(masked_cube.data) == 4 + ), f"basic cube masking weighting failed test - expected 4 got {np.sum(masked_cube.data)}" - def testTransverseMercator(self): - path = tests.get_data_path( - ["NetCDF", "transverse_mercator", "tmean_1910_1910.nc"] + def testBasicCubeIntersectHighWeight(self): + shape = shapely.geometry.box(0.1, 0.6, 1, 1) + masked_cube = apply_shapefile( + shape, self.basic_cube, minimum_weight=0.5 ) - 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 = apply_shapefile(ne_uk, test_transverse) - assert math.isclose( - np.sum(masked_test.data), 90740.25, rel_tol=0.00001 - ), "transverse mercator UK data with UK mask failed test" + assert ( + np.sum(masked_cube.data) == 7 + ), f"basic cube masking weighting failed test- expected 7 got {np.sum(masked_cube.data)}" - def testRotatedWeighted(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 = apply_shapefile( - ne_germany, test_rotated, minimum_weight=0.9 - ) - assert math.isclose( - np.sum(masked_test.data), 125.60199, rel_tol=0.00001 - ), "rotated europe data with 0.9 weight germany mask failed test" + def testCubeListError(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" + ): + apply_shapefile(shape, cubelist) + + def testNonCubeError(self): + fake = None + shape = shapely.geometry.box(1, 1, 2, 2) + with pytest.raises(TypeError, match="Received non-Cube object"): + apply_shapefile(shape, fake) + + def testLineShapeWarning(self): + shape = shapely.geometry.LineString([(0, 0.75), (2, 0.75)]) + with pytest.warns(IrisUserWarning, match="invalid type"): + masked_cube = apply_shapefile( + shape, self.basic_cube, minimum_weight=0.1 + ) + assert ( + np.sum(masked_cube.data) == 7 + ), f"basic cube masking against line failed test - expected 7 got {np.sum(masked_cube.data)}" + + def testShapeInvalid(self): + shape = None + with pytest.raises(TypeError): + apply_shapefile(shape, self.basic_cube) diff --git a/lib/iris/util.py b/lib/iris/util.py index 7b6cf553fa..70b3d59acd 100644 --- a/lib/iris/util.py +++ b/lib/iris/util.py @@ -2190,12 +2190,6 @@ def apply_shapefile(shape, cube, minimum_weight=0.0, in_place=False): >>> ... """ - try: - if shape.is_valid is False: - raise TypeError("Geometry is not a valid Shapely object") - except Exception: - raise TypeError("Geometry is not a valid Shapely object") - shapefile_mask = create_shapefile_mask(shape, cube, minimum_weight) masked_cube = mask_cube(cube, shapefile_mask, in_place=in_place) return masked_cube