diff --git a/lib/iris/tests/integration/test_trajectory.py b/lib/iris/tests/integration/test_trajectory.py index 6c54d15a42..bc38a4d280 100644 --- a/lib/iris/tests/integration/test_trajectory.py +++ b/lib/iris/tests/integration/test_trajectory.py @@ -27,41 +27,11 @@ import biggus import numpy as np -import iris.analysis.trajectory -import iris.tests.stock +import iris +import iris.tests.stock as istk - -@tests.skip_data -class TestSimple(tests.IrisTest): - def test_invalid_coord(self): - cube = iris.tests.stock.realistic_4d() - sample_pts = [('altitude', [0, 10, 50])] - with self.assertRaises(ValueError): - iris.analysis.trajectory.interpolate(cube, sample_pts, 'nearest') - - -class TestTrajectory(tests.IrisTest): - def test_trajectory_definition(self): - # basic 2-seg line along x - waypoints = [{'lat': 0, 'lon': 0}, {'lat': 0, 'lon': 1}, - {'lat': 0, 'lon': 2}] - trajectory = iris.analysis.trajectory.Trajectory(waypoints, - sample_count=21) - - self.assertEqual(trajectory.length, 2.0) - self.assertEqual(trajectory.sampled_points[19], - {'lat': 0.0, 'lon': 1.9000000000000001}) - - # 4-seg m-shape - waypoints = [{'lat': 0, 'lon': 0}, {'lat': 1, 'lon': 1}, - {'lat': 0, 'lon': 2}, {'lat': 1, 'lon': 3}, - {'lat': 0, 'lon': 4}] - trajectory = iris.analysis.trajectory.Trajectory(waypoints, - sample_count=33) - - self.assertEqual(trajectory.length, 5.6568542494923806) - self.assertEqual(trajectory.sampled_points[31], - {'lat': 0.12499999999999989, 'lon': 3.875}) +from iris.analysis.trajectory import (Trajectory, + interpolate as traj_interpolate) @tests.skip_data @@ -80,11 +50,10 @@ def setUp(self): def test_trajectory_extraction(self): # Pull out a single point - no interpolation required - single_point = iris.analysis.trajectory.interpolate( - self.cube, [('grid_latitude', [-0.1188]), - ('grid_longitude', [359.57958984])]) + single_point = traj_interpolate(self.cube, + [('grid_latitude', [-0.1188]), + ('grid_longitude', [359.57958984])]) expected = self.cube[..., 10, 0].data - self.assertArrayAllClose(single_point[..., 0].data, expected, rtol=2.0e-7) self.assertCML(single_point, ('trajectory', 'single_point.cml'), @@ -100,8 +69,7 @@ def test_trajectory_extraction_calc(self): y0 = scube.data[0, 0] y1 = scube.data[0, 1] expected = y0 + ((y1 - y0) * ((359.584090412 - x0)/(x1 - x0))) - trajectory_cube = iris.analysis.trajectory.interpolate(scube, - single_point) + trajectory_cube = traj_interpolate(scube, single_point) self.assertArrayAllClose(trajectory_cube.data, expected, rtol=2.0e-7) def _traj_to_sample_points(self, trajectory): @@ -121,12 +89,9 @@ def test_trajectory_extraction_axis_aligned(self): {'grid_latitude': -0.1188, 'grid_longitude': 359.57958984}, {'grid_latitude': -0.1188, 'grid_longitude': 359.66870117} ] - trajectory = iris.analysis.trajectory.Trajectory(waypoints, - sample_count=100) - + trajectory = Trajectory(waypoints, sample_count=100) sample_points = self._traj_to_sample_points(trajectory) - trajectory_cube = iris.analysis.trajectory.interpolate(self.cube, - sample_points) + trajectory_cube = traj_interpolate(self.cube, sample_points) self.assertCML(trajectory_cube, ('trajectory', 'constant_latitude.cml')) @@ -137,10 +102,9 @@ def test_trajectory_extraction_zigzag(self): {'grid_latitude': -0.0828, 'grid_longitude': 359.6606}, {'grid_latitude': -0.0468, 'grid_longitude': 359.6246}, ] - trajectory = iris.analysis.trajectory.Trajectory(waypoints, - sample_count=20) + trajectory = Trajectory(waypoints, sample_count=20) sample_points = self._traj_to_sample_points(trajectory) - trajectory_cube = iris.analysis.trajectory.interpolate( + trajectory_cube = traj_interpolate( self.cube[0, 0], sample_points) expected = np.array([287.95953369, 287.9190979, 287.95550537, 287.93240356, 287.83850098, 287.87869263, @@ -154,6 +118,39 @@ def test_trajectory_extraction_zigzag(self): checksum=False) self.assertArrayAllClose(trajectory_cube.data, expected, rtol=2.0e-7) + def test_colpex__nearest(self): + # Check a smallish nearest-neighbour interpolation against a result + # snapshot. + test_cube = self.cube[0][0] + # Test points on a regular grid, a bit larger than the source region. + xmin, xmax = [fn(test_cube.coord(axis='x').points) + for fn in (np.min, np.max)] + ymin, ymax = [fn(test_cube.coord(axis='x').points) + for fn in (np.min, np.max)] + fractions = [-0.23, -0.01, 0.27, 0.624, 0.983, 1.052, 1.43] + x_points = [xmin + frac * (xmax - xmin) for frac in fractions] + y_points = [ymin + frac * (ymax - ymin) for frac in fractions] + x_points, y_points = np.meshgrid(x_points, y_points) + sample_points = [('grid_longitude', x_points.flatten()), + ('grid_latitude', y_points.flatten())] + result = traj_interpolate(test_cube, sample_points, + method='nearest') + expected = [ + 288.07168579, 288.07168579, 287.9367981, 287.82736206, + 287.78564453, 287.8374939, 287.8374939, 288.07168579, + 288.07168579, 287.9367981, 287.82736206, 287.78564453, + 287.8374939, 287.8374939, 288.07168579, 288.07168579, + 287.9367981, 287.82736206, 287.78564453, 287.8374939, + 287.8374939, 288.07168579, 288.07168579, 287.9367981, + 287.82736206, 287.78564453, 287.8374939, 287.8374939, + 288.07168579, 288.07168579, 287.9367981, 287.82736206, + 287.78564453, 287.8374939, 287.8374939, 288.07168579, + 288.07168579, 287.9367981, 287.82736206, 287.78564453, + 287.8374939, 287.8374939, 288.07168579, 288.07168579, + 287.9367981, 287.82736206, 287.78564453, 287.8374939, + 287.8374939] + self.assertArrayAllClose(result.data, expected) + @tests.skip_data class TestTriPolar(tests.IrisTest): @@ -176,10 +173,8 @@ def setUp(self): ('latitude', latitudes)] def test_tri_polar(self): - # extract - sampled_cube = iris.analysis.trajectory.interpolate(self.cube, - self.sample_points) + sampled_cube = traj_interpolate(self.cube, self.sample_points) self.assertCML(sampled_cube, ('trajectory', 'tri_polar_latitude_slice.cml')) @@ -187,25 +182,62 @@ def test_tri_polar_method_linear_fails(self): # Try to request linear interpolation. # Not allowed, as we have multi-dimensional coords. self.assertRaises(iris.exceptions.CoordinateMultiDimError, - iris.analysis.trajectory.interpolate, self.cube, + traj_interpolate, self.cube, self.sample_points, method="linear") def test_tri_polar_method_unknown_fails(self): # Try to request unknown interpolation. - self.assertRaises(ValueError, iris.analysis.trajectory.interpolate, + self.assertRaises(ValueError, traj_interpolate, self.cube, self.sample_points, method="linekar") - -class TestHybridHeight(tests.IrisTest): + def test_tri_polar__nearest(self): + # Check a smallish nearest-neighbour interpolation against a result + # snapshot. + test_cube = self.cube + # Use just one 2d layer, just to be faster. + test_cube = test_cube[0][0] + # Fix the fill value of the data to zero, just so that we get the same + # result under numpy < 1.11 as with 1.11. + # NOTE: numpy<1.11 *used* to assign missing data points into an + # unmasked array as =0.0, now =fill-value. + # TODO: arguably, we should support masked data properly in the + # interpolation routine. In the legacy code, that is unfortunately + # just not the case. + test_cube.data.set_fill_value(0.0) + + # Test points on a regular global grid, with unrelated steps + offsets + # and an extended range of longitude values. + x_points = np.arange(-185.23, +360.0, 73.123) + y_points = np.arange(-89.12, +90.0, 42.847) + x_points, y_points = np.meshgrid(x_points, y_points) + sample_points = [('longitude', x_points.flatten()), + ('latitude', y_points.flatten())] + result = traj_interpolate(test_cube, sample_points, + method='nearest') + expected = [ + 0., 0., 0., 0., + 0., 0., 0., 0., + 12.13186264, 10.69991493, 9.86881161, 7.08723927, + 9.04308414, 12.56258678, 10.63761806, 9.19426727, + 28.93525505, 23.85289955, 26.94649506, 0., + 27.88831711, 28.65439224, 23.39414215, 26.78363228, + 13.53453922, 0., 17.41485596, 0., + 0., 13.0413475, 0., 17.10849571, + -1.67040622, -1.64783156, 0., -1.97898054, + -1.67642927, -1.65173221, -1.623945, 0.] + + self.assertArrayAllClose(result.data, expected) + + +class TestLazyData(tests.IrisTest): def test_hybrid_height(self): - cube = tests.stock.simple_4d_with_hybrid_height() + cube = istk.simple_4d_with_hybrid_height() # Put a biggus array on the cube so we can test deferred loading. cube.lazy_data(biggus.NumpyArrayAdapter(cube.data)) traj = (('grid_latitude', [20.5, 21.5, 22.5, 23.5]), ('grid_longitude', [31, 32, 33, 34])) - xsec = iris.analysis.trajectory.interpolate(cube, traj, - method='nearest') + xsec = traj_interpolate(cube, traj, method='nearest') # Check that creating the trajectory hasn't led to the original # data being loaded. diff --git a/lib/iris/tests/unit/analysis/interpolate_private/__init__.py b/lib/iris/tests/unit/analysis/interpolate_private/__init__.py new file mode 100644 index 0000000000..bb8f1ebc4a --- /dev/null +++ b/lib/iris/tests/unit/analysis/interpolate_private/__init__.py @@ -0,0 +1,20 @@ +# (C) British Crown Copyright 2016, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +"""Unit tests for the :mod:`iris.analysis._interpolate_private` module.""" + +from __future__ import (absolute_import, division, print_function) +from six.moves import (filter, input, map, range, zip) # noqa diff --git a/lib/iris/tests/unit/analysis/interpolate_private/test__nearest_neighbour_indices_ndcoords.py b/lib/iris/tests/unit/analysis/interpolate_private/test__nearest_neighbour_indices_ndcoords.py new file mode 100644 index 0000000000..6996a2bad5 --- /dev/null +++ b/lib/iris/tests/unit/analysis/interpolate_private/test__nearest_neighbour_indices_ndcoords.py @@ -0,0 +1,208 @@ +# (C) British Crown Copyright 2016, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +""" +Unit tests for +:meth:`iris.analysis._interpolate_private._nearest_neighbour_indices_ndcoords`. + +""" + +from __future__ import (absolute_import, division, print_function) +from six.moves import (filter, input, map, range, zip) # noqa + +# Import iris.tests first so that some things can be initialised before +# importing anything else. +import iris.tests as tests + +import mock +import numpy as np + +from iris.cube import Cube +from iris.coords import DimCoord, AuxCoord + +from iris.analysis._interpolate_private import \ + _nearest_neighbour_indices_ndcoords as nn_ndinds + + +class Test2d(tests.IrisTest): + def test_nonlatlon_simple_2d(self): + co_y = DimCoord([10.0, 20.0], long_name='y') + co_x = DimCoord([1.0, 2.0, 3.0], long_name='x') + cube = Cube(np.zeros((2, 3))) + cube.add_dim_coord(co_y, 0) + cube.add_dim_coord(co_x, 1) + sample_point = [('x', 2.8), ('y', 18.5)] + result = nn_ndinds(cube, sample_point) + self.assertEqual(result, (1, 2)) + + def test_latlon_simple_2d(self): + co_y = DimCoord([10.0, 20.0], + standard_name='latitude', units='degrees') + co_x = DimCoord([1.0, 2.0, 3.0], + standard_name='longitude', units='degrees') + cube = Cube(np.zeros((2, 3))) + cube.add_dim_coord(co_y, 0) + cube.add_dim_coord(co_x, 1) + sample_point = [('longitude', 2.8), ('latitude', 18.5)] + result = nn_ndinds(cube, sample_point) + self.assertEqual(result, (1, 2)) + + +class Test1d(tests.IrisTest): + def test_nonlatlon_simple_1d(self): + co_x = AuxCoord([1.0, 2.0, 3.0, 1.0, 2.0, 3.0], long_name='x') + co_y = AuxCoord([10.0, 10.0, 10.0, 20.0, 20.0, 20.0], long_name='y') + cube = Cube(np.zeros(6)) + cube.add_aux_coord(co_y, 0) + cube.add_aux_coord(co_x, 0) + sample_point = [('x', 2.8), ('y', 18.5)] + result = nn_ndinds(cube, sample_point) + self.assertEqual(result, (5,)) + + def test_latlon_simple_1d(self): + cube = Cube([11.0, 12.0, 13.0, 21.0, 22.0, 23.0]) + co_x = AuxCoord([1.0, 2.0, 3.0, 1.0, 2.0, 3.0], + standard_name='longitude', units='degrees') + co_y = AuxCoord([10.0, 10.0, 10.0, 20.0, 20.0, 20.0], + standard_name='latitude', units='degrees') + cube.add_aux_coord(co_y, 0) + cube.add_aux_coord(co_x, 0) + sample_point = [('longitude', 2.8), ('latitude', 18.5)] + result = nn_ndinds(cube, sample_point) + self.assertEqual(result, (5,)) + + +class TestApiExtras(tests.IrisTest): + # Check operation with alternative calling setups. + def test_no_y_dim(self): + # Operate in X only, returned slice should be [:, ix]. + co_x = DimCoord([1.0, 2.0, 3.0], long_name='x') + co_y = DimCoord([10.0, 20.0], long_name='y') + cube = Cube(np.zeros((2, 3))) + cube.add_dim_coord(co_y, 0) + cube.add_dim_coord(co_x, 1) + sample_point = [('x', 2.8)] + result = nn_ndinds(cube, sample_point) + self.assertEqual(result, (slice(None), 2)) + + def test_no_x_dim(self): + # Operate in Y only, returned slice should be [iy, :]. + co_x = DimCoord([1.0, 2.0, 3.0], long_name='x') + co_y = DimCoord([10.0, 20.0], long_name='y') + cube = Cube(np.zeros((2, 3))) + cube.add_dim_coord(co_y, 0) + cube.add_dim_coord(co_x, 1) + sample_point = [('y', 18.5)] + result = nn_ndinds(cube, sample_point) + self.assertEqual(result, (1, slice(None))) + + def test_sample_dictionary(self): + # Pass sample_point arg as a dictionary: this usage mode is deprecated. + co_x = AuxCoord([1.0, 2.0, 3.0], long_name='x') + co_y = AuxCoord([10.0, 20.0], long_name='y') + cube = Cube(np.zeros((2, 3))) + cube.add_aux_coord(co_y, 0) + cube.add_aux_coord(co_x, 1) + sample_point = {'x': 2.8, 'y': 18.5} + warn_call = self.patch( + 'iris.analysis._interpolate_private.warn_deprecated') + result = nn_ndinds(cube, sample_point) + self.assertEqual(result, (1, 2)) + self.assertEqual(warn_call.call_count, 1) + self.assertIn('dictionary to specify points is deprecated', + warn_call.call_args[0][0]) + + +class TestLatlon(tests.IrisTest): + # Check correct calculations on lat-lon points. + def _testcube_latlon_1d(self, lats, lons): + cube = Cube(np.zeros(len(lons))) + co_x = AuxCoord(lons, standard_name='longitude', units='degrees') + co_y = AuxCoord(lats, standard_name='latitude', units='degrees') + cube.add_aux_coord(co_y, 0) + cube.add_aux_coord(co_x, 0) + return cube + + def _check_latlon_1d(self, lats, lons, sample_point, expect): + cube = self._testcube_latlon_1d(lats, lons) + result = nn_ndinds(cube, sample_point) + self.assertEqual(result, (expect,)) + + def test_lat_scaling(self): + # Check that (88, 25) is closer to (88, 0) than to (87, 25) + self._check_latlon_1d( + lats=[88, 87], + lons=[0, 25], + sample_point=[('latitude', 88), ('longitude', 25)], + expect=0) + + def test_alternate_latlon_names_okay(self): + # Check that (88, 25) is **STILL** closer to (88, 0) than to (87, 25) + # ... when coords have odd, but still recognisable, latlon names. + cube = self._testcube_latlon_1d(lats=[88, 87], + lons=[0, 25]) + cube.coord('latitude').rename('y_latitude_y') + cube.coord('longitude').rename('x_longitude_x') + sample_point = [('y_latitude_y', 88), ('x_longitude_x', 25)] + result = nn_ndinds(cube, sample_point) + self.assertEqual(result, (0,)) + + def test_alternate_nonlatlon_names_different(self): + # Check that (88, 25) is **NOT** closer to (88, 0) than to (87, 25) + # ... by plain XY euclidean-distance, if coords have non-latlon names. + cube = self._testcube_latlon_1d(lats=[88, 87], + lons=[0, 25]) + cube.coord('latitude').rename('y') + cube.coord('longitude').rename('x') + sample_point = [('y', 88), ('x', 25)] + result = nn_ndinds(cube, sample_point) + self.assertEqual(result, (1,)) + + def test_lons_wrap_359_0(self): + # Check that (0, 359) is closer to (0, 0) than to (0, 350) + self._check_latlon_1d( + lats=[0, 0], + lons=[0, 350], + sample_point=[('latitude', 0), ('longitude', 359)], + expect=0) + + def test_lons_wrap_359_neg1(self): + # Check that (0, 359) is closer to (0, -1) than to (0, 350) + self._check_latlon_1d( + lats=[0, 0], + lons=[350, -1], + sample_point=[('latitude', 0), ('longitude', 359)], + expect=1) + + def test_lons_wrap_neg179_plus179(self): + # Check that (0, -179) is closer to (0, 179) than to (0, -170) + self._check_latlon_1d( + lats=[0, 0], + lons=[-170, 179], + sample_point=[('latitude', 0), ('longitude', -179)], + expect=1) + + def test_lons_over_pole(self): + # Check that (89, 0) is closer to (89, 180) than to (85, 0) + self._check_latlon_1d( + lats=[85, 89], + lons=[0, 180], + sample_point=[('latitude', 89), ('longitude', 0)], + expect=1) + + +if __name__ == "__main__": + tests.main() diff --git a/lib/iris/tests/unit/analysis/trajectory/__init__.py b/lib/iris/tests/unit/analysis/trajectory/__init__.py new file mode 100644 index 0000000000..da78cf3095 --- /dev/null +++ b/lib/iris/tests/unit/analysis/trajectory/__init__.py @@ -0,0 +1,20 @@ +# (C) British Crown Copyright 2016, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +"""Unit tests for the :mod:`iris.analysis.trajectory` module.""" + +from __future__ import (absolute_import, division, print_function) +from six.moves import (filter, input, map, range, zip) # noqa diff --git a/lib/iris/tests/unit/analysis/trajectory/test_Trajectory.py b/lib/iris/tests/unit/analysis/trajectory/test_Trajectory.py new file mode 100644 index 0000000000..85ca32e979 --- /dev/null +++ b/lib/iris/tests/unit/analysis/trajectory/test_Trajectory.py @@ -0,0 +1,74 @@ +# (C) British Crown Copyright 2016, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +""" +Unit tests for :class:`iris.analysis.trajectory.Trajectory`. + +""" + +from __future__ import (absolute_import, division, print_function) +from six.moves import (filter, input, map, range, zip) # noqa + +# Import iris.tests first so that some things can be initialised before +# importing anything else. +import iris.tests as tests + +import numpy as np + +from iris.analysis.trajectory import Trajectory + + +class Test___init__(tests.IrisTest): + def test_2_points(self): + # basic 2-seg line along x + waypoints = [{'lat': 0, 'lon': 0}, {'lat': 1, 'lon': 2}] + trajectory = Trajectory(waypoints, sample_count=5) + + self.assertEqual(trajectory.length, np.sqrt(5)) + self.assertEqual(trajectory.sample_count, 5) + self.assertEqual(trajectory.sampled_points, + [{'lat': 0.0, 'lon': 0.0}, + {'lat': 0.25, 'lon': 0.5}, + {'lat': 0.5, 'lon': 1.0}, + {'lat': 0.75, 'lon': 1.5}, + {'lat': 1.0, 'lon': 2.0}]) + + def test_3_points(self): + # basic 2-seg line along x + waypoints = [{'lat': 0, 'lon': 0}, {'lat': 0, 'lon': 1}, + {'lat': 0, 'lon': 2}] + trajectory = Trajectory(waypoints, sample_count=21) + + self.assertEqual(trajectory.length, 2.0) + self.assertEqual(trajectory.sample_count, 21) + self.assertEqual(trajectory.sampled_points[19], + {'lat': 0.0, 'lon': 1.9000000000000001}) + + def test_zigzag(self): + # 4-seg m-shape + waypoints = [{'lat': 0, 'lon': 0}, {'lat': 1, 'lon': 1}, + {'lat': 0, 'lon': 2}, {'lat': 1, 'lon': 3}, + {'lat': 0, 'lon': 4}] + trajectory = Trajectory(waypoints, sample_count=33) + + self.assertEqual(trajectory.length, 5.6568542494923806) + self.assertEqual(trajectory.sample_count, 33) + self.assertEqual(trajectory.sampled_points[31], + {'lat': 0.12499999999999989, 'lon': 3.875}) + + +if __name__ == "__main__": + tests.main() diff --git a/lib/iris/tests/unit/analysis/trajectory/test_interpolate.py b/lib/iris/tests/unit/analysis/trajectory/test_interpolate.py new file mode 100644 index 0000000000..b04d422b1e --- /dev/null +++ b/lib/iris/tests/unit/analysis/trajectory/test_interpolate.py @@ -0,0 +1,184 @@ +# (C) British Crown Copyright 2016, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +""" +Unit tests for :meth:`iris.analysis.trajectory.interpolate`. + +""" + +from __future__ import (absolute_import, division, print_function) +from six.moves import (filter, input, map, range, zip) # noqa + +# Import iris.tests first so that some things can be initialised before +# importing anything else. +import iris.tests as tests + +import numpy as np + +from iris.coords import AuxCoord, DimCoord +import iris.tests.stock + +from iris.analysis.trajectory import interpolate + + +class TestFailCases(tests.IrisTest): + @tests.skip_data + def test_derived_coord(self): + cube = iris.tests.stock.realistic_4d() + sample_pts = [('altitude', [0, 10, 50])] + msg = "'altitude'.*derived coordinates are not allowed" + with self.assertRaisesRegexp(ValueError, msg): + interpolate(cube, sample_pts, 'nearest') + + # Try to request unknown interpolation method. + def test_unknown_method(self): + cube = iris.tests.stock.simple_2d() + sample_point = [('x', 2.8)] + msg = 'Unhandled interpolation.*linekar' + with self.assertRaisesRegexp(ValueError, msg): + interpolate(cube, sample_point, method="linekar") + + +class TestNearest(tests.IrisTest): + # Test interpolation with 'nearest' method. + # This is basically a wrapper to the routine: + # 'analysis._interpolate_private._nearest_neighbour_indices_ndcoords'. + # That has its own test, so we don't test the basic calculation + # exhaustively here. Instead we check the way it handles the source and + # result cubes (especially coordinates). + + def setUp(self): + cube = iris.tests.stock.simple_3d() + # Actually, this cube *isn't* terribly realistic, as the lat+lon coords + # have integer type, which in this case produces some peculiar results. + # Let's fix that (and not bother to test the peculiar behaviour). + for coord_name in ('longitude', 'latitude'): + coord = cube.coord(coord_name) + coord.points = coord.points.astype(float) + self.test_cube = cube + # Define coordinates for a single-point testcase. + y_val, x_val = 0, -90 + # Work out cube indices of the testpoint. + self.single_point_iy = \ + np.where(cube.coord('latitude').points == y_val)[0][0] + self.single_point_ix = \ + np.where(cube.coord('longitude').points == x_val)[0][0] + # Use slightly-different values to test nearest-neighbour operation. + self.single_sample_point = [('latitude', [y_val + 19.23]), + ('longitude', [x_val - 17.54])] + + def test_single_point_same_cube(self): + # Check exact result matching for a single point. + cube = self.test_cube + result = interpolate(cube, self.single_sample_point, method='nearest') + # Check that the result is a single trajectory point, exactly equal to + # the expected part of the original data. + self.assertEqual(result.shape[-1], 1) + result = result[..., 0] + expected = cube[:, self.single_point_iy, self.single_point_ix] + self.assertEqual(result, expected) + + def test_multi_point_same_cube(self): + # Check an exact result for multiple points. + cube = self.test_cube + # Use latitude selection to recreate a whole row of the original cube. + sample_points = [('longitude', [-180, -90, 0, 90]), + ('latitude', [0, 0, 0, 0])] + result = interpolate(cube, sample_points, method='nearest') + + # The result should be identical to a single latitude section of the + # original, but with modified coords (latitude has 4 repeated zeros). + expected = cube[:, 1, :] + # Result 'longitude' is now an aux coord. + co_x = expected.coord('longitude') + expected.remove_coord(co_x) + expected.add_aux_coord(co_x, 1) + # Result 'latitude' is now an aux coord containing 4*[0]. + expected.remove_coord('latitude') + co_y = AuxCoord([0, 0, 0, 0], + standard_name='latitude', units='degrees') + expected.add_aux_coord(co_y, 1) + self.assertEqual(result, expected) + + def test_aux_coord_noninterpolation_dim(self): + # Check exact result with an aux-coord mapped to an uninterpolated dim. + cube = self.test_cube + cube.add_aux_coord(DimCoord([17, 19], long_name='aux0'), 0) + + # The result cube should exactly equal a single source point. + result = interpolate(cube, self.single_sample_point, method='nearest') + self.assertEqual(result.shape[-1], 1) + result = result[..., 0] + expected = cube[:, self.single_point_iy, self.single_point_ix] + self.assertEqual(result, expected) + + def test_aux_coord_one_interp_dim(self): + # Check exact result with an aux-coord over one interpolation dims. + cube = self.test_cube + cube.add_aux_coord(AuxCoord([11, 12, 13, 14], long_name='aux_x'), + 2) + + # The result cube should exactly equal a single source point. + result = interpolate(cube, self.single_sample_point, method='nearest') + self.assertEqual(result.shape[-1], 1) + result = result[..., 0] + expected = cube[:, self.single_point_iy, self.single_point_ix] + self.assertEqual(result, expected) + + def test_aux_coord_both_interp_dims(self): + # Check exact result with an aux-coord over both interpolation dims. + cube = self.test_cube + cube.add_aux_coord(AuxCoord([[11, 12, 13, 14], + [21, 22, 23, 24], + [31, 32, 33, 34]], long_name='aux_xy'), + (1, 2)) + + # The result cube should exactly equal a single source point. + result = interpolate(cube, self.single_sample_point, method='nearest') + self.assertEqual(result.shape[-1], 1) + result = result[..., 0] + expected = cube[:, self.single_point_iy, self.single_point_ix] + self.assertEqual(result, expected) + + def test_aux_coord_fail_mixed_dims(self): + # Check behaviour with an aux-coord mapped over both interpolation and + # non-interpolation dims : not supported. + cube = self.test_cube + cube.add_aux_coord(AuxCoord([[111, 112, 113, 114], + [211, 212, 213, 214]], + long_name='aux_0x'), + (0, 2)) + msg = 'Expected to find exactly one point.*Found 2' + with self.assertRaisesRegexp(Exception, msg): + interpolate(cube, self.single_sample_point, method='nearest') + + def test_metadata(self): + # Check exact result matching for a single point, with additional + # attributes and cell-methods. + cube = self.test_cube + cube.attributes['ODD_ATTR'] = 'string-value-example' + cube.add_cell_method(iris.coords.CellMethod('mean', 'area')) + result = interpolate(cube, self.single_sample_point, method='nearest') + # Check that the result is a single trajectory point, exactly equal to + # the expected part of the original data. + self.assertEqual(result.shape[-1], 1) + result = result[..., 0] + expected = cube[:, self.single_point_iy, self.single_point_ix] + self.assertEqual(result, expected) + + +if __name__ == "__main__": + tests.main()