Skip to content

Commit

Permalink
Added more tests and split into integration and unit tests. Testing w…
Browse files Browse the repository at this point in the history
…arnings/errors aswell.
  • Loading branch information
acchamber committed Nov 9, 2023
1 parent 6decedf commit d168f24
Show file tree
Hide file tree
Showing 4 changed files with 162 additions and 76 deletions.
18 changes: 14 additions & 4 deletions lib/iris/_shapefiles.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,8 @@
import shapely.ops
import shapely.prepared as prepped

from iris.exceptions import IrisUserWarning

# ------------------------------------------------------------------------------
# GLOBAL VARIABLES
# ------------------------------------------------------------------------------
Expand All @@ -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"
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
79 changes: 79 additions & 0 deletions lib/iris/tests/integration/test_shapefiles.py
Original file line number Diff line number Diff line change
@@ -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"
135 changes: 69 additions & 66 deletions lib/iris/tests/test_shapefiles.py
Original file line number Diff line number Diff line change
@@ -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)
6 changes: 0 additions & 6 deletions lib/iris/util.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

0 comments on commit d168f24

Please sign in to comment.