- Iris v2.1
+ Iris v2.2
A powerful, format-agnostic, community-driven Python library for analysing and
diff --git a/docs/iris/src/userguide/interpolation_and_regridding.rst b/docs/iris/src/userguide/interpolation_and_regridding.rst
index 4f2b06de44..e3cd622541 100644
--- a/docs/iris/src/userguide/interpolation_and_regridding.rst
+++ b/docs/iris/src/userguide/interpolation_and_regridding.rst
@@ -5,6 +5,8 @@
import numpy as np
import iris
+ import warnings
+ warnings.simplefilter('ignore')
=================================
Cube interpolation and regridding
@@ -137,10 +139,9 @@ This cube has a "hybrid-height" vertical coordinate system, meaning that the ver
coordinate is unevenly spaced in altitude:
>>> print(column.coord('altitude').points)
- [418.6983642578125 434.57049560546875 456.79278564453125 485.3664855957031
- 520.2932739257812 561.5751953125 609.2144775390625 663.214111328125
- 723.5769653320312 790.306640625 863.4072265625 942.88232421875
- 1028.737060546875 1120.9764404296875 1219.6051025390625]
+ [ 418.69836 434.5705 456.7928 485.3665 520.2933 561.5752
+ 609.2145 663.2141 723.57697 790.30664 863.4072 942.8823
+ 1028.737 1120.9764 1219.6051 ]
We could regularise the vertical coordinate by defining 10 equally spaced altitude
sample points between 400 and 1250 and interpolating our vertical coordinate onto
@@ -184,9 +185,8 @@ For example, to mask values that lie beyond the range of the original data:
>>> scheme = iris.analysis.Linear(extrapolation_mode='mask')
>>> new_column = column.interpolate(sample_points, scheme)
>>> print(new_column.coord('altitude').points)
- [nan 494.44451904296875 588.888916015625 683.333251953125 777.77783203125
- 872.2222290039062 966.666748046875 1061.111083984375 1155.555419921875
- nan]
+ [ nan 494.44452 588.8889 683.33325 777.77783 872.2222
+ 966.66675 1061.1111 1155.5554 nan]
.. _caching_an_interpolator:
diff --git a/docs/iris/src/whatsnew/2.2.rst b/docs/iris/src/whatsnew/2.2.rst
new file mode 100644
index 0000000000..4edc305025
--- /dev/null
+++ b/docs/iris/src/whatsnew/2.2.rst
@@ -0,0 +1,36 @@
+What's New in Iris 2.2.0
+************************
+
+:Release: 2.2.0a0
+:Date:
+
+
+This document explains the new/changed features of Iris in the alpha release
+of version 2.2.0
+(:doc:`View all changes `).
+
+
+Iris 2.2.0 Features
+===================
+.. _showcase:
+
+.. admonition:: 2-Dimensional Coordinate Plotting
+
+ The iris plot functions :func:`~iris.plot.pcolor` and
+ :func:`~iris.plot.pcolormesh` now accommodate the plotting of 2-dimensional
+ coordinates as well as 1-dimensional coordinates.
+
+ To enable this feature, each coordinate passed in for plotting will be
+ automatically checked for contiguity. Coordinate bounds must either be
+ contiguous, or the cube's data must be masked at the discontiguities in
+ order to avoid plotting errors.
+
+ The iris plot functions :func:`iris.plot.quiver` and
+ :func:`iris.plot.streamplot` have been added, and these also work with
+ 2-dimensional plot coordinates.
+
+.. admonition:: 2-Dimensional Grid Vectors
+
+ The iris functions :func:`iris.analysis.cartography.gridcell_angles` and
+ :func:`iris.analysis.cartography.rotate_grid_vectors` have been added,
+ allowing you to convert gridcell-oriented vectors to true-North/East ones.
diff --git a/docs/iris/src/whatsnew/index.rst b/docs/iris/src/whatsnew/index.rst
index 104c5074ca..c3a34303f0 100644
--- a/docs/iris/src/whatsnew/index.rst
+++ b/docs/iris/src/whatsnew/index.rst
@@ -9,6 +9,7 @@ Iris versions.
.. toctree::
:maxdepth: 2
+ 2.2.rst
2.1.rst
2.0.rst
1.13.rst
diff --git a/lib/iris/__init__.py b/lib/iris/__init__.py
index 2e7eba5892..882b52aaea 100644
--- a/lib/iris/__init__.py
+++ b/lib/iris/__init__.py
@@ -121,7 +121,7 @@ def callback(cube, field, filename):
# Iris revision.
-__version__ = '2.2.0dev0'
+__version__ = '2.2.0a0'
# Restrict the names imported when using "from iris import *"
__all__ = ['load', 'load_cube', 'load_cubes', 'load_raw',
diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py
new file mode 100644
index 0000000000..90876840af
--- /dev/null
+++ b/lib/iris/analysis/_grid_angles.py
@@ -0,0 +1,451 @@
+# (C) British Crown Copyright 2010 - 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 .
+"""
+Code to implement vector rotation by angles, and inferring gridcell angles
+from coordinate points and bounds.
+
+"""
+from __future__ import (absolute_import, division, print_function)
+from six.moves import (filter, input, map, range, zip) # noqa
+
+import numpy as np
+
+import cartopy.crs as ccrs
+import iris
+
+
+def _3d_xyz_from_latlon(lon, lat):
+ """
+ Return locations of (lon, lat) in 3D space.
+
+ Args:
+
+ * lon, lat: (float array)
+ Arrays of longitudes and latitudes, in degrees.
+ Both the same shape.
+
+ Returns:
+
+ * xyz : (array, dtype=float64)
+ Cartesian coordinates on a unit sphere.
+ Shape is (3, ).
+ The x / y / z coordinates are in xyz[0 / 1 / 2].
+
+ """
+ lon1 = np.deg2rad(lon).astype(np.float64)
+ lat1 = np.deg2rad(lat).astype(np.float64)
+
+ x = np.cos(lat1) * np.cos(lon1)
+ y = np.cos(lat1) * np.sin(lon1)
+ z = np.sin(lat1)
+
+ result = np.concatenate([array[np.newaxis] for array in (x, y, z)])
+
+ return result
+
+
+def _latlon_from_xyz(xyz):
+ """
+ Return arrays of lons+lats angles from xyz locations.
+
+ Args:
+
+ * xyz: (array)
+ Array of 3-D cartesian coordinates.
+ Shape (3, ).
+ x / y / z values are in xyz[0 / 1 / 2],
+
+ Returns:
+
+ * lonlat : (array)
+ longitude and latitude position angles, in degrees.
+ Shape (2, ).
+ The longitudes / latitudes are in lonlat[0 / 1].
+
+ """
+ lons = np.rad2deg(np.arctan2(xyz[1], xyz[0]))
+ radii = np.sqrt(np.sum(xyz * xyz, axis=0))
+ lats = np.rad2deg(np.arcsin(xyz[2] / radii))
+ return np.array([lons, lats])
+
+
+def _angle(p, q, r):
+ """
+ Estimate grid-angles to true-Eastward direction from positions in the same
+ grid row, but at increasing column (grid-Eastward) positions.
+
+ {P, Q, R} are locations of consecutive points in the same grid row.
+ These could be successive points in a single grid,
+ e.g. {T(i-1,j), T(i,j), T(i+1,j)}
+ or a mixture of U/V and T gridpoints if row positions are aligned,
+ e.g. {v(i,j), f(i,j), v(i+1,j)}.
+
+ Method:
+
+ Calculate dot product of a unit-vector parallel to P-->R, unit-scaled,
+ with the unit eastward (true longitude) vector at Q.
+ This value is cos(required angle).
+ Discriminate between +/- angles by comparing latitudes of P and R.
+ Return NaN where any P-->R are zero.
+
+ .. NOTE::
+
+ This method assumes that the vector PR is parallel to the surface
+ at the longitude of Q, as it uses the length of PR as the basis for
+ the cosine ratio.
+ That is only exact when Q is at the same longitude as the midpoint
+ of PR, and this typically causes errors which grow with increasing
+ gridcell angle.
+ However, we retain this method because it reproduces the "standard"
+ gridcell-orientation-angle arrays found in files output by the CICE
+ model, which presumably uses an equivalent calculation.
+
+ Args:
+
+ * p, q, r : (float array)
+ Arrays of angles, in degrees.
+ All the same shape.
+ Shape is (2, ).
+ Longitudes / latitudes are in array[0 / 1].
+
+ Returns:
+
+ * angle : (float array)
+ Grid angles relative to true-East, in degrees.
+ Positive when grid-East is anticlockwise from true-East.
+ Shape is same as .
+
+ """
+ mid_lons = np.deg2rad(q[0])
+
+ 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)
+
+ index = pr_norm == 0
+ pr_norm[index] = 1
+
+ cosine = np.maximum(np.minimum(pr_top / pr_norm, 1), -1)
+ cosine[index] = 0
+
+ psi = np.arccos(cosine) * np.sign(r[1] - p[1])
+ psi[index] = np.nan
+
+ return np.rad2deg(psi)
+
+
+def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'):
+ """
+ Calculate gridcell orientations for an arbitrary 2-dimensional grid.
+
+ The input grid is defined by two 2-dimensional coordinate arrays with the
+ same dimensions (ny, nx), specifying the geolocations of a 2D mesh.
+
+ Input values may be coordinate points (ny, nx) or bounds (ny, nx, 4).
+ However, if points, the edges in the X direction are assumed to be
+ connected by wraparound.
+
+ Input can be either two arrays, two coordinates, or a single cube
+ containing two suitable coordinates identified with the 'x' and'y' axes.
+
+ Args:
+
+ The inputs (x [,y]) can be any of the folliwing :
+
+ * x (:class:`~iris.cube.Cube`):
+ a grid cube with 2D X and Y coordinates, identified by 'axis'.
+ The coordinates must be 2-dimensional with the same shape.
+ The two dimensions represent grid dimensions in the order Y, then X.
+
+ * x, y (:class:`~iris.coords.Coord`):
+ X and Y coordinates, specifying grid locations on the globe.
+ The coordinates must be 2-dimensional with the same shape.
+ The two dimensions represent grid dimensions in the order Y, then X.
+ If there is no coordinate system, they are assumed to be true
+ longitudes and latitudes. Units must convertible to 'degrees'.
+
+ * x, y (2-dimensional arrays of same shape (ny, nx)):
+ longitude and latitude cell center locations, in degrees.
+ The two dimensions represent grid dimensions in the order Y, then X.
+
+ * x, y (3-dimensional arrays of same shape (ny, nx, 4)):
+ longitude and latitude cell bounds, in degrees.
+ The first two dimensions are grid dimensions in the order Y, then X.
+ The last index maps cell corners anticlockwise from bottom-left.
+
+ Optional Args:
+
+ * cell_angle_boundpoints (string):
+ Controls which gridcell bounds locations are used to calculate angles,
+ if the inputs are bounds or bounded coordinates.
+ Valid values are 'lower-left, lower-right', which takes the angle from
+ the lower left to the lower right corner, and 'mid-lhs, mid-rhs' which
+ takes an angles between the average of the left-hand and right-hand
+ pairs of corners. The default is 'mid-lhs, mid-rhs'.
+
+ Returns:
+
+ angles : (2-dimensional cube)
+
+ Cube of angles of grid-x vector from true Eastward direction for
+ each gridcell, in degrees.
+ It also has "true" longitude and latitude coordinates, with no
+ coordinate system.
+ When the input has coords, then the output ones are identical if
+ the inputs are true-latlons, otherwise they are transformed
+ true-latlon versions.
+ When the input has bounded coords, then the output coords have
+ matching bounds and centrepoints (possibly transformed).
+ When the input is 2d arrays, or has unbounded coords, then the
+ output coords have matching points and no bounds.
+ When the input is 3d arrays, then the output coords have matching
+ bounds, and the centrepoints are an average of the 4 boundpoints.
+
+ """
+ cube = None
+ if hasattr(x, 'add_aux_coord'):
+ # Passed a cube : extract 'x' and ;'y' axis coordinates.
+ cube = x # Save for later checking.
+ x, y = cube.coord(axis='x'), cube.coord(axis='y')
+
+ # Now should have either 2 coords or 2 arrays.
+ if not hasattr(x, 'shape') or not hasattr(y, 'shape'):
+ msg = ('Inputs (x,y) must have array shape property.'
+ 'Got type(x)={} and type(y)={}.')
+ raise ValueError(msg.format(type(x), type(y)))
+
+ x_coord, y_coord = None, None
+ if hasattr(x, 'bounds') and hasattr(y, 'bounds'):
+ # x and y are Coords.
+ x_coord, y_coord = x.copy(), y.copy()
+
+ # They must be angles : convert into degrees
+ for coord in (x_coord, y_coord):
+ if not coord.units.is_convertible('degrees'):
+ msg = ('Input X and Y coordinates must have angular '
+ 'units. Got units of "{!s}" and "{!s}".')
+ raise ValueError(msg.format(x_coord.units, y_coord.units))
+ coord.convert_units('degrees')
+
+ if x_coord.ndim != 2 or y_coord.ndim != 2:
+ msg = ('Coordinate inputs must have 2-dimensional shape. '
+ 'Got x-shape of {} and y-shape of {}.')
+ raise ValueError(msg.format(x_coord.shape, y_coord.shape))
+ if x_coord.shape != y_coord.shape:
+ msg = ('Coordinate inputs must have same shape. '
+ 'Got x-shape of {} and y-shape of {}.')
+ raise ValueError(msg.format(x_coord.shape, y_coord.shape))
+ if cube:
+ x_dims, y_dims = (cube.coord_dims(co) for co in (x, y))
+ if x_dims != y_dims:
+ msg = ('X and Y coordinates must have the same cube '
+ 'dimensions. Got x-dims = {} and y-dims = {}.')
+ raise ValueError(msg.format(x_dims, y_dims))
+ cs = x_coord.coord_system
+ if y_coord.coord_system != cs:
+ msg = ('Coordinate inputs must have same coordinate system. '
+ 'Got x of {} and y of {}.')
+ raise ValueError(msg.format(cs, y_coord.coord_system))
+
+ # Base calculation on bounds if we have them, or points as a fallback.
+ if x_coord.has_bounds() and y_coord.has_bounds():
+ x, y = x_coord.bounds, y_coord.bounds
+ else:
+ x, y = x_coord.points, y_coord.points
+
+ # Make sure these arrays are ordinary lats+lons, in degrees.
+ if cs is not None:
+ # Transform points into true lats + lons.
+ crs_src = cs.as_cartopy_crs()
+ crs_pc = ccrs.PlateCarree()
+
+ def transform_xy_arrays(x, y):
+ # Note: flatten, as transform_points is limited to 2D arrays.
+ shape = x.shape
+ x, y = (arr.flatten() for arr in (x, y))
+ pts = crs_pc.transform_points(crs_src, x, y)
+ x = pts[..., 0].reshape(shape)
+ y = pts[..., 1].reshape(shape)
+ return x, y
+
+ # Transform the main reference points into standard lats+lons.
+ x, y = transform_xy_arrays(x, y)
+
+ # Likewise replace the original coordinates with transformed ones,
+ # because the calculation also needs the centrepoint values.
+ xpts, ypts = (coord.points for coord in (x_coord, y_coord))
+ xbds, ybds = (coord.bounds for coord in (x_coord, y_coord))
+ xpts, ypts = transform_xy_arrays(xpts, ypts)
+ xbds, ybds = transform_xy_arrays(xbds, ybds)
+ x_coord = iris.coords.AuxCoord(
+ points=xpts, bounds=xbds,
+ standard_name='longitude', units='degrees')
+ y_coord = iris.coords.AuxCoord(
+ points=ypts, bounds=ybds,
+ standard_name='latitude', units='degrees')
+
+ elif hasattr(x, 'bounds') or hasattr(y, 'bounds'):
+ # One was a Coord, and the other not ?
+ is_and_not = ('x', 'y')
+ if hasattr(y, 'bounds'):
+ is_and_not = reversed(is_and_not)
+ msg = 'Input {!r} is a Coordinate, but {!r} is not.'
+ raise ValueError(msg.format(*is_and_not))
+
+ # Now have either 2 points arrays (ny, nx) or 2 bounds arrays (ny, nx, 4).
+ # Construct (lhs, mid, rhs) where these represent 3 points at increasing
+ # grid-x indices (columns).
+ # Also make suitable X and Y coordinates for the result cube.
+ if x.ndim == 2:
+ # Data is points arrays.
+ # Use previous + subsequent points along grid-x-axis as references.
+
+ # PROBLEM: we assume that the rhs connects to the lhs, so we should
+ # really only use this if data is full-longitudes (as a 'circular'
+ # coordinate).
+ # This is mentioned in the docstring, but we have no general means of
+ # checking it.
+
+ # NOTE: we take the 2d grid as presented, so the second dimension is
+ # the 'X' dim. Again, that is implicit + can't be checked.
+ mid = np.array([x, y])
+ lhs = np.roll(mid, 1, 2)
+ rhs = np.roll(mid, -1, 2)
+ if not x_coord:
+ # Create coords for result cube : with no bounds.
+ y_coord = iris.coords.AuxCoord(x, standard_name='latitude',
+ units='degrees')
+ x_coord = iris.coords.AuxCoord(y, standard_name='longitude',
+ units='degrees')
+ else:
+ # Data is bounds arrays.
+ # Use gridcell corners at different grid-x positions as references.
+ # NOTE: so with bounds, we *don't* need full circular longitudes.
+ xyz = _3d_xyz_from_latlon(x, y)
+ # Support two different choices of reference points locations.
+ angle_boundpoints_vals = {'mid-lhs, mid-rhs': '03_to_12',
+ 'lower-left, lower-right': '0_to_1'}
+ bounds_pos = angle_boundpoints_vals.get(cell_angle_boundpoints)
+ if bounds_pos == '0_to_1':
+ lhs_xyz = xyz[..., 0]
+ rhs_xyz = xyz[..., 1]
+ elif bounds_pos == '03_to_12':
+ lhs_xyz = 0.5 * (xyz[..., 0] + xyz[..., 3])
+ rhs_xyz = 0.5 * (xyz[..., 1] + xyz[..., 2])
+ else:
+ msg = ('unrecognised cell_angle_boundpoints of "{}", '
+ 'must be one of {}')
+ raise ValueError(msg.format(cell_angle_boundpoints,
+ list(angle_boundpoints_vals.keys())))
+ if not x_coord:
+ # Create bounded coords for result cube.
+ # Use average of lhs+rhs points in 3d to get 'mid' points,
+ # as coords without points are not allowed.
+ mid_xyz = 0.5 * (lhs_xyz + rhs_xyz)
+ mid_latlons = _latlon_from_xyz(mid_xyz)
+ # Create coords with given bounds, and averaged centrepoints.
+ x_coord = iris.coords.AuxCoord(
+ points=mid_latlons[0], bounds=x,
+ standard_name='longitude', units='degrees')
+ y_coord = iris.coords.AuxCoord(
+ points=mid_latlons[1], bounds=y,
+ standard_name='latitude', units='degrees')
+
+ # Convert lhs and rhs points back to latlon form -- IN DEGREES !
+ lhs = _latlon_from_xyz(lhs_xyz)
+ rhs = _latlon_from_xyz(rhs_xyz)
+ # 'mid' is coord.points, whether from input or just made up.
+ mid = np.array([x_coord.points, y_coord.points])
+
+ # Do the angle calcs, and return as a suitable cube.
+ angles = _angle(lhs, mid, rhs)
+ result = iris.cube.Cube(angles,
+ long_name='gridcell_angle_from_true_east',
+ units='degrees')
+ result.add_aux_coord(x_coord, (0, 1))
+ result.add_aux_coord(y_coord, (0, 1))
+ return result
+
+
+def rotate_grid_vectors(u_cube, v_cube, grid_angles_cube=None,
+ grid_angles_kwargs=None):
+ """
+ Rotate distance vectors from grid-oriented to true-latlon-oriented.
+
+ Can also rotate by arbitrary angles, if they are passed in.
+
+ .. Note::
+
+ This operation overlaps somewhat in function with
+ :func:`iris.analysis.cartography.rotate_winds`.
+ However, that routine only rotates vectors according to transformations
+ between coordinate systems.
+ This function, by contrast, can rotate vectors by arbitrary angles.
+ Most commonly, the angles are estimated solely from grid sampling
+ points, using :func:`gridcell_angles` : This allows operation on
+ complex meshes defined by two-dimensional coordinates, such as most
+ ocean grids.
+
+ Args:
+
+ * u_cube, v_cube : (cube)
+ Cubes of grid-u and grid-v vector components.
+ Units should be differentials of true-distance, e.g. 'm/s'.
+
+ Optional args:
+
+ * grid_angles_cube : (cube)
+ gridcell orientation angles.
+ Units must be angular, i.e. can be converted to 'radians'.
+ If not provided, grid angles are estimated from 'u_cube' using the
+ :func:`gridcell_angles` method.
+
+ * grid_angles_kwargs : (dict or None)
+ Additional keyword args to be passed to the :func:`gridcell_angles`
+ method, if it is used.
+
+ Returns:
+
+ true_u, true_v : (cube)
+ Cubes of true-north oriented vector components.
+ Units are same as inputs.
+
+ .. Note::
+
+ Vector magnitudes will always be the same as the inputs.
+
+ """
+ u_out, v_out = (cube.copy() for cube in (u_cube, v_cube))
+ if not grid_angles_cube:
+ grid_angles_kwargs = grid_angles_kwargs or {}
+ grid_angles_cube = gridcell_angles(u_cube, **grid_angles_kwargs)
+ gridangles = grid_angles_cube.copy()
+ gridangles.convert_units('radians')
+ uu, vv, aa = (cube.data for cube in (u_out, v_out, gridangles))
+ mags = np.sqrt(uu*uu + vv*vv)
+ angs = np.arctan2(vv, uu) + aa
+ uu, vv = mags * np.cos(angs), mags * np.sin(angs)
+
+ # Promote all to masked arrays, and also apply mask at bad (NaN) angles.
+ mask = np.isnan(aa)
+ for cube in (u_out, v_out, aa):
+ if hasattr(cube.data, 'mask'):
+ mask |= cube.data.mask
+ u_out.data = np.ma.masked_array(uu, mask=mask)
+ v_out.data = np.ma.masked_array(vv, mask=mask)
+
+ return u_out, v_out
diff --git a/lib/iris/analysis/_regrid.py b/lib/iris/analysis/_regrid.py
index eb7a16d075..71304a5dde 100644
--- a/lib/iris/analysis/_regrid.py
+++ b/lib/iris/analysis/_regrid.py
@@ -231,8 +231,9 @@ def _regrid(src_data, x_dim, y_dim,
data = np.empty(shape, dtype=dtype)
# The interpolation class requires monotonically increasing
- # coordinates, so flip the coordinate(s) and data if the aren't.
- reverse_x = src_x_coord.points[0] > src_x_coord.points[1]
+ # coordinates, so flip the coordinate(s) and data if they aren't.
+ reverse_x = (src_x_coord.points[0] > src_x_coord.points[1] if
+ src_x_coord.points.size > 1 else False)
reverse_y = src_y_coord.points[0] > src_y_coord.points[1]
flip_index = [slice(None)] * src_data.ndim
if reverse_x:
diff --git a/lib/iris/analysis/cartography.py b/lib/iris/analysis/cartography.py
index 4b68dcc949..a778f5f846 100644
--- a/lib/iris/analysis/cartography.py
+++ b/lib/iris/analysis/cartography.py
@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2010 - 2017, Met Office
+# (C) British Crown Copyright 2010 - 2018, Met Office
#
# This file is part of Iris.
#
@@ -37,7 +37,24 @@
import iris.coord_systems
import iris.exceptions
from iris.util import _meshgrid
-
+from ._grid_angles import gridcell_angles, rotate_grid_vectors
+
+# List of contents to control Sphinx autodocs.
+# Unfortunately essential to get docs for the grid_angles functions.
+__all__ = [
+ 'area_weights',
+ 'cosine_latitude_weights',
+ 'get_xy_contiguous_bounded_grids',
+ 'get_xy_grids',
+ 'gridcell_angles',
+ 'project',
+ 'rotate_grid_vectors',
+ 'rotate_pole',
+ 'rotate_winds',
+ 'unrotate_pole',
+ 'wrap_lons',
+ 'DistanceDifferential',
+ 'PartialDifferential']
# This value is used as a fall-back if the cube does not define the earth
DEFAULT_SPHERICAL_EARTH_RADIUS = 6367470
@@ -332,7 +349,7 @@ def _quadrant_area(radian_lat_bounds, radian_lon_bounds, radius_of_earth):
def area_weights(cube, normalize=False):
- """
+ r"""
Returns an array of area weights, with the same dimensions as the cube.
This is a 2D lat/lon area weights array, repeated over the non lat/lon
@@ -440,7 +457,7 @@ def area_weights(cube, normalize=False):
def cosine_latitude_weights(cube):
- """
+ r"""
Returns an array of latitude weights, with the same dimensions as
the cube. The weights are the cosine of latitude.
@@ -950,7 +967,7 @@ def _transform_distance_vectors_tolerance_mask(src_crs, x, y, tgt_crs,
def rotate_winds(u_cube, v_cube, target_cs):
- """
+ r"""
Transform wind vectors to a different coordinate system.
The input cubes contain U and V components parallel to the local X and Y
diff --git a/lib/iris/analysis/geometry.py b/lib/iris/analysis/geometry.py
index a415e10091..cca5d836ec 100644
--- a/lib/iris/analysis/geometry.py
+++ b/lib/iris/analysis/geometry.py
@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2010 - 2015, Met Office
+# (C) British Crown Copyright 2010 - 2018, Met Office
#
# This file is part of Iris.
#
@@ -223,7 +223,7 @@ def geometry_area_weights(cube, geometry, normalize=False):
else:
slices.append(slice(None))
- weights[slices] = subweights
+ weights[tuple(slices)] = subweights
# Fix for the limitation of iris.analysis.MEAN weights handling.
# Broadcast the array to the full shape of the cube
diff --git a/lib/iris/coords.py b/lib/iris/coords.py
index e4a8eb35f2..ff2f170a2f 100644
--- a/lib/iris/coords.py
+++ b/lib/iris/coords.py
@@ -143,6 +143,85 @@ def __new__(cls, name_or_coord, minimum, maximum,
'groupby_point, groupby_slice')
+def _discontiguity_in_2d_bounds(bds, abs_tol=1e-4):
+ """
+ Check bounds of a 2-dimensional coordinate are contiguous
+ Args:
+ bds: Array of bounds of shape (X,Y,4)
+ abs_tol: tolerance
+
+ Returns:
+ Bool, if there are no discontinuities
+ absolute difference along the x axis
+ absolute difference along the y axis
+
+ """
+ # Check bds has the shape (ny, nx, 4)
+ if not bds.ndim == 3 and bds.shape[2] == 4:
+ raise ValueError('2D coordinates must have 4 bounds per point '
+ 'for 2D coordinate plotting')
+
+ # Check ordering:
+ # i i+1
+ # j @0 @1
+ # j+1 @3 @2
+ def mod360_diff(x1, x2):
+ diff = x1 - x2
+ diff = (diff + 360.0 + 180.0) % 360.0 - 180.0
+ return diff
+
+ # Compare cell with the cell next to it (i+1)
+ diffs_along_x = mod360_diff(bds[:, :-1, 1], bds[:, 1:, 0])
+ # Compare cell with the cell above it (j+1)
+ diffs_along_y = mod360_diff(bds[:-1, :, 3], bds[1:, :, 0])
+
+ def eq_diffs(x1):
+ return np.all(np.abs(x1) < abs_tol)
+
+ match_y0_x1 = eq_diffs(diffs_along_x)
+ match_y1_x0 = eq_diffs(diffs_along_y)
+
+ all_eq = match_y0_x1 and match_y1_x0
+
+ return all_eq, np.abs(diffs_along_x), np.abs(diffs_along_y)
+
+
+def _get_2d_coord_bound_grid(bds):
+ """
+ Function used that takes a bounds array for a 2-D coordinate variable with
+ 4 sides and returns the bounds grid.
+
+ Cf standards requires the four vertices of the cell to be traversed
+ anti-clockwise if the coordinates are defined in a right handed coordinate
+ system.
+
+ selects the zeroth vertex of each cell and then adds the column the first
+ vertex at the end. For the top row it uses the thirs vertex, with the
+ second added on to the end.
+ e.g.
+ # 0-0-0-0-1
+ # 0-0-0-0-1
+ # 3-3-3-3-2
+
+
+ Args:
+ bounds: array of shape (X,Y,4)
+
+ Returns:
+ array of shape (X+1, Y+1)
+
+ """
+ bds_shape = bds.shape
+ result = np.zeros((bds_shape[0] + 1, bds_shape[1] + 1))
+
+ result[:-1, :-1] = bds[:, :, 0]
+ result[:-1, -1] = bds[:, -1, 1]
+ result[-1, :-1] = bds[-1, :, 3]
+ result[-1, -1] = bds[-1, -1, 2]
+
+ return result
+
+
class Cell(collections.namedtuple('Cell', ['point', 'bound'])):
"""
An immutable representation of a single cell of a coordinate, including the
@@ -951,15 +1030,22 @@ def cells(self):
return _CellIterator(self)
def _sanity_check_contiguous(self):
- if self.ndim != 1:
- raise iris.exceptions.CoordinateMultiDimError(
- 'Invalid operation for {!r}. Contiguous bounds are not defined'
- ' for multi-dimensional coordinates.'.format(self.name()))
- if self.nbounds != 2:
- raise ValueError(
- 'Invalid operation for {!r}, with {} bounds. Contiguous bounds'
- ' are only defined for coordinates with 2 bounds.'.format(
- self.name(), self.nbounds))
+ if self.ndim == 1:
+ if self.nbounds != 2:
+ raise ValueError('Invalid operation for {!r}, with {} bounds. '
+ 'Contiguous bounds are only defined for 1D '
+ 'coordinates with 2 bounds.'.format
+ (self.name(), self.nbounds))
+ elif self.ndim == 2:
+ if self.nbounds != 4:
+ raise ValueError('Invalid operation for {!r}, with {} bounds. '
+ 'Contiguous bounds are only defined for 2D '
+ 'coordinates with 4 bounds.'.format
+ (self.name(), self.nbounds))
+ else:
+ raise ValueError('Invalid operation for {!r}. Contiguous bounds '
+ 'are not defined for coordinates with more than '
+ '2 dimensions.'.format(self.name()))
def is_contiguous(self, rtol=1e-05, atol=1e-08):
"""
@@ -979,19 +1065,30 @@ def is_contiguous(self, rtol=1e-05, atol=1e-08):
"""
if self.has_bounds():
self._sanity_check_contiguous()
- return np.allclose(self.bounds[1:, 0], self.bounds[:-1, 1],
- rtol=rtol, atol=atol)
+ if self.ndim == 1:
+ contiguous = np.allclose(self.bounds[1:, 0],
+ self.bounds[:-1, 1],
+ rtol=rtol, atol=atol)
+ elif self.ndim == 2:
+ contiguous, _, _ = _discontiguity_in_2d_bounds(self.bounds,
+ abs_tol=atol)
else:
- return False
+ contiguous = False
+ return contiguous
def contiguous_bounds(self):
"""
- Returns the N+1 bound values for a contiguous bounded coordinate
+ Returns the N+1 bound values for a contiguous bounded 1D coordinate
of length N.
+ Returns the (N+1, M+1) bound values for a contiguous bounded 2D
+ coordinate of shape (N, M).
+
+ Assumes input is contiguous.
+
.. note::
- If the coordinate is does not have bounds, this method will
+ If the coordinate does not have bounds, this method will
return bounds positioned halfway between the coordinate's points.
"""
@@ -1003,8 +1100,11 @@ def contiguous_bounds(self):
self._sanity_check_contiguous()
bounds = self.bounds
- c_bounds = np.resize(bounds[:, 0], bounds.shape[0] + 1)
- c_bounds[-1] = bounds[-1, 1]
+ if self.ndim == 1:
+ c_bounds = np.resize(bounds[:, 0], bounds.shape[0] + 1)
+ c_bounds[-1] = bounds[-1, 1]
+ elif self.ndim == 2:
+ c_bounds = _get_2d_coord_bound_grid(bounds)
return c_bounds
def is_monotonic(self):
diff --git a/lib/iris/pandas.py b/lib/iris/pandas.py
index b4e1e9787a..244a7c011b 100644
--- a/lib/iris/pandas.py
+++ b/lib/iris/pandas.py
@@ -130,18 +130,7 @@ def _as_pandas_coord(coord):
def _assert_shared(np_obj, pandas_obj):
"""Ensure the pandas object shares memory."""
- if hasattr(pandas_obj, 'base'):
- base = pandas_obj.base
- else:
- base = pandas_obj[0].base
-
- # Prior to Pandas 0.17, when pandas_obj is a Series, pandas_obj.values
- # returns a view of the underlying array, and pandas_obj.base, which calls
- # pandas_obj.values.base, returns the underlying array. In 0.17 and 0.18
- # pandas_obj.values returns the underlying array, so base may be None even
- # if the array is shared.
- if base is None:
- base = pandas_obj.values
+ values = pandas_obj.values
def _get_base(array):
# Chase the stack of NumPy `base` references back to the original array
@@ -149,7 +138,7 @@ def _get_base(array):
array = array.base
return array
- base = _get_base(base)
+ base = _get_base(values)
np_base = _get_base(np_obj)
if base is not np_base:
msg = ('Pandas {} does not share memory'
diff --git a/lib/iris/plot.py b/lib/iris/plot.py
index 64ce54c04c..7bf8ea07f3 100644
--- a/lib/iris/plot.py
+++ b/lib/iris/plot.py
@@ -31,13 +31,13 @@
import cartopy.crs as ccrs
import cartopy.mpl.geoaxes
from cartopy.geodesic import Geodesic
+import cftime
import matplotlib.axes
import matplotlib.collections as mpl_collections
import matplotlib.dates as mpl_dates
import matplotlib.pyplot as plt
from matplotlib.offsetbox import AnchoredText
import matplotlib.ticker as mpl_ticker
-import cftime
import matplotlib.transforms as mpl_transforms
import numpy as np
import numpy.ma as ma
@@ -51,7 +51,6 @@
import iris.palette
from iris.util import _meshgrid
-
# Cynthia Brewer citation text.
BREWER_CITE = 'Colours based on ColorBrewer.org'
@@ -95,15 +94,17 @@ def get_span(coord):
else:
span = set(cube.coord_dims(coord))
return span
+
spans = list(map(get_span, coords))
for span, coord in zip(spans, coords):
if not span:
msg = 'The coordinate {!r} doesn\'t span a data dimension.'
raise ValueError(msg.format(coord.name()))
- if mode == iris.coords.BOUND_MODE and len(span) != 1:
- raise ValueError('The coordinate {!r} is multi-dimensional and'
- ' cannot be used in a cell-based plot.'
- .format(coord.name()))
+ if mode == iris.coords.BOUND_MODE and len(span) not in [1, 2]:
+ raise ValueError('The coordinate {!r} has {} dimensions.'
+ 'Cell-based plotting is only supporting for'
+ 'coordinates with one or two dimensions.'
+ .format(coord.name()), len(span))
# Check the combination of coordinates spans enough (ndims) data
# dimensions.
@@ -129,7 +130,7 @@ def get_span(coord):
return PlotDefn(plot_coords, transpose)
-def _valid_bound_coord(coord):
+def _valid_bound_dim_coord(coord):
result = None
if coord and coord.ndim == 1 and coord.nbounds:
result = coord
@@ -154,7 +155,7 @@ def _get_plot_defn(cube, mode, ndims=2):
# When appropriate, restrict to 1D with bounds.
if mode == iris.coords.BOUND_MODE:
- coords = list(map(_valid_bound_coord, coords))
+ coords = list(map(_valid_bound_dim_coord, coords))
def guess_axis(coord):
axis = None
@@ -172,6 +173,19 @@ def guess_axis(coord):
aux_coords.sort(key=lambda coord: coord._as_defn())
coords[dim] = aux_coords[0]
+ # If plotting a 2 dimensional plot, check for 2d coordinates
+ if ndims == 2:
+ missing_dims = [dim for dim, coord in enumerate(coords) if
+ coord is None]
+ if missing_dims:
+ # Note that this only picks up coordinates that span the dims
+ two_dim_coords = cube.coords(dimensions=missing_dims)
+ two_dim_coords = [coord for coord in two_dim_coords
+ if coord.ndim == 2]
+ if len(two_dim_coords) >= 2:
+ two_dim_coords.sort(key=lambda coord: coord._as_defn())
+ coords = two_dim_coords[:2]
+
if mode == iris.coords.POINT_MODE:
# Allow multi-dimensional aux_coords to override the dim_coords
# along the Z axis. This results in a preference for using the
@@ -195,6 +209,7 @@ def sort_key(coord):
return (order.get(axis, 0),
coords.index(coord),
coord and coord.name())
+
sorted_coords = sorted(coords, key=sort_key)
transpose = (sorted_coords != coords)
@@ -256,6 +271,56 @@ def _invert_yaxis(v_coord, axes=None):
axes.invert_yaxis()
+def _check_contiguity_and_bounds(coord, data, abs_tol=1e-4, transpose=False):
+ """
+ Check that the discontinuous bounds occur where the data is masked.
+
+ If discontinuity occurs but data is masked, raise warning
+ If discontinuity occurs and data is NOT masked, raise error
+
+ Args:
+ coords:
+ Array of the bounds of a 2D coord, of shape (X,Y,4)
+ data:
+ data of the the cube we are plotting
+ abs_tol:
+ tolerance when checking the contiguity
+
+ """
+ if transpose:
+ data = data.T
+
+ bounds = coord.bounds
+
+ both_dirs_contiguous, diffs_along_x, diffs_along_y = \
+ iris.coords._discontiguity_in_2d_bounds(bounds, abs_tol=abs_tol)
+
+ if not both_dirs_contiguous:
+
+ # True where data exists
+ mask_invert = np.logical_not(data.mask)
+
+ # Where a discontinuity occurs the grid will not be created correctly.
+ # This does not matter if the data is masked as this is not plotted.
+ # So for places where data exists (opposite of the mask) AND a
+ # diff exists. If any exist, raise a warning
+ not_masked_at_discontinuity_along_x = np.any(
+ np.logical_and(mask_invert[:, :-1], diffs_along_x))
+
+ not_masked_at_discontinuity_along_y = np.any(
+ np.logical_and(mask_invert[:-1, ], diffs_along_y))
+
+ # If discontinuity occurs but not masked, any grid will be created
+ # incorrectly, so raise a warning
+ if not_masked_at_discontinuity_along_x or \
+ not_masked_at_discontinuity_along_y:
+ raise ValueError('The bounds of the {} coordinate are not'
+ ' contiguous and data is not masked where the'
+ ' discontiguity occurs. Not able to create a'
+ ' suitable mesh to give to'
+ ' Matplotlib'.format(coord.name()))
+
+
def _draw_2d_from_bounds(draw_method_name, cube, *args, **kwargs):
# NB. In the interests of clarity we use "u" and "v" to refer to the
# horizontal and vertical axes on the matplotlib plot.
@@ -263,10 +328,25 @@ def _draw_2d_from_bounds(draw_method_name, cube, *args, **kwargs):
# Get & remove the coords entry from kwargs.
coords = kwargs.pop('coords', None)
if coords is not None:
- plot_defn = _get_plot_defn_custom_coords_picked(cube, coords, mode)
+ plot_defn = _get_plot_defn_custom_coords_picked(cube, coords, mode,
+ ndims=2)
else:
plot_defn = _get_plot_defn(cube, mode, ndims=2)
+ twodim_contig_atol = kwargs.pop('two_dim_coord_contiguity_atol',
+ 1e-4)
+ for coord in plot_defn.coords:
+ if hasattr(coord, 'has_bounds'):
+ if coord.ndim == 2 and coord.has_bounds():
+ try:
+ _check_contiguity_and_bounds(coord, data=cube.data,
+ abs_tol=twodim_contig_atol)
+ except ValueError:
+ if _check_contiguity_and_bounds(coord, data=cube.data,
+ abs_tol=twodim_contig_atol,
+ transpose=True) is True:
+ plot_defn.transpose = True
+
if _can_draw_map(plot_defn.coords):
result = _map_common(draw_method_name, None, iris.coords.BOUND_MODE,
cube, plot_defn, *args, **kwargs)
@@ -309,7 +389,16 @@ def _draw_2d_from_bounds(draw_method_name, cube, *args, **kwargs):
plot_arrays.append(values)
u, v = plot_arrays
- u, v = _broadcast_2d(u, v)
+
+ # If the data is tranposed, 2D coordinates will also need to be
+ # tranposed.
+ if u.ndim == v.ndim == 2 and plot_defn.transpose is True:
+ u = u.T
+ v = v.T
+
+ if u.ndim == v.ndim == 1:
+ u, v = _broadcast_2d(u, v)
+
axes = kwargs.pop('axes', None)
draw_method = getattr(axes if axes else plt, draw_method_name)
result = draw_method(u, v, data, *args, **kwargs)
@@ -434,8 +523,8 @@ def _fixup_dates(coord, values):
raise IrisError(msg)
r = [nc_time_axis.CalendarDateTime(
- cftime.datetime(*date), coord.units.calendar)
- for date in dates]
+ cftime.datetime(*date), coord.units.calendar)
+ for date in dates]
values = np.empty(len(r), dtype=object)
values[:] = r
return values
@@ -675,8 +764,8 @@ def _ensure_cartopy_axes_and_determine_kwargs(x_coord, y_coord, kwargs):
axes = kwargs.get('axes')
if axes is None:
if (isinstance(cs, iris.coord_systems.RotatedGeogCS) and
- x_coord.points.max() > 180 and x_coord.points.max() < 360 and
- x_coord.points.min() > 0):
+ x_coord.points.max() > 180 and x_coord.points.max() < 360 and
+ x_coord.points.min() > 0):
# The RotatedGeogCS has 0 - 360 extent, different from the
# assumptions made by Cartopy: rebase longitudes for the map axes
# to set the datum longitude to the International Date Line.
@@ -721,17 +810,22 @@ def _map_common(draw_method_name, arg_func, mode, cube, plot_defn,
else:
raise ValueError("Expected 1D or 2D XY coords")
else:
- try:
- x, y = _meshgrid(x_coord.contiguous_bounds(),
- y_coord.contiguous_bounds())
- # Exception translation.
- except iris.exceptions.CoordinateMultiDimError:
- raise ValueError("Could not get XY grid from bounds. "
- "X or Y coordinate not 1D.")
- except ValueError:
- raise ValueError("Could not get XY grid from bounds. "
- "X or Y coordinate doesn't have 2 bounds "
- "per point.")
+ if not x_coord.ndim == y_coord.ndim == 2:
+ try:
+ x, y = _meshgrid(x_coord.contiguous_bounds(),
+ y_coord.contiguous_bounds())
+ # Exception translation.
+ except iris.exceptions.CoordinateMultiDimError:
+ raise ValueError("Expected two 1D coords. Could not get XY"
+ " grid from bounds. X or Y coordinate not"
+ " 1D.")
+ except ValueError:
+ raise ValueError("Could not get XY grid from bounds. "
+ "X or Y coordinate doesn't have 2 bounds "
+ "per point.")
+ else:
+ x = x_coord.contiguous_bounds()
+ y = y_coord.contiguous_bounds()
# Obtain the data array.
data = cube.data
@@ -1005,7 +1099,11 @@ def outline(cube, coords=None, color='k', linewidth=None, axes=None):
def pcolor(cube, *args, **kwargs):
"""
- Draws a pseudocolor plot based on the given Cube.
+ Draws a pseudocolor plot based on the given 2-dimensional Cube.
+
+ The cube must have either two 1-dimensional coordinates or two
+ 2-dimensional coordinates with contiguous bounds to plot against each
+ other.
Kwargs:
@@ -1019,6 +1117,11 @@ def pcolor(cube, *args, **kwargs):
* axes: the :class:`matplotlib.axes.Axes` to use for drawing.
Defaults to the current axes if none provided.
+ * two_dim_coord_contiguity_atol: absolute tolerance when checking for
+ contiguity between cells in a two dimensional coordinate.
+
+
+
See :func:`matplotlib.pyplot.pcolor` for details of other valid
keyword arguments.
@@ -1031,7 +1134,11 @@ def pcolor(cube, *args, **kwargs):
def pcolormesh(cube, *args, **kwargs):
"""
- Draws a pseudocolor plot based on the given Cube.
+ Draws a pseudocolor plot based on the given 2-dimensional Cube.
+
+ The cube must have either two 1-dimensional coordinates or two
+ 2-dimensional coordinates with contiguous bounds to plot against each
+ other.
Kwargs:
@@ -1045,6 +1152,9 @@ def pcolormesh(cube, *args, **kwargs):
* axes: the :class:`matplotlib.axes.Axes` to use for drawing.
Defaults to the current axes if none provided.
+ * two_dim_coord_contiguity_atol: absolute tolerance when checking for
+ contiguity between cells in a two dimensional coordinate.
+
See :func:`matplotlib.pyplot.pcolormesh` for details of other
valid keyword arguments.
@@ -1073,6 +1183,7 @@ def points(cube, *args, **kwargs):
keyword arguments.
"""
+
def _scatter_args(u, v, data, *args, **kwargs):
return ((u, v) + args, kwargs)
@@ -1080,6 +1191,139 @@ def _scatter_args(u, v, data, *args, **kwargs):
*args, **kwargs)
+def _vector_component_args(x_points, y_points, u_data, *args, **kwargs):
+ """
+ Callback from _draw_2d_from_points for 'quiver' and 'streamlines'.
+
+ Returns arguments (x, y, u, v), to be passed to the underlying matplotlib
+ call.
+
+ "u_data" will always be "u_cube.data".
+ The matching "v_cube.data" component is stored in kwargs['_v_data'].
+
+ """
+ v_data = kwargs.pop('_v_data')
+
+ # Rescale u+v values for plot distortion.
+ crs = kwargs.get('transform', None)
+ if crs:
+ if not isinstance(crs, (ccrs.PlateCarree, ccrs.RotatedPole)):
+ msg = ('Can only plot vectors provided in a lat-lon '
+ 'projection, i.e. "cartopy.crs.PlateCarree" or '
+ '"cartopy.crs.RotatedPole". This '
+ "cubes coordinate system is {}.")
+ raise ValueError(msg.format(crs))
+ # Given the above check, the Y points must be latitudes.
+ # We therefore **assume** they are in degrees : I'm not sure this
+ # is wise, but all the rest of this plot code does that, e.g. in
+ # _map_common.
+ # TODO: investigate degree units assumptions, here + elsewhere.
+
+ # Implement a latitude scaling, but preserve the given magnitudes.
+ v_data = v_data.copy()
+ mags = np.sqrt(u_data * u_data + v_data * v_data)
+ v_data *= np.cos(np.deg2rad(y_points))
+ scales = mags / np.sqrt(u_data * u_data + v_data * v_data)
+ u_data *= scales
+ v_data *= scales
+
+ return ((x_points, y_points, u_data, v_data), kwargs)
+
+
+def quiver(u_cube, v_cube, *args, **kwargs):
+ """
+ Draws an arrow plot from two vector component cubes.
+
+ Args:
+
+ * u_cube, v_cube : (:class:`~iris.cube.Cube`)
+ u and v vector components. Must have same shape and units.
+ If the cubes have geographic coordinates, the values are treated as
+ true distance differentials, e.g. windspeeds, and *not* map coordinate
+ vectors. The components are aligned with the North and East of the
+ cube coordinate system.
+
+ .. Note:
+
+ At present, if u_cube and v_cube have geographic coordinates, then they
+ must be in a lat-lon coordinate system, though it may be a rotated one.
+ To transform wind values between coordinate systems, use
+ :func:`iris.analysis.cartography.rotate_vectors`.
+ To transform coordinate grid points, you will need to create
+ 2-dimensional arrays of x and y values. These can be transformed with
+ :meth:`cartopy.crs.CRS.transform_points`.
+
+ Kwargs:
+
+ * coords: (list of :class:`~iris.coords.Coord` or string)
+ Coordinates or coordinate names. Use the given coordinates as the axes
+ for the plot. The order of the given coordinates indicates which axis
+ to use for each, where the first element is the horizontal
+ axis of the plot and the second element is the vertical axis
+ of the plot.
+
+ * axes: the :class:`matplotlib.axes.Axes` to use for drawing.
+ Defaults to the current axes if none provided.
+
+ See :func:`matplotlib.pyplot.quiver` for details of other valid
+ keyword arguments.
+
+ """
+ #
+ # TODO: check u + v cubes for compatibility.
+ #
+ kwargs['_v_data'] = v_cube.data
+ return _draw_2d_from_points('quiver', _vector_component_args, u_cube,
+ *args, **kwargs)
+
+
+def streamplot(u_cube, v_cube, *args, **kwargs):
+ """
+ Draws a streamline plot from two vector component cubes.
+
+ Args:
+
+ * u_cube, v_cube : (:class:`~iris.cube.Cube`)
+ u and v vector components. Must have same shape and units.
+ If the cubes have geographic coordinates, the values are treated as
+ true distance differentials, e.g. windspeeds, and *not* map coordinate
+ vectors. The components are aligned with the North and East of the
+ cube coordinate system.
+
+ .. Note:
+
+ At present, if u_cube and v_cube have geographic coordinates, then they
+ must be in a lat-lon coordinate system, though it may be a rotated one.
+ To transform wind values between coordinate systems, use
+ :func:`iris.analysis.cartography.rotate_vectors`.
+ To transform coordinate grid points, you will need to create
+ 2-dimensional arrays of x and y values. These can be transformed with
+ :meth:`cartopy.crs.CRS.transform_points`.
+
+ Kwargs:
+
+ * coords: (list of :class:`~iris.coords.Coord` or string)
+ Coordinates or coordinate names. Use the given coordinates as the axes
+ for the plot. The order of the given coordinates indicates which axis
+ to use for each, where the first element is the horizontal
+ axis of the plot and the second element is the vertical axis
+ of the plot.
+
+ * axes: the :class:`matplotlib.axes.Axes` to use for drawing.
+ Defaults to the current axes if none provided.
+
+ See :func:`matplotlib.pyplot.quiver` for details of other valid
+ keyword arguments.
+
+ """
+ #
+ # TODO: check u + v cubes for compatibility.
+ #
+ kwargs['_v_data'] = v_cube.data
+ return _draw_2d_from_points('streamplot', _vector_component_args, u_cube,
+ *args, **kwargs)
+
+
def plot(*args, **kwargs):
"""
Draws a line plot based on the given cube(s) or coordinate(s).
diff --git a/lib/iris/tests/__init__.py b/lib/iris/tests/__init__.py
index 8c64dc15b1..7550fce1ca 100644
--- a/lib/iris/tests/__init__.py
+++ b/lib/iris/tests/__init__.py
@@ -663,7 +663,23 @@ def assertArrayAllClose(self, a, b, rtol=1.0e-7, atol=0.0, **kwargs):
For full details see underlying routine numpy.testing.assert_allclose.
"""
- np.testing.assert_allclose(a, b, rtol=rtol, atol=atol, **kwargs)
+ ok = np.allclose(a, b, atol=atol, rtol=rtol, equal_nan=True)
+ if not ok:
+ # Calculate errors above a pointwise tolerance : The method is
+ # taken from "numpy.core.numeric.isclose".
+ errors = (np.abs(a-b) - atol + rtol * np.abs(b))
+ worst_inds = np.unravel_index(np.argmax(errors.flat), errors.shape)
+ # Build a more useful message than from np.testing.assert_allclose.
+ msg = ('\nARRAY CHECK FAILED "assertArrayAllClose" :'
+ '\n with shapes={} {}, atol={}, rtol={}'
+ '\n worst at element {} : a={} b={}'
+ '\n absolute error ~{:.3g}, equivalent to rtol ~{:.3e}')
+ aval, bval = a[worst_inds], b[worst_inds]
+ absdiff = np.abs(aval - bval)
+ equiv_rtol = absdiff / bval
+ raise AssertionError(msg.format(
+ a.shape, b.shape, atol, rtol, worst_inds, aval, bval, absdiff,
+ equiv_rtol))
@contextlib.contextmanager
def temp_filename(self, suffix=''):
@@ -1196,6 +1212,12 @@ class MyPlotTests(test.GraphicsTest):
'Test(s) require "python-stratify", which is not available.')
+SKIP_2D_TESTS = True
+skip_2d = unittest.skipIf(
+ SKIP_2D_TESTS,
+ 'Test(s) broken by WIP on 2d coords support -- temporarily disabled.')
+
+
def no_warnings(func):
"""
Provides a decorator to ensure that there are no warnings raised
diff --git a/lib/iris/tests/integration/test_regridding.py b/lib/iris/tests/integration/test_regridding.py
index ef521aa395..47d155f3f2 100644
--- a/lib/iris/tests/integration/test_regridding.py
+++ b/lib/iris/tests/integration/test_regridding.py
@@ -109,5 +109,116 @@ def test_nearest(self):
self.assertArrayShapeStats(res, (1, 6, 3, 4), 315.890808, 11.000724)
+class TestZonalMean_global(tests.IrisTest):
+ def setUp(self):
+ np.random.seed(0)
+ self.src = iris.cube.Cube(np.random.random_integers(0, 10, (140, 1)))
+ s_crs = iris.coord_systems.GeogCS(6371229.0)
+ sy_coord = iris.coords.DimCoord(
+ np.linspace(-90, 90, 140), standard_name='latitude',
+ units='degrees', coord_system=s_crs)
+ sx_coord = iris.coords.DimCoord(
+ -180, bounds=[-180, 180], standard_name='longitude',
+ units='degrees', circular=True, coord_system=s_crs)
+ self.src.add_dim_coord(sy_coord, 0)
+ self.src.add_dim_coord(sx_coord, 1)
+
+ def test_linear_same_crs_global(self):
+ # Regrid the zonal mean onto an identical coordinate system target, but
+ # on a different set of longitudes - which should result in no change.
+ points = [-150, -90, -30, 30, 90, 150]
+ bounds = [[-180, -120], [-120, -60], [-60, 0], [0, 60], [60, 120],
+ [120, 180]]
+ sx_coord = self.src.coord(axis='x')
+ sy_coord = self.src.coord(axis='y')
+ x_coord = sx_coord.copy(points, bounds=bounds)
+ grid = iris.cube.Cube(
+ np.zeros([sy_coord.points.size, x_coord.points.size]))
+ grid.add_dim_coord(sy_coord, 0)
+ grid.add_dim_coord(x_coord, 1)
+
+ res = self.src.regrid(grid, iris.analysis.Linear())
+
+ # Ensure data remains unchanged.
+ # (the same along each column)
+ self.assertTrue(
+ np.array([(res.data[:, 0]-res.data[:, i]).max() for i in
+ range(1, res.shape[1])]).max() < 1e-10)
+ self.assertArrayAlmostEqual(res.data[:, 0], self.src.data.reshape(-1))
+
+
+class TestZonalMean_regional(TestZonalMean_global, tests.IrisTest):
+ def setUp(self):
+ super(TestZonalMean_regional, self).setUp()
+
+ # Define a target grid and a target result (what we expect the
+ # regridder to return).
+ sx_coord = self.src.coord(axis='x')
+ sy_coord = self.src.coord(axis='y')
+ grid_crs = iris.coord_systems.RotatedGeogCS(
+ 37.5, 177.5, ellipsoid=iris.coord_systems.GeogCS(6371229.0))
+ grid_x = sx_coord.copy(np.linspace(350, 370, 100))
+ grid_x.circular = False
+ grid_x.coord_system = grid_crs
+ grid_y = sy_coord.copy(np.linspace(-10, 10, 100))
+ grid_y.coord_system = grid_crs
+ grid = iris.cube.Cube(
+ np.zeros([grid_y.points.size, grid_x.points.size]))
+ grid.add_dim_coord(grid_y, 0)
+ grid.add_dim_coord(grid_x, 1)
+
+ # The target result is derived by regridding a multi-column version of
+ # the source to the target (i.e. turning a zonal mean regrid into a
+ # conventional regrid).
+ self.tar = self.zonal_mean_as_multi_column(self.src).regrid(
+ grid, iris.analysis.Linear())
+ self.grid = grid
+
+ def zonal_mean_as_multi_column(self, src_cube):
+ # Munge the source (duplicate source latitudes) so that we can
+ # utilise linear regridding as a conventional problem (that is, to
+ # duplicate columns so that it is no longer a zonal mean problem).
+ src_cube2 = src_cube.copy()
+ src_cube2.coord(axis='x').points = -90
+ src_cube2.coord(axis='x').bounds = [-180, 0]
+ src_cube.coord(axis='x').points = 90
+ src_cube.coord(axis='x').bounds = [0, 180]
+ src_cubes = iris.cube.CubeList([src_cube, src_cube2])
+ return src_cubes.concatenate_cube()
+
+ def test_linear_rotated_regional(self):
+ # Ensure that zonal mean source data is linearly interpolated onto a
+ # high resolution target.
+ regridder = iris.analysis.Linear()
+ res = self.src.regrid(self.grid, regridder)
+ self.assertArrayAlmostEqual(res.data, self.tar.data)
+
+ def test_linear_rotated_regional_no_extrapolation(self):
+ # Capture the case where our source remains circular but we don't use
+ # extrapolation.
+ regridder = iris.analysis.Linear(extrapolation_mode='nan')
+ res = self.src.regrid(self.grid, regridder)
+ self.assertArrayAlmostEqual(res.data, self.tar.data)
+
+ def test_linear_rotated_regional_not_circular(self):
+ # Capture the case where our source is not circular but we utilise
+ # extrapolation.
+ regridder = iris.analysis.Linear()
+ self.src.coord(axis='x').circular = False
+ res = self.src.regrid(self.grid, regridder)
+ self.assertArrayAlmostEqual(res.data, self.tar.data)
+
+ def test_linear_rotated_regional_no_extrapolation_not_circular(self):
+ # Confirm how zonal mean actually works in so far as, that
+ # extrapolation and circular source handling is the means by which
+ # these usecases are supported.
+ # In the case where the source is neither using extrapolation and is
+ # not circular, then 'nan' values will result (as we would expect).
+ regridder = iris.analysis.Linear(extrapolation_mode='nan')
+ self.src.coord(axis='x').circular = False
+ res = self.src.regrid(self.grid, regridder)
+ self.assertTrue(np.isnan(res.data).all())
+
+
if __name__ == "__main__":
tests.main()
diff --git a/lib/iris/tests/stock.py b/lib/iris/tests/stock/__init__.py
similarity index 74%
rename from lib/iris/tests/stock.py
rename to lib/iris/tests/stock/__init__.py
index 25df8eda33..5965d5a208 100644
--- a/lib/iris/tests/stock.py
+++ b/lib/iris/tests/stock/__init__.py
@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2010 - 2017, Met Office
+# (C) British Crown Copyright 2010 - 2018, Met Office
#
# This file is part of Iris.
#
@@ -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():
"""
@@ -46,12 +47,12 @@ def lat_lon_cube():
"""
cube = Cube(np.arange(12, dtype=np.int32).reshape((3, 4)))
cs = GeogCS(6371229)
- coord = iris.coords.DimCoord(points=np.array([-1, 0, 1], dtype=np.int32),
+ coord = DimCoord(points=np.array([-1, 0, 1], dtype=np.int32),
standard_name='latitude',
units='degrees',
coord_system=cs)
cube.add_dim_coord(coord, 0)
- coord = iris.coords.DimCoord(points=np.array([-1, 0, 1, 2], dtype=np.int32),
+ coord = DimCoord(points=np.array([-1, 0, 1, 2], dtype=np.int32),
standard_name='longitude',
units='degrees',
coord_system=cs)
@@ -99,7 +100,7 @@ def simple_1d(with_bounds=True):
cube.units = '1'
points = np.arange(11, dtype=np.int32) + 1
bounds = np.column_stack([np.arange(11, dtype=np.int32), np.arange(11, dtype=np.int32) + 1])
- coord = iris.coords.DimCoord(points, long_name='foo', units='1', bounds=bounds)
+ coord = DimCoord(points, long_name='foo', units='1', bounds=bounds)
cube.add_dim_coord(coord, 0)
return cube
@@ -124,12 +125,15 @@ def simple_2d(with_bounds=True):
cube = Cube(np.arange(12, dtype=np.int32).reshape((3, 4)))
cube.long_name = 'thingness'
cube.units = '1'
- y_points = np.array([ 2.5, 7.5, 12.5])
+ y_points = np.array([2.5, 7.5, 12.5])
y_bounds = np.array([[0, 5], [5, 10], [10, 15]], dtype=np.int32)
- y_coord = iris.coords.DimCoord(y_points, long_name='bar', units='1', bounds=y_bounds if with_bounds else None)
+ y_coord = DimCoord(y_points, long_name='bar', units='1',
+ bounds=y_bounds if with_bounds else None)
x_points = np.array([ -7.5, 7.5, 22.5, 37.5])
- x_bounds = np.array([[-15, 0], [0, 15], [15, 30], [30, 45]], dtype=np.int32)
- x_coord = iris.coords.DimCoord(x_points, long_name='foo', units='1', bounds=x_bounds if with_bounds else None)
+ x_bounds = np.array([[-15, 0], [0, 15], [15, 30], [30, 45]],
+ dtype=np.int32)
+ x_coord = DimCoord(x_points, long_name='foo', units='1',
+ bounds=x_bounds if with_bounds else None)
cube.add_dim_coord(y_coord, 0)
cube.add_dim_coord(x_coord, 1)
@@ -191,9 +195,8 @@ def simple_3d_w_multidim_coords(with_bounds=True):
[[5, 15], [15, 20], [20, 35], [35, 50]],
[[10, 20], [20, 25], [25, 40], [40, 60]]],
dtype=np.int32)
- y_coord = iris.coords.AuxCoord(points=y_points, long_name='bar',
- units='1',
- bounds=y_bounds if with_bounds else None)
+ y_coord = AuxCoord(points=y_points, long_name='bar', units='1',
+ bounds=y_bounds if with_bounds else None)
x_points = np.array([[-7.5, 7.5, 22.5, 37.5],
[-12.5, 4., 26.5, 47.5],
[2.5, 14., 36.5, 44.]])
@@ -201,12 +204,10 @@ def simple_3d_w_multidim_coords(with_bounds=True):
[[-25, 0], [0, 8], [8, 45], [45, 50]],
[[-5, 10], [10, 18], [18, 55], [18, 70]]],
dtype=np.int32)
- x_coord = iris.coords.AuxCoord(points=x_points, long_name='foo',
- units='1',
- bounds=x_bounds if with_bounds else None)
- wibble_coord = iris.coords.DimCoord(np.array([10., 30.],
- dtype=np.float32),
- long_name='wibble', units='1')
+ x_coord = AuxCoord(points=x_points, long_name='foo', units='1',
+ bounds=x_bounds if with_bounds else None)
+ wibble_coord = DimCoord(np.array([10., 30.], dtype=np.float32),
+ long_name='wibble', units='1')
cube.add_dim_coord(wibble_coord, [0])
cube.add_aux_coord(y_coord, [1, 2])
@@ -238,13 +239,13 @@ def simple_3d():
cube = Cube(np.arange(24, dtype=np.int32).reshape((2, 3, 4)))
cube.long_name = 'thingness'
cube.units = '1'
- wibble_coord = iris.coords.DimCoord(np.array([10., 30.],
+ wibble_coord = DimCoord(np.array([10., 30.],
dtype=np.float32),
long_name='wibble', units='1')
- lon = iris.coords.DimCoord([-180, -90, 0, 90],
+ lon = DimCoord([-180, -90, 0, 90],
standard_name='longitude',
units='degrees', circular=True)
- lat = iris.coords.DimCoord([90, 0, -90],
+ lat = DimCoord([90, 0, -90],
standard_name='latitude', units='degrees')
cube.add_dim_coord(wibble_coord, [0])
cube.add_dim_coord(lat, [1])
@@ -294,39 +295,46 @@ def track_1d(duplicate_x=False):
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
"""
- cube = Cube(np.arange(11, dtype=np.int32), standard_name='air_temperature', units='K')
- bounds = np.column_stack([np.arange(11, dtype=np.int32), np.arange(11, dtype=np.int32) + 1])
+ cube = Cube(np.arange(11, dtype=np.int32), standard_name='air_temperature',
+ units='K')
+ bounds = np.column_stack([np.arange(11, dtype=np.int32),
+ np.arange(11, dtype=np.int32) + 1])
pts = bounds[:, 1]
- coord = iris.coords.AuxCoord(pts, 'projection_x_coordinate', units='1', bounds=bounds)
+ coord = AuxCoord(pts, 'projection_x_coordinate', units='1', bounds=bounds)
cube.add_aux_coord(coord, [0])
if duplicate_x:
- coord = iris.coords.AuxCoord(pts, 'projection_x_coordinate', units='1', bounds=bounds)
+ coord = AuxCoord(pts, 'projection_x_coordinate', units='1',
+ bounds=bounds)
cube.add_aux_coord(coord, [0])
- coord = iris.coords.AuxCoord(pts * 2, 'projection_y_coordinate', units='1', bounds=bounds * 2)
+ coord = AuxCoord(pts * 2, 'projection_y_coordinate', units='1',
+ bounds=bounds * 2)
cube.add_aux_coord(coord, 0)
return cube
def simple_2d_w_multidim_and_scalars():
data = np.arange(50, dtype=np.int32).reshape((5, 10))
- cube = iris.cube.Cube(data, long_name='test 2d dimensional cube', units='meters')
+ cube = iris.cube.Cube(data, long_name='test 2d dimensional cube',
+ units='meters')
# DimCoords
- dim1 = iris.coords.DimCoord(np.arange(5, dtype=np.float32) * 5.1 + 3.0, long_name='dim1', units='meters')
- dim2 = iris.coords.DimCoord(np.arange(10, dtype=np.int32), long_name='dim2', units='meters',
- bounds=np.arange(20, dtype=np.int32).reshape(10, 2))
+ dim1 = DimCoord(np.arange(5, dtype=np.float32) * 5.1 + 3.0,
+ long_name='dim1', units='meters')
+ dim2 = DimCoord(np.arange(10, dtype=np.int32),
+ long_name='dim2', units='meters',
+ bounds=np.arange(20, dtype=np.int32).reshape(10, 2))
# Scalars
- an_other = iris.coords.AuxCoord(3.0, long_name='an_other', units='meters')
- yet_an_other = iris.coords.DimCoord(23.3, standard_name='air_temperature',
- long_name='custom long name',
- var_name='custom_var_name',
- units='K')
+ an_other = AuxCoord(3.0, long_name='an_other', units='meters')
+ yet_an_other = DimCoord(23.3, standard_name='air_temperature',
+ long_name='custom long name',
+ var_name='custom_var_name', units='K')
# Multidim
- my_multi_dim_coord = iris.coords.AuxCoord(np.arange(50, dtype=np.int32).reshape(5, 10),
- long_name='my_multi_dim_coord', units='1',
- bounds=np.arange(200, dtype=np.int32).reshape(5, 10, 4))
+ my_multi_dim_coord = AuxCoord(np.arange(50, dtype=np.int32).reshape(5, 10),
+ long_name='my_multi_dim_coord', units='1',
+ bounds=np.arange(200, dtype=np.int32).
+ reshape(5, 10, 4))
cube.add_dim_coord(dim1, 0)
cube.add_dim_coord(dim2, 1)
@@ -361,18 +369,21 @@ def hybrid_height():
"""
data = np.arange(12, dtype='i8').reshape((3, 4))
- orography = icoords.AuxCoord([10, 25, 50, 5], standard_name='surface_altitude', units='m')
- model_level = icoords.AuxCoord([2, 1, 0], standard_name='model_level_number')
- level_height = icoords.DimCoord([100, 50, 10], long_name='level_height',
- units='m', attributes={'positive': 'up'},
- bounds=[[150, 75], [75, 20], [20, 0]])
- sigma = icoords.AuxCoord([0.8, 0.9, 0.95], long_name='sigma',
- bounds=[[0.7, 0.85], [0.85, 0.97], [0.97, 1.0]])
- hybrid_height = iris.aux_factory.HybridHeightFactory(level_height, sigma, orography)
+ orography = AuxCoord([10, 25, 50, 5], standard_name='surface_altitude',
+ units='m')
+ model_level = AuxCoord([2, 1, 0], standard_name='model_level_number')
+ level_height = DimCoord([100, 50, 10], long_name='level_height',
+ units='m', attributes={'positive': 'up'},
+ bounds=[[150, 75], [75, 20], [20, 0]])
+ sigma = AuxCoord([0.8, 0.9, 0.95], long_name='sigma',
+ bounds=[[0.7, 0.85], [0.85, 0.97], [0.97, 1.0]])
+ hybrid_height = iris.aux_factory.HybridHeightFactory(level_height,
+ sigma, orography)
cube = iris.cube.Cube(data, standard_name='air_temperature', units='K',
dim_coords_and_dims=[(level_height, 0)],
- aux_coords_and_dims=[(orography, 1), (model_level, 0), (sigma, 0)],
+ aux_coords_and_dims=[(orography, 1),
+ (model_level, 0), (sigma, 0)],
aux_factories=[hybrid_height])
return cube
@@ -381,25 +392,25 @@ def simple_4d_with_hybrid_height():
cube = iris.cube.Cube(np.arange(3*4*5*6, dtype='i8').reshape(3, 4, 5, 6),
"air_temperature", units="K")
- cube.add_dim_coord(iris.coords.DimCoord(np.arange(3, dtype='i8'), "time",
- units="hours since epoch"), 0)
- cube.add_dim_coord(iris.coords.DimCoord(np.arange(4, dtype='i8')+10,
- "model_level_number", units="1"), 1)
- cube.add_dim_coord(iris.coords.DimCoord(np.arange(5, dtype='i8')+20,
- "grid_latitude",
- units="degrees"), 2)
- cube.add_dim_coord(iris.coords.DimCoord(np.arange(6, dtype='i8')+30,
- "grid_longitude",
- units="degrees"), 3)
-
- cube.add_aux_coord(iris.coords.AuxCoord(np.arange(4, dtype='i8')+40,
- long_name="level_height",
- units="m"), 1)
- cube.add_aux_coord(iris.coords.AuxCoord(np.arange(4, dtype='i8')+50,
- long_name="sigma", units="1"), 1)
- cube.add_aux_coord(iris.coords.AuxCoord(np.arange(5*6, dtype='i8').reshape(5, 6)+100,
- long_name="surface_altitude",
- units="m"), [2, 3])
+ cube.add_dim_coord(DimCoord(np.arange(3, dtype='i8'), "time",
+ units="hours since epoch"), 0)
+ cube.add_dim_coord(DimCoord(np.arange(4, dtype='i8')+10,
+ "model_level_number", units="1"), 1)
+ cube.add_dim_coord(DimCoord(np.arange(5, dtype='i8')+20,
+ "grid_latitude",
+ units="degrees"), 2)
+ cube.add_dim_coord(DimCoord(np.arange(6, dtype='i8')+30,
+ "grid_longitude",
+ units="degrees"), 3)
+
+ cube.add_aux_coord(AuxCoord(np.arange(4, dtype='i8')+40,
+ long_name="level_height",
+ units="m"), 1)
+ cube.add_aux_coord(AuxCoord(np.arange(4, dtype='i8')+50,
+ long_name="sigma", units="1"), 1)
+ cube.add_aux_coord(AuxCoord(np.arange(5*6, dtype='i8').reshape(5, 6)+100,
+ long_name="surface_altitude",
+ units="m"), [2, 3])
cube.add_aux_factory(iris.aux_factory.HybridHeightFactory(
delta=cube.coord("level_height"),
@@ -449,7 +460,8 @@ def realistic_4d():
Returns a realistic 4d cube.
>>> print(repr(realistic_4d()))
-
+
"""
# the stock arrays were created in Iris 0.8 with:
@@ -477,7 +489,8 @@ def realistic_4d():
if not os.path.isfile(data_path):
raise IOError('Test data is not available at {}.'.format(data_path))
r = np.load(data_path)
- # sort the arrays based on the order they were originally given. The names given are of the form 'arr_1' or 'arr_10'
+ # sort the arrays based on the order they were originally given.
+ # The names given are of the form 'arr_1' or 'arr_10'
_, arrays = zip(*sorted(six.iteritems(r), key=lambda item: int(item[0][4:])))
lat_pts, lat_bnds, lon_pts, lon_bnds, level_height_pts, \
@@ -486,25 +499,37 @@ def realistic_4d():
ll_cs = RotatedGeogCS(37.5, 177.5, ellipsoid=GeogCS(6371229.0))
- lat = icoords.DimCoord(lat_pts, standard_name='grid_latitude', units='degrees',
- bounds=lat_bnds, coord_system=ll_cs)
- lon = icoords.DimCoord(lon_pts, standard_name='grid_longitude', units='degrees',
- bounds=lon_bnds, coord_system=ll_cs)
+ lat = icoords.DimCoord(lat_pts, standard_name='grid_latitude',
+ units='degrees', bounds=lat_bnds,
+ coord_system=ll_cs)
+ lon = icoords.DimCoord(lon_pts, standard_name='grid_longitude',
+ units='degrees', bounds=lon_bnds,
+ coord_system=ll_cs)
level_height = icoords.DimCoord(level_height_pts, long_name='level_height',
units='m', bounds=level_height_bnds,
attributes={'positive': 'up'})
- model_level = icoords.DimCoord(model_level_pts, standard_name='model_level_number',
+ model_level = icoords.DimCoord(model_level_pts,
+ standard_name='model_level_number',
units='1', attributes={'positive': 'up'})
- sigma = icoords.AuxCoord(sigma_pts, long_name='sigma', units='1', bounds=sigma_bnds)
- orography = icoords.AuxCoord(orography, standard_name='surface_altitude', units='m')
- time = icoords.DimCoord(time_pts, standard_name='time', units='hours since 1970-01-01 00:00:00')
- forecast_period = icoords.DimCoord(forecast_period_pts, standard_name='forecast_period', units='hours')
+ sigma = icoords.AuxCoord(sigma_pts, long_name='sigma', units='1',
+ bounds=sigma_bnds)
+ orography = icoords.AuxCoord(orography, standard_name='surface_altitude',
+ units='m')
+ time = icoords.DimCoord(time_pts, standard_name='time',
+ units='hours since 1970-01-01 00:00:00')
+ forecast_period = icoords.DimCoord(forecast_period_pts,
+ standard_name='forecast_period',
+ units='hours')
- hybrid_height = iris.aux_factory.HybridHeightFactory(level_height, sigma, orography)
+ hybrid_height = iris.aux_factory.HybridHeightFactory(level_height,
+ sigma, orography)
- cube = iris.cube.Cube(data, standard_name='air_potential_temperature', units='K',
- dim_coords_and_dims=[(time, 0), (model_level, 1), (lat, 2), (lon, 3)],
- aux_coords_and_dims=[(orography, (2, 3)), (level_height, 1), (sigma, 1),
+ cube = iris.cube.Cube(data, standard_name='air_potential_temperature',
+ units='K',
+ dim_coords_and_dims=[(time, 0), (model_level, 1),
+ (lat, 2), (lon, 3)],
+ aux_coords_and_dims=[(orography, (2, 3)),
+ (level_height, 1), (sigma, 1),
(forecast_period, None)],
attributes={'source': 'Iris test case'},
aux_factories=[hybrid_height])
@@ -516,7 +541,8 @@ def realistic_4d_no_derived():
Returns a realistic 4d cube without hybrid height
>>> print(repr(realistic_4d()))
-
+
"""
cube = realistic_4d()
@@ -534,24 +560,30 @@ def realistic_4d_w_missing_data():
data_archive = np.load(data_path)
data = ma.masked_array(data_archive['arr_0'], mask=data_archive['arr_1'])
- # sort the arrays based on the order they were originally given. The names given are of the form 'arr_1' or 'arr_10'
+ # sort the arrays based on the order they were originally given.
+ # The names given are of the form 'arr_1' or 'arr_10'
ll_cs = GeogCS(6371229)
- lat = iris.coords.DimCoord(np.arange(20, dtype=np.float32), standard_name='grid_latitude',
- units='degrees', coord_system=ll_cs)
- lon = iris.coords.DimCoord(np.arange(20, dtype=np.float32), standard_name='grid_longitude',
- units='degrees', coord_system=ll_cs)
- time = iris.coords.DimCoord([1000., 1003., 1006.], standard_name='time',
- units='hours since 1970-01-01 00:00:00')
- forecast_period = iris.coords.DimCoord([0.0, 3.0, 6.0], standard_name='forecast_period', units='hours')
- pressure = iris.coords.DimCoord(np.array([ 800., 900., 1000.], dtype=np.float32),
- long_name='pressure', units='hPa')
+ lat = DimCoord(np.arange(20, dtype=np.float32),
+ standard_name='grid_latitude',
+ units='degrees', coord_system=ll_cs)
+ lon = DimCoord(np.arange(20, dtype=np.float32),
+ standard_name='grid_longitude',
+ units='degrees', coord_system=ll_cs)
+ time = DimCoord([1000., 1003., 1006.], standard_name='time',
+ units='hours since 1970-01-01 00:00:00')
+ forecast_period = DimCoord([0.0, 3.0, 6.0],
+ standard_name='forecast_period',
+ units='hours')
+ pressure = DimCoord(np.array([800., 900., 1000.], dtype=np.float32),
+ long_name='pressure', units='hPa')
cube = iris.cube.Cube(data, long_name='missing data test data', units='K',
- dim_coords_and_dims=[(time, 0), (pressure, 1), (lat, 2), (lon, 3)],
+ dim_coords_and_dims=[(time, 0), (pressure, 1),
+ (lat, 2), (lon, 3)],
aux_coords_and_dims=[(forecast_period, 0)],
- attributes={'source':'Iris test case'})
+ attributes={'source': 'Iris test case'})
return 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..395a50d8ea
--- /dev/null
+++ b/lib/iris/tests/stock/_stock_2d_latlons.py
@@ -0,0 +1,350 @@
+# (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
+
+ # Start by making a "normal" cube with separate 1-D X and Y coords.
+ if regional:
+ # Make a small regional cube.
+ cube = sample_cube(xargs=(150., 243.75, 6), yargs=(-10., 40., 5))
+
+ # Add contiguous bounds.
+ for ax in ('x', 'y'):
+ cube.coord(axis=ax).guess_bounds()
+ else:
+ # Global data, but at a drastically reduced resolution.
+ cube = sample_cube(xargs=(37.5, 318.75, 6), yargs=(-85., 65., 5))
+
+ # Make contiguous bounds and adjust outer edges to ensure it is 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
+
+ # Now convert the 1-d coords to 2-d equivalents.
+ # Get original 1-d coords.
+ co_1d_x, co_1d_y = [cube.coord(axis=ax).copy() for ax in ('x', 'y')]
+
+ # Calculate 2-d equivalents.
+ 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:
+ # Put the cube locations into 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
new file mode 100644
index 0000000000..4a1344e83d
--- /dev/null
+++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
@@ -0,0 +1,327 @@
+# (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
+
+from cf_units import Unit
+from iris.cube import Cube
+from iris.coords import AuxCoord
+
+from iris.analysis.cartography import gridcell_angles
+from iris.tests.stock import sample_2d_latlons, lat_lon_cube
+
+
+def _2d_multicells_testcube(cellsize_degrees=1.0):
+ """
+ Create a test cube with a grid of X and Y points, where each gridcell
+ is independent (disjoint), arranged at an angle == the x-coord point.
+
+ """
+ # Setup np.linspace arguments to make the coordinate points.
+ x0, x1, nx = -164, 164, 9
+ y0, y1, ny = -75, 75, 7
+
+ lats = np.linspace(y0, y1, ny, endpoint=True)
+ lons_angles = np.linspace(x0, x1, nx, endpoint=True)
+ x_pts_2d, y_pts_2d = np.meshgrid(lons_angles, lats)
+
+ # Make gridcells rectangles surrounding these centrepoints, but also
+ # tilted at various angles (= same as x-point lons, as that's easy).
+
+ # Calculate centrepoint lons+lats : in radians, and shape (ny, nx, 1).
+ xangs, yangs = np.deg2rad(x_pts_2d), np.deg2rad(y_pts_2d)
+ xangs, yangs = [arr[..., None] for arr in (xangs, yangs)]
+ # Program which corners are up+down on each gridcell axis.
+ dx_corners = [[[-1, 1, 1, -1]]]
+ dy_corners = [[[-1, -1, 1, 1]]]
+ # Calculate the relative offsets in x+y at the 4 corners.
+ x_ofs_2d = cellsize_degrees * np.cos(xangs) * dx_corners
+ x_ofs_2d -= cellsize_degrees * np.sin(xangs) * dy_corners
+ y_ofs_2d = cellsize_degrees * np.cos(xangs) * dy_corners
+ y_ofs_2d += cellsize_degrees * np.sin(xangs) * dx_corners
+ # Apply a latitude stretch to make correct angles on the globe.
+ y_ofs_2d *= np.cos(yangs)
+ # Make bounds arrays by adding the corner offsets to the centrepoints.
+ x_bds_2d = x_pts_2d[..., None] + x_ofs_2d
+ y_bds_2d = y_pts_2d[..., None] + y_ofs_2d
+
+ # Create a cube with these points + bounds in its 'X' and 'Y' coords.
+ co_x = AuxCoord(points=x_pts_2d, bounds=x_bds_2d,
+ standard_name='longitude', units='degrees')
+ co_y = AuxCoord(points=y_pts_2d, bounds=y_bds_2d,
+ standard_name='latitude', units='degrees')
+ cube = Cube(np.zeros((ny, nx)))
+ cube.add_aux_coord(co_x, (0, 1))
+ cube.add_aux_coord(co_y, (0, 1))
+ return cube
+
+
+class TestGridcellAngles(tests.IrisTest):
+ def setUp(self):
+ # Make a small "normal" contiguous-bounded cube to test on.
+ # This one is regional.
+ self.standard_regional_cube = sample_2d_latlons(
+ regional=True, transformed=True)
+ # Record the standard correct angle answers.
+ result_cube = gridcell_angles(self.standard_regional_cube)
+ result_cube.convert_units('degrees')
+ self.standard_result_cube = result_cube
+ self.standard_small_cube_results = result_cube.data
+
+ def _check_multiple_orientations_and_latitudes(self,
+ method='mid-lhs, mid-rhs',
+ atol_degrees=0.005,
+ cellsize_degrees=1.0):
+
+ cube = _2d_multicells_testcube(cellsize_degrees=cellsize_degrees)
+
+ # Calculate gridcell angles at each point.
+ angles_cube = gridcell_angles(cube, cell_angle_boundpoints=method)
+
+ # Check that the results are a close match to the original intended
+ # gridcell orientation angles.
+ # NOTE: neither the above gridcell construction nor the calculation
+ # itself are exact : Errors scale as the square of gridcell sizes.
+ angles_cube.convert_units('degrees')
+ angles_calculated = angles_cube.data
+
+ # Note: the gridcell angles **should** just match the longitudes at
+ # each point
+ angles_expected = cube.coord('longitude').points
+
+ # Wrap both into standard range for comparison.
+ angles_calculated = (angles_calculated + 360.) % 360.
+ angles_expected = (angles_expected + 360.) % 360.
+
+ # Assert (toleranced) equality, and return results.
+ self.assertArrayAllClose(angles_calculated, angles_expected,
+ atol=atol_degrees)
+
+ return angles_calculated, angles_expected
+
+ def test_various_orientations_and_locations(self):
+ self._check_multiple_orientations_and_latitudes()
+
+ def test_result_form(self):
+ # Check properties of the result cube *other than* the data values.
+ test_cube = self.standard_regional_cube
+ result_cube = self.standard_result_cube
+ self.assertEqual(result_cube.long_name,
+ 'gridcell_angle_from_true_east')
+ self.assertEqual(result_cube.units, Unit('degrees'))
+ self.assertEqual(len(result_cube.coords()), 2)
+ self.assertEqual(result_cube.coord(axis='x'),
+ test_cube.coord(axis='x'))
+ self.assertEqual(result_cube.coord(axis='y'),
+ test_cube.coord(axis='y'))
+
+ def test_bottom_edge_method(self):
+ # Get results with the "other" calculation method + check to tolerance.
+ # A smallish cellsize should yield similar results in both cases.
+ r1, _ = self._check_multiple_orientations_and_latitudes()
+ r2, _ = self._check_multiple_orientations_and_latitudes(
+ method='lower-left, lower-right',
+ cellsize_degrees=0.1, atol_degrees=0.1)
+
+ # Not *exactly* the same : this checks we tested the 'other' method !
+ self.assertFalse(np.allclose(r1, r2))
+ # Note: results are a bit different in places. This is acceptable.
+ self.assertArrayAllClose(r1, r2, atol=0.1)
+
+ def test_bounded_coord_args(self):
+ # Check that passing the coords gives the same result as the cube.
+ co_x, co_y = (self.standard_regional_cube.coord(axis=ax)
+ for ax in ('x', 'y'))
+ result = gridcell_angles(co_x, co_y)
+ self.assertArrayAllClose(result.data,
+ self.standard_small_cube_results)
+
+ def test_coords_radians_args(self):
+ # Check it still works with coords converted to radians.
+ co_x, co_y = (self.standard_regional_cube.coord(axis=ax)
+ for ax in ('x', 'y'))
+ for coord in (co_x, co_y):
+ coord.convert_units('radians')
+ result = gridcell_angles(co_x, co_y)
+ self.assertArrayAllClose(result.data,
+ self.standard_small_cube_results)
+
+ def test_bounds_array_args(self):
+ # Check we can calculate from bounds values alone.
+ co_x, co_y = (self.standard_regional_cube.coord(axis=ax)
+ for ax in ('x', 'y'))
+ # Results drawn from coord bounds should be nearly the same,
+ # but not exactly, because of the different 'midpoint' values.
+ result = gridcell_angles(co_x.bounds, co_y.bounds)
+ self.assertArrayAllClose(result.data,
+ self.standard_small_cube_results, atol=0.1)
+
+ def test_unbounded_regional_coord_args(self):
+ # Remove the coord bounds to check points-based calculation.
+ co_x, co_y = (self.standard_regional_cube.coord(axis=ax)
+ for ax in ('x', 'y'))
+ for coord in (co_x, co_y):
+ coord.bounds = None
+ result = gridcell_angles(co_x, co_y)
+ # Note: in this case, we can expect the leftmost and rightmost columns
+ # to be rubbish, because the data is not global.
+ # But the rest should match okay.
+ self.assertArrayAllClose(result.data[:, 1:-1],
+ self.standard_small_cube_results[:, 1:-1])
+
+ def test_points_array_args(self):
+ # Check we can calculate from points arrays alone (no coords).
+ co_x, co_y = (self.standard_regional_cube.coord(axis=ax)
+ for ax in ('x', 'y'))
+ # As previous, the leftmost and rightmost columns are not good.
+ result = gridcell_angles(co_x.points, co_y.points)
+ self.assertArrayAllClose(result.data[:, 1:-1],
+ self.standard_small_cube_results[:, 1:-1])
+
+ def test_unbounded_global(self):
+ # For a contiguous global grid, a result based on points, i.e. with the
+ # bounds removed, should be a reasonable match for the 'ideal' one
+ # based on the bounds.
+
+ # Make a global cube + calculate ideal bounds-based results.
+ global_cube = sample_2d_latlons(transformed=True)
+ result_cube = gridcell_angles(global_cube)
+ result_cube.convert_units('degrees')
+ global_cube_results = result_cube.data
+
+ # Check a points-based calculation on the same basic grid.
+ co_x, co_y = (global_cube.coord(axis=ax)
+ for ax in ('x', 'y'))
+ for coord in (co_x, co_y):
+ coord.bounds = None
+ result = gridcell_angles(co_x, co_y)
+ # In this case, the match is actually rather poor (!).
+ self.assertArrayAllClose(result.data,
+ global_cube_results,
+ atol=7.5)
+ # Leaving off first + last columns again gives a decent result.
+ self.assertArrayAllClose(result.data[:, 1:-1],
+ global_cube_results[:, 1:-1])
+
+ # NOTE: although this looks just as bad as 'test_points_array_args',
+ # maximum errors there in the end columns are actually > 100 degrees !
+
+ def test_nonlatlon_coord_system(self):
+ # Check with points specified in an unexpected coord system.
+ cube = sample_2d_latlons(regional=True, rotated=True)
+ result = gridcell_angles(cube)
+ self.assertArrayAllClose(result.data,
+ self.standard_small_cube_results)
+ # Check that the result has transformed (true-latlon) coordinates.
+ self.assertEqual(len(result.coords()), 2)
+ x_coord = result.coord(axis='x')
+ y_coord = result.coord(axis='y')
+ self.assertEqual(x_coord.shape, cube.shape)
+ self.assertEqual(y_coord.shape, cube.shape)
+ self.assertIsNotNone(cube.coord_system)
+ self.assertIsNone(x_coord.coord_system)
+ self.assertIsNone(y_coord.coord_system)
+
+ def test_fail_coords_bad_units(self):
+ # Check error with bad coords units.
+ co_x, co_y = (self.standard_regional_cube.coord(axis=ax)
+ for ax in ('x', 'y'))
+ co_y.units = 'm'
+ with self.assertRaisesRegexp(ValueError, 'must have angular units'):
+ gridcell_angles(co_x, co_y)
+
+ def test_fail_nonarraylike(self):
+ # Check error with bad args.
+ co_x, co_y = 1, 2
+ with self.assertRaisesRegexp(ValueError,
+ 'must have array shape property'):
+ gridcell_angles(co_x, co_y)
+
+ def test_fail_non2d_coords(self):
+ # Check error with bad args.
+ cube = lat_lon_cube()
+ with self.assertRaisesRegexp(ValueError,
+ 'inputs must have 2-dimensional shape'):
+ gridcell_angles(cube)
+
+ def test_fail_different_shapes(self):
+ # Check error with mismatched shapes.
+ co_x, co_y = (self.standard_regional_cube.coord(axis=ax)
+ for ax in ('x', 'y'))
+ co_y = co_y[1:]
+ with self.assertRaisesRegexp(ValueError, 'must have same shape'):
+ gridcell_angles(co_x, co_y)
+
+ def test_fail_different_coord_system(self):
+ # Check error with mismatched coord systems.
+ cube = sample_2d_latlons(regional=True, rotated=True)
+ cube.coord(axis='x').coord_system = None
+ with self.assertRaisesRegexp(ValueError,
+ 'must have same coordinate system'):
+ gridcell_angles(cube)
+
+ def test_fail_cube_dims(self):
+ # Check error with mismatched cube dims.
+ cube = self.standard_regional_cube
+ # Make 5x6 into 5x5.
+ cube = cube[:, :-1]
+ co_x = cube.coord(axis='x')
+ pts, bds = co_x.points, co_x.bounds
+ co_new_x = co_x.copy(points=pts.transpose((1, 0)),
+ bounds=bds.transpose((1, 0, 2)))
+ cube.remove_coord(co_x)
+ cube.add_aux_coord(co_new_x, (1, 0))
+ with self.assertRaisesRegexp(ValueError,
+ 'must have the same cube dimensions'):
+ gridcell_angles(cube)
+
+ def test_fail_coord_noncoord(self):
+ # Check that passing a coord + an array gives an error.
+ co_x, co_y = (self.standard_regional_cube.coord(axis=ax)
+ for ax in ('x', 'y'))
+ with self.assertRaisesRegexp(ValueError,
+ 'is a Coordinate, but .* is not'):
+ gridcell_angles(co_x, co_y.bounds)
+
+ def test_fail_noncoord_coord(self):
+ # Check that passing an array + a coord gives an error.
+ co_x, co_y = (self.standard_regional_cube.coord(axis=ax)
+ for ax in ('x', 'y'))
+ with self.assertRaisesRegexp(ValueError,
+ 'is a Coordinate, but .* is not'):
+ gridcell_angles(co_x.points, co_y)
+
+ def test_fail_bad_method(self):
+ with self.assertRaisesRegexp(ValueError,
+ 'unrecognised cell_angle_boundpoints'):
+ self._check_multiple_orientations_and_latitudes(
+ method='something_unknown')
+
+
+if __name__ == "__main__":
+ tests.main()
diff --git a/lib/iris/tests/unit/cube/test_Cube.py b/lib/iris/tests/unit/cube/test_Cube.py
index 9c96cf6b48..5633417c80 100644
--- a/lib/iris/tests/unit/cube/test_Cube.py
+++ b/lib/iris/tests/unit/cube/test_Cube.py
@@ -1430,7 +1430,9 @@ def test__masked_scalar_arraymask(self):
self._check_copy(cube, cube.copy())
def test__lazy(self):
- cube = Cube(as_lazy_data(np.array([1, 0])))
+ # Note: multiple chunks added as a workaround suggested to dask#3751,
+ # which is fixed in dask#3754.
+ cube = Cube(as_lazy_data(np.array([1, 0]), chunks=(1, 1)))
self._check_copy(cube, cube.copy())
diff --git a/lib/iris/tests/unit/plot/test_2d_coords.py b/lib/iris/tests/unit/plot/test_2d_coords.py
new file mode 100644
index 0000000000..7397921861
--- /dev/null
+++ b/lib/iris/tests/unit/plot/test_2d_coords.py
@@ -0,0 +1,167 @@
+# (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 handling and plotting of 2-dimensional coordinates"""
+
+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.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 sample_2d_latlons
+from iris.tests.stock import make_bounds_discontiguous_at_point
+
+
+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()
+ self.overspan_cube = cube3d_2dcoords()
+
+ # 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 = 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:
+ # single_dims = lat_lon_cube()
+ # lon = single_dims.coord('longitude')
+ # lat = single_dims.coord('latitude')
+ # big_lon, big_lat = testdata.grid_coords_2d_from_1d(lon, lat)
+ # mixed_dims = single_dims.copy()
+ # mixed_dims.remove_coord(lon)
+ # # TODO Fix this coord addition:
+ # # When adding an aux_coord, the function '_check_multidim_metadata'
+ # # throws an error as it requires coord.shape to be (1, ) instead of
+ # # (3, 4) or whatever.
+ # mixed_dims.add_aux_coord(big_lon)
+ #
+ # # mixed_dims is now a cube with 2 1D dim coords and an additional
+ # # 2D aux coord.
+ # self.mixed_dims = mixed_dims
+
+ self.mode = coords.BOUND_MODE
+
+ def test_2d_coords_identified(self):
+ # Test that 2d coords are identified in the plot definition without
+ # having to be specified in the args without any errors.
+ cube = self.multidim_cube
+ defn = iplt._get_plot_defn(cube, mode=self.mode)
+ self.assertEqual([coord.name() for coord in defn.coords],
+ ['bar', 'foo'])
+
+ def test_2d_coords_custom_picked(self):
+ # Test that 2d coords which are specified in the args will be
+ # accepted without any errors.
+ cube = self.multidim_cube
+ defn = iplt._get_plot_defn_custom_coords_picked(cube, ('foo', 'bar'),
+ self.mode)
+ self.assertEqual([coord.name() for coord in defn.coords],
+ ['bar', 'foo'])
+
+ def test_2d_coords_as_integers(self):
+ # Test that if you pass in 2d coords as args in the form of integers,
+ # they will still be correctly identified without any errors.
+ cube = self.multidim_cube
+ defn = iplt._get_plot_defn_custom_coords_picked(cube, (0, 1),
+ self.mode)
+ self.assertEqual([coord for coord in defn.coords],
+ [1, 0])
+
+ def test_total_span_check(self):
+ # Test that an error is raised if a user tries to plot a 2d coord
+ # against a different coord, making total number of dimensions 3.
+ cube = self.overspan_cube
+ with self.assertRaises(ValueError):
+ iplt._get_plot_defn_custom_coords_picked(cube, ('wibble', 'foo'),
+ self.mode)
+
+ # def test_2dcoord_with_1dcoord(self):
+ # # TODO Generate a cube with one 2d coord and one 1d coord
+ # # TODO Try and plot them against each other
+ # # 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)
+
+ def test_map_common_not_enough_bounds(self):
+ # Test that a lat-lon cube with 2d coords and 2 bounds per point
+ # throws an error in contiguity checks.
+ cube = self.multidim_cube
+ cube.coord('foo').rename('longitude')
+ cube.coord('bar').rename('latitude')
+ with self.assertRaises(ValueError):
+ plot_defn = iplt._get_plot_defn(cube, self.mode)
+ iplt._map_common('pcolor', None, self.mode, cube, plot_defn)
+
+ def test_map_common_2d(self):
+ # Test that a cube with 2d coords can be plotted as a map.
+ cube = self.latlon_2d
+ # Get necessary variables from _get_plot_defn to check that the test
+ # case will be accepted by _map_common.
+ plot_defn = iplt._get_plot_defn(cube, self.mode)
+ 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_unmasked(self):
+ # Check that an error occurs when the contiguity check finds
+ # discontiguous bounds but unmasked data.
+ cube = self.latlon_2d
+ cube.data.mask = ma.nomask
+ coord = cube.coord('longitude')
+ with self.assertRaises(ValueError):
+ iplt._check_contiguity_and_bounds(coord, cube.data)
+
+ def test_draw_2d_from_bounds(self):
+ # Test this function will not raise an error even with our most
+ # awkward but supported cube.
+ cube = self.latlon_2d
+ result = iplt._draw_2d_from_bounds('pcolormesh', cube)
+ self.assertTrue(result)
+
+
+if __name__ == '__main__':
+ tests.main()
diff --git a/requirements/core.txt b/requirements/core.txt
index c28cdb04a1..997c91c665 100644
--- a/requirements/core.txt
+++ b/requirements/core.txt
@@ -6,7 +6,7 @@
cartopy
cf_units>=2
cftime
-dask[array]>=0.17.1 #conda: dask>=0.17.1
+dask[array]==0.18.1 #conda: dask==0.18.1
matplotlib>=2
netcdf4
numpy>=1.14