diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py
index 1746834a65..ffe6d51c7f 100644
--- a/lib/iris/analysis/_grid_angles.py
+++ b/lib/iris/analysis/_grid_angles.py
@@ -92,7 +92,7 @@ def _angle(p, q, r):
if old_style:
mid_lons = np.deg2rad(q[0])
- pr = _3d_xyz_from_latlon(r[0], r[1]) - _3d_xyz_from_latlon(p[0], p[1])
+ pr = _3d_xyz_from_latlon(r[0], r[1]) - _3d_xyz_from_latlon(p[0], p[1])
pr_norm = np.sqrt(np.sum(pr**2, axis=0))
pr_top = pr[1] * np.cos(mid_lons) - pr[0] * np.sin(mid_lons)
@@ -124,7 +124,7 @@ def _angle(p, q, r):
lmb_hatvec_y,
lmb_hatvec_z)])
- pr = _3d_xyz_from_latlon(r[0], r[1]) - _3d_xyz_from_latlon(p[0], p[1])
+ pr = _3d_xyz_from_latlon(r[0], r[1]) - _3d_xyz_from_latlon(p[0], p[1])
# Dot products to form true-northward / true-eastward projections.
pr_cmpt_e = np.sum(pr * lmb_hatvec, axis=0)
@@ -141,7 +141,7 @@ def _angle(p, q, r):
rtol = 1.e-3
check = np.allclose(mag_rot, mag_orig, rtol=rtol)
if not check:
- print (mag_rot, mag_orig)
+ print(mag_rot, mag_orig)
assert np.allclose(mag_rot, mag_orig, rtol=rtol)
return psi
diff --git a/lib/iris/coords.py b/lib/iris/coords.py
index 8f1e48fd18..ff2f170a2f 100644
--- a/lib/iris/coords.py
+++ b/lib/iris/coords.py
@@ -1071,7 +1071,7 @@ def is_contiguous(self, rtol=1e-05, atol=1e-08):
rtol=rtol, atol=atol)
elif self.ndim == 2:
contiguous, _, _ = _discontiguity_in_2d_bounds(self.bounds,
- abs_tol=atol)
+ abs_tol=atol)
else:
contiguous = False
return contiguous
diff --git a/lib/iris/plot.py b/lib/iris/plot.py
index 8f11c8d116..7bf8ea07f3 100644
--- a/lib/iris/plot.py
+++ b/lib/iris/plot.py
@@ -334,7 +334,7 @@ def _draw_2d_from_bounds(draw_method_name, cube, *args, **kwargs):
plot_defn = _get_plot_defn(cube, mode, ndims=2)
twodim_contig_atol = kwargs.pop('two_dim_coord_contiguity_atol',
- 1e-4)
+ 1e-4)
for coord in plot_defn.coords:
if hasattr(coord, 'has_bounds'):
if coord.ndim == 2 and coord.has_bounds():
diff --git a/lib/iris/tests/stock.py b/lib/iris/tests/stock/__init__.py
similarity index 99%
rename from lib/iris/tests/stock.py
rename to lib/iris/tests/stock/__init__.py
index 1b55510d2b..5965d5a208 100644
--- a/lib/iris/tests/stock.py
+++ b/lib/iris/tests/stock/__init__.py
@@ -36,7 +36,8 @@
from iris.coords import DimCoord, AuxCoord
import iris.tests as tests
from iris.coord_systems import GeogCS, RotatedGeogCS
-
+from ._stock_2d_latlons import (sample_2d_latlons,
+ make_bounds_discontiguous_at_point)
def lat_lon_cube():
"""
diff --git a/lib/iris/tests/stock/_stock_2d_latlons.py b/lib/iris/tests/stock/_stock_2d_latlons.py
new file mode 100644
index 0000000000..75b62f31ba
--- /dev/null
+++ b/lib/iris/tests/stock/_stock_2d_latlons.py
@@ -0,0 +1,346 @@
+# (C) British Crown Copyright 2018, 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 .
+"""
+Extra stock routines for making and manipulating cubes with 2d coordinates,
+to mimic ocean grid data.
+
+"""
+from __future__ import (absolute_import, division, print_function)
+from six.moves import (filter, input, map, range, zip) # noqa
+
+import numpy as np
+import numpy.ma as ma
+
+from iris.analysis.cartography import unrotate_pole
+from iris.cube import Cube
+from iris.coords import AuxCoord, DimCoord
+from iris.coord_systems import RotatedGeogCS
+
+
+def expand_1d_x_and_y_bounds_to_2d_xy(x_bounds_1d, y_bounds_1d):
+ """
+ Convert bounds for separate 1-D X and Y coords into bounds for the
+ equivalent 2D coordinates.
+
+ The output arrays have 4 points per cell, for 4 'corners' of a gridcell,
+ in the usual anticlockwise order
+ (bottom-left, bottom-right, top-right, top-left).
+
+ If 1-dimensional X and Y coords have shapes nx and ny, then their bounds
+ have shapes (nx, 2) and (ny, 2).
+ The equivalent 2d coordinates would have values which are a "meshgrid" of
+ the original 1-D points, both having the shape (ny, ny).
+ The outputs are 2d bounds arrays suitable for such 2d coordinates.
+
+ Args:
+
+ * x_bounds_1d, y_bounds_1d : (array)
+ Coordinate bounds arrays, with shapes (nx, 2) and (ny, 2).
+
+ Result:
+
+ * x_bounds_2d, y_bounds_2d : (array)
+ Expanded 2d bounds arrays, both of shape (ny, nx, 4).
+
+ """
+ shapes = [bds.shape for bds in (x_bounds_1d, y_bounds_1d)]
+ for shape in shapes:
+ if len(shape) != 2 or shape[1] != 2:
+ msg = ('One-dimensional bounds arrays must have shapes (ny, 2) '
+ 'and (nx, 2). Got {} and {}.')
+ raise ValueError(msg.format(*shapes))
+
+ # Construct output arrays, which are both (ny, nx, 4).
+ nx, ny = [shape[0] for shape in shapes]
+ bds_2d_x = np.zeros((ny, nx, 4))
+ bds_2d_y = bds_2d_x.copy()
+
+ # Expand x bounds to 2D array (ny, nx, 4) : the same over all 'Y'.
+ # Bottom left+right corners are the original 1-D x bounds.
+ bds_2d_x[:, :, 0] = x_bounds_1d[:, 0].reshape((1, nx))
+ bds_2d_x[:, :, 1] = x_bounds_1d[:, 1].reshape((1, nx))
+ # Top left+right corners are the same as bottom left+right.
+ bds_2d_x[:, :, 2] = bds_2d_x[:, :, 1].copy()
+ bds_2d_x[:, :, 3] = bds_2d_x[:, :, 0].copy()
+
+ # Expand y bounds to 2D array (ny, nx, 4) : the same over all 'X'.
+ # Left-hand lower+upper corners are the original 1-D y bounds.
+ bds_2d_y[:, :, 0] = y_bounds_1d[:, 0].reshape((ny, 1))
+ bds_2d_y[:, :, 3] = y_bounds_1d[:, 1].reshape((ny, 1))
+ # Right-hand lower+upper corners are the same as the left-hand ones.
+ bds_2d_y[:, :, 1] = bds_2d_y[:, :, 0].copy()
+ bds_2d_y[:, :, 2] = bds_2d_y[:, :, 3].copy()
+
+ return bds_2d_x, bds_2d_y
+
+
+def grid_coords_2d_from_1d(x_coord_1d, y_coord_1d):
+ """
+ Calculate a pair of 2d X+Y coordinates from 1d ones, in a "meshgrid" style.
+ If the inputs are bounded, the outputs have 4 points per bounds in the
+ usual way, i.e. points 0,1,2,3 are the gridcell corners anticlockwise from
+ the bottom left.
+
+ """
+ for coord in (x_coord_1d, y_coord_1d):
+ if coord.ndim != 1:
+ msg = ('Input coords must be one-dimensional. '
+ 'Coordinate "{}" has shape {}.')
+ raise ValueError(msg.format(coord.name(), coord.shape))
+
+ # Calculate centre-points as a mesh of the 2 inputs.
+ pts_2d_x, pts_2d_y = np.meshgrid(x_coord_1d.points, y_coord_1d.points)
+ if not x_coord_1d.has_bounds() or not y_coord_1d.has_bounds():
+ bds_2d_x = None
+ bds_2d_y = None
+ else:
+ bds_2d_x, bds_2d_y = expand_1d_x_and_y_bounds_to_2d_xy(
+ x_coord_1d.bounds, y_coord_1d.bounds)
+
+ # Make two new coords + return them.
+ result = []
+ for pts, bds, name in zip((pts_2d_x, pts_2d_y),
+ (bds_2d_x, bds_2d_y),
+ ('longitude', 'latitude')):
+ coord = AuxCoord(pts, bounds=bds, standard_name=name, units='degrees')
+ result.append(coord)
+
+ return result
+
+
+def sample_2d_latlons(regional=False, rotated=False, transformed=False):
+ """
+ Construct small 2d cubes with 2d X and Y coordinates.
+
+ This makes cubes with 'expanded' coordinates (4 bounds per cell), analagous
+ to ORCA data.
+ The coordinates are always geographical, so either it has a coord system
+ or they are "true" lats + lons.
+ ( At present, they are always latitudes and longitudes, but maybe in a
+ rotated system. )
+ The results always have fully contiguous bounds.
+
+ Kwargs:
+ * regional (bool):
+ If False (default), results cover the whole globe, and there is
+ implicit connectivity between rhs + lhs of the array.
+ If True, coverage is regional and edges do not connect.
+ * rotated (bool):
+ If False, X and Y coordinates are true-latitudes and longitudes, with
+ an implicit coordinate system (i.e. None).
+ If True, the X and Y coordinates are lats+lons in a selected
+ rotated-latlon coordinate system.
+ * transformed (bool):
+ Build coords from rotated coords as for 'rotated', but then replace
+ their values with the equivalent "true" lats + lons, and no
+ coord-system (defaults to true-latlon).
+ In this case, the X and Y coords are no longer 'meshgrid' style,
+ i.e. the points + bounds values vary in *both* dimensions.
+
+ .. note::
+
+ 'transformed' is an alternative to 'rotated' : when 'transformed' is
+ set, then 'rotated' has no effect.
+
+ .. Some sample results printouts ::
+
+ >>> print(sample_2d_latlons())
+ test_data / (unknown) (-- : 5; -- : 6)
+ Auxiliary coordinates:
+ latitude x x
+ longitude x x
+ >>>
+ >>> print(sample_2d_latlons().coord(axis='x')[0, :2])
+ AuxCoord(array([ 37.5 , 93.75]),
+ bounds=array([[ 0. , 65.625, 65.625, 0. ],
+ [ 65.625, 121.875, 121.875, 65.625]]),
+ standard_name='longitude', units=Unit('degrees'))
+ >>> print(np.round(sample_2d_latlons().coord(axis='x').points, 3))
+ [[ 37.5 93.75 150. 206.25 262.5 318.75]
+ [ 37.5 93.75 150. 206.25 262.5 318.75]
+ [ 37.5 93.75 150. 206.25 262.5 318.75]
+ [ 37.5 93.75 150. 206.25 262.5 318.75]
+ [ 37.5 93.75 150. 206.25 262.5 318.75]]
+ >>> print(np.round(sample_2d_latlons().coord(axis='y').points, 3))
+ [[-85. -85. -85. -85. -85. -85. ]
+ [-47.5 -47.5 -47.5 -47.5 -47.5 -47.5]
+ [-10. -10. -10. -10. -10. -10. ]
+ [ 27.5 27.5 27.5 27.5 27.5 27.5]
+ [ 65. 65. 65. 65. 65. 65. ]]
+
+
+ >>> print(np.round(
+ sample_2d_latlons(rotated=True).coord(axis='x').points, 3))
+ [[ 37.5 93.75 150. 206.25 262.5 318.75]
+ [ 37.5 93.75 150. 206.25 262.5 318.75]
+ [ 37.5 93.75 150. 206.25 262.5 318.75]
+ [ 37.5 93.75 150. 206.25 262.5 318.75]
+ [ 37.5 93.75 150. 206.25 262.5 318.75]]
+ >>> print(sample_2d_latlons(rotated=True).coord(axis='y').coord_system)
+ RotatedGeogCS(75.0, 120.0)
+
+
+ >>> print(
+ sample_2d_latlons(transformed=True).coord(axis='y').coord_system)
+ None
+ >>> print(np.round(
+ sample_2d_latlons(transformed=True).coord(axis='x').points, 3))
+ [[ -50.718 -40.983 -46.74 -71.938 -79.293 -70.146]
+ [ -29.867 17.606 77.936 157.145 -141.037 -93.172]
+ [ -23.139 31.007 87.699 148.322 -154.639 -100.505]
+ [ -16.054 41.218 92.761 143.837 -164.738 -108.105]
+ [ 10.86 61.78 100.236 137.285 175.511 -135.446]]
+ >>> print(np.round(
+ sample_2d_latlons(transformed=True).coord(axis='y').points, 3))
+ [[-70.796 -74.52 -79.048 -79.26 -74.839 -70.96 ]
+ [-34.99 -46.352 -59.721 -60.34 -47.305 -35.499]
+ [ 1.976 -10.626 -22.859 -23.349 -11.595 1.37 ]
+ [ 38.914 25.531 14.312 13.893 24.585 38.215]
+ [ 74.197 60.258 51.325 51.016 59.446 73.268]]
+ >>>
+
+ """
+ def sample_cube(xargs, yargs):
+ # Make a test cube with given latitude + longitude coordinates.
+ # xargs/yargs are args for np.linspace (start, stop, N), to make the X
+ # and Y coordinate points.
+ x0, x1, nx = xargs
+ y0, y1, ny = yargs
+ # Data has cycling values, staggered a bit in successive rows.
+ data = np.zeros((ny, nx))
+ data.flat[:] = np.arange(ny * nx) % (nx + 2)
+ # Build a 2d cube with longitude + latitude coordinates.
+ cube = Cube(data, long_name='test_data')
+ x_pts = np.linspace(x0, x1, nx, endpoint=True)
+ y_pts = np.linspace(y0, y1, ny, endpoint=True)
+ co_x = DimCoord(x_pts, standard_name='longitude', units='degrees')
+ co_y = DimCoord(y_pts, standard_name='latitude', units='degrees')
+ cube.add_dim_coord(co_y, 0)
+ cube.add_dim_coord(co_x, 1)
+ return cube
+
+ if regional:
+ # Extract small region.
+ cube = sample_cube(xargs=(150., 243.75, 6), yargs=(-10., 40., 5))
+
+ # Make contiguous bounds.
+ for ax in ('x', 'y'):
+ cube.coord(axis=ax).guess_bounds()
+ else:
+ # Global data, but drastically reduced resolution.
+ cube = sample_cube(xargs=(37.5, 318.75, 6), yargs=(-85., 65., 5))
+
+ # Patch bounds to ensure it is still contiguous + global.
+ for name in ('longitude', 'latitude'):
+ coord = cube.coord(name)
+ coord.guess_bounds()
+ bds = coord.bounds.copy()
+ # Make bounds global, by fixing lowest and uppermost values.
+ if name == 'longitude':
+ bds[0, 0] = 0.0
+ bds[-1, 1] = 360.0
+ else:
+ bds[0, 0] = -90.0
+ bds[-1, 1] = 90.0
+ coord.bounds = bds
+
+ # Get 1d coordinate points + bounds + calculate 2d equivalents.
+ co_1d_x, co_1d_y = [cube.coord(axis=ax).copy() for ax in ('x', 'y')]
+ co_2d_x, co_2d_y = grid_coords_2d_from_1d(co_1d_x, co_1d_y)
+
+ # Remove the old grid coords.
+ for coord in (co_1d_x, co_1d_y):
+ cube.remove_coord(coord)
+
+ # Add the new grid coords.
+ for coord in (co_2d_x, co_2d_y):
+ cube.add_aux_coord(coord, (0, 1))
+
+ if transformed or rotated:
+ # Take the lats + lons as being in a rotated coord system.
+ pole_lat, pole_lon = 75.0, 120.0
+ if transformed:
+ # Reproject coordinate values from rotated to true lat-lons.
+ co_x, co_y = [cube.coord(axis=ax) for ax in ('x', 'y')]
+ # Unrotate points.
+ lons, lats = co_x.points, co_y.points
+ lons, lats = unrotate_pole(lons, lats, pole_lon, pole_lat)
+ co_x.points, co_y.points = lons, lats
+ # Unrotate bounds.
+ lons, lats = co_x.bounds, co_y.bounds
+ # Note: save the shape, flatten + then re-apply the shape, because
+ # "unrotate_pole" uses "cartopy.crs.CRS.transform_points", which
+ # only works on arrays of 1 or 2 dimensions.
+ shape = lons.shape
+ lons, lats = unrotate_pole(lons.flatten(), lats.flatten(),
+ pole_lon, pole_lat)
+ co_x.bounds, co_y.bounds = lons.reshape(shape), lats.reshape(shape)
+ else:
+ # "Just" rotate operation : add a coord-system to each coord.
+ cs = RotatedGeogCS(pole_lat, pole_lon)
+ for coord in cube.coords():
+ coord.coord_system = cs
+
+ return cube
+
+
+def make_bounds_discontiguous_at_point(cube, at_iy, at_ix):
+ """
+ Meddle with the XY grid bounds of a cube to make the grid discontiguous.
+
+ Changes the points and bounds of a single gridcell, so that it becomes
+ discontinuous with the next gridcell to its right.
+ Also masks the cube data at the given point.
+
+ The cube must have bounded 2d 'x' and 'y' coordinates.
+
+ TODO: add a switch to make a discontiguity in the *y*-direction instead ?
+
+ """
+ x_coord = cube.coord(axis='x')
+ y_coord = cube.coord(axis='y')
+ assert x_coord.shape == y_coord.shape
+ assert (coord.bounds.ndim == 3 and coord.shape[-1] == 4
+ for coord in (x_coord, y_coord))
+
+ # For both X and Y coord, move points + bounds to create a discontinuity.
+ def adjust_coord(coord):
+ pts, bds = coord.points, coord.bounds
+ # Fetch the 4 bounds (bottom-left, bottom-right, top-right, top-left)
+ bds_bl, bds_br, bds_tr, bds_tl = bds[at_iy, at_ix]
+ # Make a discontinuity "at" (iy, ix), by moving the right-hand edge of
+ # the cell to the midpoint of the existing left+right bounds.
+ new_bds_br = 0.5 * (bds_bl + bds_br)
+ new_bds_tr = 0.5 * (bds_tl + bds_tr)
+ bds_br, bds_tr = new_bds_br, new_bds_tr
+ bds[at_iy, at_ix] = [bds_bl, bds_br, bds_tr, bds_tl]
+ # Also reset the cell midpoint to the middle of the 4 new corners,
+ # in case having a midpoint outside the corners might cause a problem.
+ new_pt = 0.25 * sum([bds_bl, bds_br, bds_tr, bds_tl])
+ pts[at_iy, at_ix] = new_pt
+ # Write back the coord points+bounds (can only assign whole arrays).
+ coord.points, coord.bounds = pts, bds
+
+ adjust_coord(x_coord)
+ adjust_coord(y_coord)
+ # Also mask the relevant data point.
+ data = cube.data # N.B. fetch all the data.
+ if not ma.isMaskedArray(data):
+ # Promote to masked array, to avoid converting mask to NaN.
+ data = ma.masked_array(data)
+ data[at_iy, at_ix] = ma.masked
+ cube.data = data
diff --git a/lib/iris/tests/test_coding_standards.py b/lib/iris/tests/test_coding_standards.py
index 39c73e8285..c24ac46cb3 100644
--- a/lib/iris/tests/test_coding_standards.py
+++ b/lib/iris/tests/test_coding_standards.py
@@ -93,7 +93,7 @@ class StandardReportWithExclusions(pep8.StandardReport):
'*/iris/io/format_picker.py',
'*/iris/tests/__init__.py',
'*/iris/tests/pp.py',
- '*/iris/tests/stock.py',
+ '*/iris/tests/stock/__init__.py',
'*/iris/tests/system_test.py',
'*/iris/tests/test_analysis.py',
'*/iris/tests/test_analysis_calculus.py',
diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
deleted file mode 100644
index de00d2bc76..0000000000
--- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
+++ /dev/null
@@ -1,392 +0,0 @@
-# (C) British Crown Copyright 2018, 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 function
-:func:`iris.analysis.cartography.gridcell_angles`.
-
-"""
-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
-import numpy.ma as ma
-
-import cartopy.crs as ccrs
-from iris.cube import Cube
-from iris.coords import DimCoord, AuxCoord
-import iris.coord_systems
-from iris.analysis.cartography import unrotate_pole
-
-from iris.analysis.cartography import (gridcell_angles,
- rotate_grid_vectors)
-
-import matplotlib.pyplot as plt
-from orca_utils.plot_testing.blockplot_from_bounds import blockplot_2dll
-
-
-def _rotated_grid_sample(pole_lat=15, pole_lon=-180,
- lon_bounds=np.linspace(-30, 30, 6, endpoint=True),
- lat_bounds=np.linspace(-30, 30, 6, endpoint=True)):
- # Calculate *true* lat_bounds+lon_bounds for the rotated grid.
- lon_bounds = np.array(lon_bounds, dtype=float)
- lat_bounds = np.array(lat_bounds, dtype=float)
- # Construct centrepoints.
- lons = 0.5 * (lon_bounds[:-1] + lon_bounds[1:])
- lats = 0.5 * (lat_bounds[:-1] + lat_bounds[1:])
- # Convert all to full 2d arrays.
- lon_bounds, lat_bounds = np.meshgrid(lon_bounds, lat_bounds)
- lons, lats = np.meshgrid(lons, lats)
- # Calculate true lats+lons for all points.
- lons_true_bds, lats_true_bds = unrotate_pole(lon_bounds, lat_bounds,
- pole_lon, pole_lat)
- lons_true, lats_true = unrotate_pole(lons, lats, pole_lon, pole_lat)
- # Make the 'unified' bounds into contiguous (ny, nx, 4) arrays.
- def expand_unified_bds(bds):
- ny, nx = bds.shape
- bds_4 = np.zeros((ny - 1, nx - 1, 4))
- bds_4[:, :, 0] = bds[:-1, :-1]
- bds_4[:, :, 1] = bds[:-1, 1:]
- bds_4[:, :, 2] = bds[1:, 1:]
- bds_4[:, :, 3] = bds[1:, :-1]
- return bds_4
-
- lon_true_bds4, lat_true_bds4 = (expand_unified_bds(bds)
- for bds in (lons_true_bds, lats_true_bds))
- # Make these into a 2d-latlon grid for a cube
- cube = Cube(np.zeros(lon_true_bds4.shape[:-1]))
- co_x = AuxCoord(lons_true, bounds=lon_true_bds4,
- standard_name='longitude', units='degrees')
- co_y = AuxCoord(lats_true, bounds=lat_true_bds4,
- standard_name='latitude', units='degrees')
- cube.add_aux_coord(co_x, (0, 1))
- cube.add_aux_coord(co_y, (0, 1))
- return cube
-
-
-class TestGridcellAngles(tests.IrisTest):
- def _singlecell_30deg_cube(self, x0=90., y0=0., dx=20., dy=10.):
- x_pts = np.array([[x0]])
- y_pts = np.array([[y0]])
- x_bds = x0 + dx * np.array([[[-1., 1, 0.5, -1.5]]])
-# self.assertArrayAllClose(x_bds, np.array([[[70., 110, 100, 60]]]))
- y_bds = y0 + dy * np.array([[[-1., 1, 3, 1]]])
-# self.assertArrayAllClose(y_bds, np.array([[[-10., 10, 30, 10]]]))
- co_x = AuxCoord(points=x_pts, bounds=x_bds,
- standard_name='longitude', units='degrees')
- co_y = AuxCoord(points=y_pts, bounds=y_bds,
- standard_name='latitude', units='degrees')
- cube = Cube(np.zeros((1, 1)))
- cube.add_aux_coord(co_x, (0, 1))
- cube.add_aux_coord(co_y, (0, 1))
- return cube
-
- def _singlecell_diamond_cube(self, x0=90., y0=0., dy=10., dx_eq=None):
- if dx_eq is None:
- dx_eq = dy
- x_pts = np.array([[x0]])
- y_pts = np.array([[y0]])
- dx = dx_eq / np.cos(np.deg2rad(y0))
- x_bds = np.array([[[x0, x0 + dx, x0, x0 - dx]]])
- y_bds = np.array([[[y0 - dy, y0, y0 + dy, y0]]])
- co_x = AuxCoord(points=x_pts, bounds=x_bds,
- standard_name='longitude', units='degrees')
- co_y = AuxCoord(points=y_pts, bounds=y_bds,
- standard_name='latitude', units='degrees')
- cube = Cube(np.zeros((1, 1)))
- cube.add_aux_coord(co_x, (0, 1))
- cube.add_aux_coord(co_y, (0, 1))
- return cube
-
- def test_single_cell_equatorial(self):
- plt.switch_backend('tkagg')
- plt.figure(figsize=(10,10))
- ax = plt.axes(projection=ccrs.Mercator())
-# ax = plt.axes(projection=ccrs.NorthPolarStereo())
-# ax = plt.axes(projection=ccrs.Orthographic(central_longitude=90.,
-# central_latitude=30.))
-
- lon0 = 90.0
- dy = 2.0
- y_0, y_n, ny = -80, 80, 9
- angles = []
- for lat in np.linspace(y_0, y_n, ny):
- cube = self._singlecell_diamond_cube(x0=lon0, y0=lat, dy=dy)
- angles_cube = gridcell_angles(cube,
-# cell_angle_boundpoints='mid-lhs, mid-rhs')
- cell_angle_boundpoints='lower-left, lower-right')
- tmp_cube = angles_cube.copy()
- tmp_cube.convert_units('degrees')
- print('')
- print(lat)
- co_x, co_y = (cube.coord(axis=ax) for ax in ('x', 'y'))
- print()
- print(' at : {}, {}'.format(co_x.points[0, 0], co_y.points[0, 0]))
- print(' x-bds:')
- print(co_x.bounds)
- print(' y-bds:')
- print(co_y.bounds)
- angle = tmp_cube.data[0, 0]
- angles.append(angle)
- print(angle)
- blockplot_2dll(cube)
-
- ax.coastlines()
- ax.set_global()
-
- # Plot constant NEly (45deg) arrows.
- xx = np.array([lon0] * ny)
- yy = np.linspace(y_0, y_n, ny) - dy
- uu = np.array([1.0] * ny)
- plt.quiver(xx, yy,
- uu, np.cos(np.deg2rad(yy)),
- zorder=2, color='red',
- scale_units='xy',
- transform=ccrs.PlateCarree())
-
- # Also plot returned angles.
- angles_arr_rad = np.deg2rad(angles)
- u_arr = uu * np.cos(angles_arr_rad)
- v_arr = uu * np.sin(angles_arr_rad) * np.cos(np.deg2rad(yy))
-
- plt.quiver(xx, yy,
- u_arr,
- v_arr,
- zorder=2, color='magenta',
- scale_units='xy',
- width=0.005,
- scale=0.2e-6,
-# width=0.5,
- transform=ccrs.PlateCarree())
-
- plt.show()
-
-
- def test_values(self):
- # Construct a rotated-pole grid and check angle calculation.
- testcube = _rotated_grid_sample()
-
- cell_angle_boundpoints = 'mid-lhs, mid-rhs'
-# cell_angle_boundpoints = 'lower-left, lower-right'
-# cell_angle_boundpoints = 'garble'
- angles_cube = gridcell_angles(
- testcube,
- cell_angle_boundpoints=cell_angle_boundpoints)
- angles_cube.convert_units('radians')
-
- # testing phase...
- print(np.rad2deg(angles_cube.data))
-
- import matplotlib.pyplot as plt
- plt.switch_backend('tkagg')
-
-# plot_map = 'north_polar_stereographic'
-# plot_map = 'plate_carree'
-# plot_map = 'mercator'
- plot_map = 'north_polar_orthographic'
- if plot_map == 'plate_carree':
- scale = 0.1
- map_proj = ccrs.PlateCarree()
- elif plot_map == 'mercator':
- scale = 3.0e-6
- map_proj = ccrs.Mercator()
- map_proj._threshold *= 0.01
- elif plot_map == 'north_polar_orthographic':
- scale = 3.0e-6
- map_proj = ccrs.Orthographic(central_longitude=0.0,
- central_latitude=90.0,)
- map_proj._threshold *= 0.01
- elif plot_map == 'north_polar_stereographic':
- scale = 3.0e-6
- map_proj = ccrs.NorthPolarStereo()
- else:
- assert 0
-
- ax = plt.axes(projection=map_proj)
- data_proj = ccrs.PlateCarree()
-
- deg_scale = 10.0
-
-# angles = 'uv'
- angles = 'xy'
-
- ax.coastlines()
- ax.gridlines()
- for i_bnd in range(4):
- color = ['black', 'red', 'blue', 'magenta'][i_bnd]
- plt.plot(testcube.coord('longitude').bounds[..., i_bnd],
- testcube.coord('latitude').bounds[..., i_bnd],
- '+', markersize=10., markeredgewidth=2.,
- markerfacecolor=color, markeredgecolor=color,
- transform=data_proj)
-
-
- # Show plain 0,1 + 1,0 (PlateCarree) vectors unrotated at the given points.
- pts_shape = testcube.coord('longitude').shape
- ny, nx = pts_shape
- u0 = np.ones(pts_shape)
- v0 = np.zeros(pts_shape)
- u1 = v0.copy()
- v1 = u0.copy()
-
- x0s = testcube.coord('longitude').points
- y0s = testcube.coord('latitude').points
- yscale = np.cos(np.deg2rad(y0s))
- plt.quiver(x0s, y0s, u0, v0 * yscale,
- color='blue', width=0.005,
- headwidth=2., # headlength=1.0, headaxislength=0.7,
- angles=angles,
- scale_units='xy', scale=scale,
- transform=data_proj)
- plt.quiver(x0s, y0s, u1, v1 * yscale,
- color='red', width=0.005,
- headwidth=2., # headlength=1.0, headaxislength=0.7,
- angles=angles,
- scale_units='xy', scale=scale,
- transform=data_proj)
-
- # Add 45deg arrows (NEly), still on a PlateCarree map.
- plt.quiver(x0s, y0s, v1, v1 * yscale,
- color='green', width=0.005,
- headwidth=2., # headlength=1.0, headaxislength=0.7,
- angles=angles,
- scale_units='xy', scale=scale,
- transform=data_proj)
-
-
-
- #
- # Repeat the above plotting short lines INSTEAD of quiver.
- #
- u0d = x0s + deg_scale * u0
- v0d = y0s + deg_scale * v0
- u1d = x0s + deg_scale * u1
- v1d = y0s + deg_scale * v1
- u2d = x0s + deg_scale * u0
- v2d = y0s + deg_scale * v1
- for iy in range(ny):
- for ix in range(nx):
- plt.plot([x0s[iy, ix], u0d[iy, ix]],
- [y0s[iy, ix], v0d[iy, ix]],
- ':', color='blue', linewidth=0.5,
- transform=data_proj)
- plt.plot([x0s[iy, ix], u1d[iy, ix]],
- [y0s[iy, ix], v1d[iy, ix]],
- ':', color='red', linewidth=0.5,
- transform=data_proj)
- plt.plot([x0s[iy, ix], u2d[iy, ix]],
- [y0s[iy, ix], v2d[iy, ix]],
- ':', color='green', linewidth=0.5,
- transform=data_proj)
-
-
- # Overplot BL-BR and BL-TL lines from the cell bounds.
- co_lon, co_lat = [testcube.coord(name).copy()
- for name in ('longitude', 'latitude')]
- for co in (co_lon, co_lat):
- co.convert_units('degrees')
- lon_bds, lat_bds = [co.bounds for co in (co_lon, co_lat)]
-# ny, nx = lon_bds.shape[:-1]
- for iy in range(ny):
- for ix in range(nx):
- x0, y0 = lon_bds[iy, ix, 0], lat_bds[iy, ix, 0]
- x1, y1 = lon_bds[iy, ix, 1], lat_bds[iy, ix, 1]
- x2, y2 = lon_bds[iy, ix, 3], lat_bds[iy, ix, 3]
- plt.plot([x0, x1], [y0, y1], 'x-',
- color='orange',
- transform=data_proj)
- plt.plot([x0, x2], [y0, y2], 'x-',
- color='orange', linestyle='--',
- transform=data_proj)
-
- # Plot U0, rotated by cell angles, also at cell bottom-lefts.
- u0_cube, u1_cube, v0_cube, v1_cube = [testcube.copy(data=aa)
- for aa in (u0, v0, u1, v1)]
- u0r_cube, v0r_cube = rotate_grid_vectors(
- u0_cube, v0_cube, grid_angles_cube=angles_cube)
- u0r, v0r = [cube.data for cube in (u0r_cube, v0r_cube)]
-
- xbl, ybl = lon_bds[..., 0], lat_bds[..., 0]
- #
- # Replace quiver here with delta-based lineplot
- #
- urd = xbl + deg_scale * u0r
- vrd = ybl + deg_scale * v0r * yscale
- for iy in range(ny):
- for ix in range(nx):
- plt.plot([xbl[iy, ix], urd[iy, ix]],
- [ybl[iy, ix], vrd[iy, ix]],
- ':', color='magenta', linewidth=2.5,
- transform=data_proj)
- # Show this is the SAME as lineplot
- plt.quiver(xbl, ybl, u0r, v0r * yscale,
- color='magenta', width=0.01,
- headwidth=1.2, # headlength=1.0, headaxislength=0.7,
- angles=angles,
- scale_units='xy', scale=scale,
- transform=data_proj)
-
- plt.suptitle('angles from "{}"'.format(cell_angle_boundpoints))
-
-# # Also draw small lines pointing at the correct (TRUE, not ) angle.
-# ny, nx = x0s.shape
-# size_degrees = 1.0
-# angles = angles_cube.copy()
-# angles.convert_units('radians')
-# angles = angles.data
-# lats = testcube.coord('latitude').copy()
-# lats.convert_units('radians')
-# lats = lats.points
-# dxs = size_degrees * u0.copy() #* np.cos(angles)
-# dys = size_degrees * u0.copy() # / np.sqrt(np.cos(lats))
-# x1s = x0s + dxs
-# y1s = y0s + dys
-## for iy in range(ny):
-## for ix in range(nx):
-## plt.plot([x0s[iy, ix], x1s[iy, ix]],
-## [y0s[iy, ix], y1s[iy, ix]],
-## 'o-', markersize=4., markeredgewidth=0.,
-## color='green', # scale_units='xy', scale=scale,
-## transform=data_proj)
-# plt.quiver(x0s, y0s, dxs, dys,
-# color='green', linewidth=0.2,
-# angles=angles,
-# scale_units='xy', scale=scale * 0.6,
-# transform=data_proj)
-
-
-
- ax.set_global()
- plt.show()
-
- angles_cube.convert_units('degrees')
-
- self.assertArrayAllClose(
- angles_cube.data,
- [[33.421, 17.928, 0., -17.928, -33.421],
- [41.981, 24.069, 0., -24.069, -41.981],
- [56.624, 37.809, 0., -37.809, -56.624],
- [79.940, 74.227, 0., -74.227, -79.940],
- [107.313, 126.361, -180., -126.361, -107.313]],
- atol=0.002)
-
-
-if __name__ == "__main__":
- tests.main()
diff --git a/lib/iris/tests/unit/plot/test_2d_coords.py b/lib/iris/tests/unit/plot/test_2d_coords.py
index bb27341baa..7397921861 100644
--- a/lib/iris/tests/unit/plot/test_2d_coords.py
+++ b/lib/iris/tests/unit/plot/test_2d_coords.py
@@ -22,20 +22,25 @@
# Import iris.tests first so that some things can be initialised before
# importing anything else.
import iris.tests as tests
-import iris.coords as coords
import numpy.ma as ma
+import iris.coords as coords
from iris.tests.stock import simple_2d_w_multidim_coords as cube_2dcoords
from iris.tests.stock import simple_3d_w_multidim_coords as cube3d_2dcoords
-from iris.tests.stock import lat_lon_cube
+from iris.tests.stock import sample_2d_latlons
+from iris.tests.stock import make_bounds_discontiguous_at_point
-from orca_utils.plot_testing import testdata_2d_coords as testdata
if tests.MPL_AVAILABLE:
import iris.plot as iplt
+def full2d_global():
+ return sample_2d_latlons(transformed=True)
+
+
+@tests.skip_data
class Test_2d_coords_plot_defn_bound_mode(tests.IrisTest):
def setUp(self):
self.multidim_cube = cube_2dcoords()
@@ -43,8 +48,8 @@ def setUp(self):
# latlon_2d is a cube with 2d coords, 4 bounds per point,
# discontiguities in the bounds but masked data at the discontiguities.
- self.latlon_2d = testdata.full2d_global()
- testdata.make_bounds_discontiguous_at_point(self.latlon_2d, 2, 2)
+ self.latlon_2d = full2d_global()
+ make_bounds_discontiguous_at_point(self.latlon_2d, 2, 2)
# # Take a latlon cube with 1D coords, broadcast the coords into 2D
# # ones, then add ONE of them back into the cube in place of original:
@@ -106,9 +111,10 @@ def test_total_span_check(self):
# # TODO Find out where I can put a catch for this (if necessary)
# cube = self.mixed_dims
# with self.assertRaises(ValueError):
- # iplt._get_plot_defn_custom_coords_picked(cube,
- # ('latitude', 'longitude'),
- # self.mode)
+ # iplt._get_plot_defn_custom_coords_picked(
+ # cube,
+ # ('latitude', 'longitude'),
+ # self.mode)
def test_map_common_not_enough_bounds(self):
# Test that a lat-lon cube with 2d coords and 2 bounds per point
@@ -129,16 +135,16 @@ def test_map_common_2d(self):
result = iplt._map_common('pcolor', None, self.mode, cube, plot_defn)
self.assertTrue(result)
- # def test_discontiguous_masked(self):
- # # Test that a contiguity check will raise a warning (not an error) for
- # # discontiguous bounds but appropriately masked data.
- # cube = self.latlon_2d
- # coord = cube.coord('longitude')
- # msg = 'The bounds of the longitude coordinate are not contiguous. ' \
- # 'However, data is masked where the discontiguity occurs so ' \
- # 'plotting anyway.'
- # with self.assertWarnsRegexp(msg):
- # iplt._check_contiguity_and_bounds(coord, cube.data)
+# def test_discontiguous_masked(self):
+# # Test that a contiguity check will raise a warning (not an error) for
+# # discontiguous bounds but appropriately masked data.
+# cube = self.latlon_2d
+# coord = cube.coord('longitude')
+# msg = 'The bounds of the longitude coordinate are not contiguous. ' \
+# 'However, data is masked where the discontiguity occurs so ' \
+# 'plotting anyway.'
+# with self.assertWarnsRegexp(msg):
+# iplt._check_contiguity_and_bounds(coord, cube.data)
def test_discontiguous_unmasked(self):
# Check that an error occurs when the contiguity check finds
@@ -155,3 +161,7 @@ def test_draw_2d_from_bounds(self):
cube = self.latlon_2d
result = iplt._draw_2d_from_bounds('pcolormesh', cube)
self.assertTrue(result)
+
+
+if __name__ == '__main__':
+ tests.main()