diff --git a/lib/ants/io/load.py b/lib/ants/io/load.py index 089e828..33c2dc3 100644 --- a/lib/ants/io/load.py +++ b/lib/ants/io/load.py @@ -57,7 +57,6 @@ """ import copy -import warnings from contextlib import contextmanager from functools import wraps @@ -344,21 +343,9 @@ def load_function(*args, **kwargs): # Ensure that we leave appropriate calling to the underlying iris load # function. - # TODO https://github.com/MetOffice/ANTS/issues/91, remove warning filter - # workaround when iris issue https://github.com/SciTools/iris/issues/5749 has - # been fixed. - with warnings.catch_warnings(): - warnings.filterwarnings( - "ignore", - "Ignoring a datum in netCDF load for consistency with existing " - "behaviour. In a future version of Iris, this datum will be applied. " - "To apply the datum when loading, use the " - "iris.FUTURE.datum_support flag.", - FutureWarning, - ) - # Use context manager to avoid permanently modifying iris behaviour. - with ants_format_agent(): - cubes = func(*args, **kwargs) + # Use context manager to avoid permanently modifying iris behaviour. + with ants_format_agent(): + cubes = func(*args, **kwargs) if cubes is not None: try: ants.utils.cube.derive_circular_status(cubes) diff --git a/lib/ants/tests/fileformats/test_integration.py b/lib/ants/tests/fileformats/test_integration.py index 1ca549c..1cc522a 100644 --- a/lib/ants/tests/fileformats/test_integration.py +++ b/lib/ants/tests/fileformats/test_integration.py @@ -5,7 +5,6 @@ import tempfile import unittest.mock as mock import warnings -from unittest import expectedFailure import ants import ants.io.save as save @@ -106,13 +105,12 @@ def test_grib(self): class TestExceptions(ants.tests.TestCase): - @expectedFailure def test_no_datum_warning(self): # Save and reload a cube without a datum. - # If this test passes unexpectedly, we can complete ticket + # This test relates to the completed ticket # https://github.com/MetOffice/ANTS/issues/91. - # Currently, iris raises a warning even if the source file does not + # Previously, iris raised a warning even if the source file did not # contain a datum. See https://github.com/SciTools/iris/issues/5749. cube = ants.tests.stock.geodetic((2, 2)) assert cube.coord_system().datum is None diff --git a/lib/ants/tests/utils/test_transform_bbox.py b/lib/ants/tests/utils/test_transform_bbox.py index bd7b177..20ec072 100644 --- a/lib/ants/tests/utils/test_transform_bbox.py +++ b/lib/ants/tests/utils/test_transform_bbox.py @@ -2,7 +2,12 @@ # # This file is part of ANTS and is released under the BSD 3-Clause license. # See LICENSE.txt in the root of the repository for full licensing details. +import re + import ants.tests +import iris.coord_systems +import numpy as np +import pytest from ants.coord_systems import OSGB, UM_SPHERE from ants.utils import transform_bbox @@ -46,9 +51,62 @@ def test_points_crossing_dateline_range2(self): self.assertArrayAlmostEqual(res.bounds, current_target) +@pytest.mark.parametrize( + "lon_lat_shift", + [ + [45, 15.0], + [90.0, 30.0], + [135.0, 45.0], + [180.0, 90.0], + ], + ids=["(45.0, 15.0)", "(90.0, 30.0)", "(135.0, 45.0)", "(180.0, 90.0)"], +) +def test_um_sphere_remains_invertible(lon_lat_shift): + """Test that polygons remain invertible for two geodetic crss.""" + + lon_shift, lat_shift = lon_lat_shift + # origin of the OSGB in lat, lon + lat_0, lon_0 = 0, 0 + + bbox_points = np.array( + [lon_0 - lon_shift, lat_0 - lat_shift, lon_0 + lon_shift, lat_0 + lat_shift] + ) + bbox = [ + (lon_0 - lon_shift, lat_0 - lat_shift), + (lon_0 + lon_shift, lat_0 - lat_shift), + (lon_0 + lon_shift, lat_0 + lat_shift), + (lon_0 - lon_shift, lat_0 + lat_shift), + ] + + geoms = transform_bbox(bbox, UM_SPHERE.crs, UM_SPHERE.crs) + bounds = geoms.bounds + assert np.array_equal(bbox_points, bounds) + + class TestDiffCS(TestCommon, ants.tests.TestCase): # Different coordinate system tests. + def setUp(self): + self.msg = ( + "Attempting to project bounding box (GeogCS(6371229.0)) beyond " + "the extent of the target coordinate system limits " + "(TransverseMercator(.*" + ) + self.msg = self.msg.replace("(", r"\(") + self.msg = self.msg.replace(")", r"\)") + super().setUp() + def test_points_inside_projected_crs(self): + """Project a bounding box from the OSGB to UM Sphere. + + Notes + ----- + The ants OSGB crs is a general transverse mercator crs in iris + which has different projection limits to the cartopy OSGB crs. The + projection limits for the ants OSGB are the general limits of a transverse + mercator and are larger than the cartopy OSGB crs which is restricted to a + valid domain over the UK. + """ + bbox_points = (-12, -12, 7e5, 13e5) bbox = self._gen_bbox(*bbox_points) res = transform_bbox(bbox, OSGB.crs, UM_SPHERE.crs) @@ -56,14 +114,83 @@ def test_points_inside_projected_crs(self): tar = [-9.49660933, 49.76607039, 3.63474423, 61.46518886] self.assertArrayAlmostEqual(res.bounds, tar) - def test_point_lie_beyond_crs_definition(self): - bbox = [(-180, -90), (180, -90), (180, 90), (-180, 90)] - msg = ( - "Attempting to project bounding box (GeogCS(6371229.0)) beyond " - "the extent of the target coordinate system limits " - "(TransverseMercator(.*" + def test_ants_osgb_not_invertible_global_poly(self): + """Test that a polygon spanning the globe is not invertible. + + The generic Cartopy Transverse Mercator CRS uses fixed + projection-domain bounds of approximately + (-2e7, -1e7, 2e7, 1e7) metres. + """ + + # origin of the OSGB in lat, lon + lat_0, lon_0 = 0, 0 + inv_msg = re.escape( + "The bounding box in the crs (TransverseMercator) is not invertible" + " to the source crs (GeogCS)." ) - msg = msg.replace("(", r"\(") - msg = msg.replace(")", r"\)") - with self.assertRaisesRegex(ValueError, msg): + bbox_points = (lon_0 - 180.0, lat_0 - 90, lon_0 + 180, lat_0 + 90) + bbox = self._gen_bbox(*bbox_points) + with self.assertRaisesRegex(ValueError, inv_msg): transform_bbox(bbox, UM_SPHERE.crs, OSGB.crs) + + def test_point_lie_beyond_crs_definition(self): + """Test invalid projection with the iris OSGB. + + As the OSGB crs is a regional crs (transverse Mercator), we only + get accurate projections within a restricted domain. The iris OSGB + returns the cartopy OSGB crs when converted to a cartopy projection. + The projection limits of this domain are (0, 0, 7e5, 13e5). + + The bounds are not clipped in this case and instead return NaN. + """ + + osgb_crs = iris.coord_systems.OSGB() + bbox_points = (-180, -90, 180, 90) + bbox = self._gen_bbox(*bbox_points) + + with self.assertRaisesRegex(ValueError, self.msg): + transform_bbox(bbox, UM_SPHERE.crs, osgb_crs) + + def test_nan_bounds(self): + """Demonstrate some out-of-domain projections produce NaN bounds + rather than clipped bounds. + + In this case we define a box which which produces bounds of NaN. Once + we are far from sensible projection limits, the behaviour becomes + inconsistent. + """ + # origin of the OSGB in lat, lon + lat_0, lon_0 = 49.0, -2.0 + + bbox_points = (lon_0 + 120, lat_0 - 70, lon_0 + 140, lat_0 - 50) + bbox = self._gen_bbox(*bbox_points) + + with self.assertRaisesRegex(ValueError, self.msg): + transform_bbox(bbox, UM_SPHERE.crs, OSGB.crs) + + +@pytest.mark.parametrize("lon_shift", [85.0, 90.0, 145.0, 180.0]) +def test_ants_osgb_not_invertible(lon_shift): + """Test polygons that are clipped to the domain boundary are + not invertible. + + For large enough polygons the bounds of the polygon become + the bounds of the domain at (-2e7, -1e7, 2e7, 1e7) metres. + """ + + # origin of the OSGB in lat, lon + lat_0, lon_0 = 49.0, -2.0 + + inv_msg = re.escape( + "The bounding box in the crs (TransverseMercator) is not invertible" + " to the source crs (GeogCS)." + ) + + bbox = [ + (lon_0 - lon_shift, lat_0 - 85), + (lon_0 + lon_shift, lat_0 - 85), + (lon_0 + lon_shift, lat_0 + 85), + (lon_0 - lon_shift, lat_0 + 85), + ] + with pytest.raises(ValueError, match=inv_msg): + transform_bbox(bbox, UM_SPHERE.crs, OSGB.crs) diff --git a/lib/ants/utils/__init__.py b/lib/ants/utils/__init__.py index aab9230..02cfa5c 100644 --- a/lib/ants/utils/__init__.py +++ b/lib/ants/utils/__init__.py @@ -96,7 +96,8 @@ def transform_bbox(points, src_crs, tgt_crs): src_geom = sgeom.Polygon(points) tgt_geom = cartopy_tgt_crs.project_geometry(src_geom, cartopy_src_crs) - for item in tgt_geom.bounds: + bounds = tgt_geom.bounds + for item in bounds: if np.isnan(item): msg = ( "Attempting to project bounding box ({}) beyond the extent of " @@ -104,6 +105,19 @@ def transform_bbox(points, src_crs, tgt_crs): ) raise ValueError(msg.format(src_crs, tgt_crs)) + # Check the bounds are invertible. + inversion = cartopy_src_crs.transform_points( + cartopy_tgt_crs, + np.array([bounds[0], bounds[2]]), + np.array([bounds[1], bounds[3]]), + ) + if np.isnan(inversion).any(): + msg = ( + "The bounding box in the crs ({}) is not invertible" + " to the source crs ({})." + ) + raise ValueError(msg.format(type(tgt_crs).__name__, type(src_crs).__name__)) + return tgt_geom