Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 3 additions & 16 deletions lib/ants/io/load.py
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,6 @@
"""

import copy
import warnings
from contextlib import contextmanager
from functools import wraps

Expand Down Expand Up @@ -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)
Expand Down
6 changes: 2 additions & 4 deletions lib/ants/tests/fileformats/test_integration.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
145 changes: 136 additions & 9 deletions lib/ants/tests/utils/test_transform_bbox.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -46,24 +51,146 @@ 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)
self.assertEqual(len(res.geoms), 1)
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)
16 changes: 15 additions & 1 deletion lib/ants/utils/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -96,14 +96,28 @@ 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 "
"the target coordinate system limits ({})."
)
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


Expand Down