Vector plots 2

lib/iris/analysis/_grid_angles.py

@@ -24,6 +24,7 @@
 
 import numpy as np
 
+import cartopy.crs as ccrs
 import iris
 
 
@@ -65,13 +66,13 @@ def _latlon_from_xyz(xyz):
     Returns:
 
         lonlat : (array)
-            spherical angles, of dims (2, <others>), in radians.
+            spherical angles, of dims (2, <others>), in degrees.
             Dim 0 maps longitude, latitude.
 
     """
-    lons = np.arctan2(xyz[1], xyz[0])
+    lons = np.rad2deg(np.arctan2(xyz[1], xyz[0]))
     axial_radii = np.sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1])
-    lats = np.arctan2(xyz[2], axial_radii)
+    lats = np.rad2deg(np.arctan2(xyz[2], axial_radii))
     return np.array([lons, lats])
 
 
@@ -87,64 +88,22 @@ def _angle(p, q, r):
     p, q, r, are all 2-element arrays [lon, lat] of angles in degrees.
 
     """
-#    old_style = True
-    old_style = False
-    if old_style:
-        mid_lons = np.deg2rad(q[0])
+    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)
+    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
+    index = pr_norm == 0
+    pr_norm[index] = 1
 
-        cosine = np.maximum(np.minimum(pr_top / pr_norm, 1), -1)
-        cosine[index] = 0
+    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
-    else:
-        # Calculate unit vectors.
-        midpt_lons, midpt_lats = q[0], q[1]
-        lmb_r, phi_r = (np.deg2rad(arr) for arr in (midpt_lons, midpt_lats))
-        phi_hatvec_x = -np.sin(phi_r) * np.cos(lmb_r)
-        phi_hatvec_y = -np.sin(phi_r) * np.sin(lmb_r)
-        phi_hatvec_z = np.cos(phi_r)
-        shape_xyz = (1,) + midpt_lons.shape
-        phi_hatvec = np.concatenate([arr.reshape(shape_xyz)
-                                     for arr in (phi_hatvec_x,
-                                                 phi_hatvec_y,
-                                                 phi_hatvec_z)])
-        lmb_hatvec_z = np.zeros(midpt_lons.shape)
-        lmb_hatvec_y = np.cos(lmb_r)
-        lmb_hatvec_x = -np.sin(lmb_r)
-        lmb_hatvec = np.concatenate([arr.reshape(shape_xyz)
-                                     for arr in (lmb_hatvec_x,
-                                                 lmb_hatvec_y,
-                                                 lmb_hatvec_z)])
-
-        pr = _3d_xyz_from_latlon(r[0], r[1]) - _3d_xyz_from_latlon(p[0], p[1])
-
-        # Dot products to form true-northward / true-eastward projections.
-        pr_cmpt_e = np.sum(pr * lmb_hatvec, axis=0)
-        pr_cmpt_n = np.sum(pr * phi_hatvec, axis=0)
-        psi = np.arctan2(pr_cmpt_n, pr_cmpt_e)
-
-        # TEMPORARY CHECKS:
-        # ensure that the two unit vectors are perpendicular.
-        dotprod = np.sum(phi_hatvec * lmb_hatvec, axis=0)
-        assert np.allclose(dotprod, 0.0)
-        # ensure that the vector components carry the original magnitude.
-        mag_orig = np.sum(pr * pr)
-        mag_rot = np.sum(pr_cmpt_e * pr_cmpt_e) + np.sum(pr_cmpt_n * pr_cmpt_n)
-        rtol = 1.e-3
-        check = np.allclose(mag_rot, mag_orig, rtol=rtol)
-        if not check:
-            print(mag_rot, mag_orig)
-            assert np.allclose(mag_rot, mag_orig, rtol=rtol)
-
-    return psi
+    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'):
@@ -201,7 +160,7 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'):
         angles : (2-dimensional cube)
 
             Cube of angles of grid-x vector from true Eastward direction for
-            each gridcell, in radians.
+            each gridcell, in degrees.
             It also has longitude and latitude coordinates.  If coordinates
             were input the output has identical ones :  If the input was 2d
             arrays, the output coords have no bounds; or, if the input was 3d
@@ -223,9 +182,17 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'):
 
     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()
-        x_coord.convert_units('degrees')
-        y_coord.convert_units('degrees')
+
+        # 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 {}.')
@@ -234,19 +201,36 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'):
             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))
-# NOTE: would like to check that dims are in correct order, but can't do that
-# if there is no cube.
-# TODO: **document** -- another input format requirement
-#        x_dims, y_dims = (cube.coord_dims(co) for co in (x_coord, y_coord))
-#        if x_dims != (0, 1) or y_dims != (0, 1):
-#            msg = ('Coordinate inputs must map to cube dimensions (0, 1). ',
-#                   'Got x-dims of {} and y-dims of {}.')
-#            raise ValueError(msg.format(x_dims, y_dims))
+        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()
+            # 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(x, y, src_crs=crs_src)
+            x = pts[..., 0].reshape(shape)
+            y = pts[..., 1].reshape(shape)
+
     elif hasattr(x, 'bounds') or hasattr(y, 'bounds'):
         # One was a Coord, and the other not ?
         is_and_not = ('x', 'y')
@@ -255,16 +239,22 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'):
         msg = 'Input {!r} is a Coordinate, but {!r} is not.'
         raise ValueError(*is_and_not)
 
-    # Now have either 2 points arrays or 2 bounds arrays.
-    # Construct (lhs, mid, rhs) where these represent 3 adjacent points with
-    # increasing longitudes.
+    # 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 X
+    # indices.
+    # Also make suitable X and Y coordinates for the result cube.
     if x.ndim == 2:
-        # PROBLEM: we can't use this if data is not full-longitudes,
-        # i.e. rhs of array must connect to lhs (aka 'circular' coordinate).
-        # But we have no means of checking that ?
-
+        # Data is points arrays.
         # Use previous + subsequent points along longitude-axis as references.
-        # NOTE: we also have no way to check that dim #2 really is the 'X' dim.
+
+        # 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)
@@ -275,9 +265,11 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'):
             x_coord = iris.coords.AuxCoord(y, standard_name='longitude',
                                            units='degrees')
     else:
-        # Get lhs and rhs locations by averaging top+bottom each side.
+        # Data is bounds arrays.
+        # Use different bounds points in the longitude-axis 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)
@@ -294,8 +286,8 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'):
                                         list(angle_boundpoints_vals.keys())))
         if not x_coord:
             # Create bounded coords for result cube.
-            # Use average lhs+rhs points in 3d to get 'mid' points, as coords
-            # with no points are not allowed.
+            # 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.
@@ -305,28 +297,30 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'):
             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 = np.rad2deg(_latlon_from_xyz(lhs_xyz))
-        rhs = np.rad2deg(_latlon_from_xyz(rhs_xyz))
-        # mid is coord.points, whether input or made up.
+        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='radians')
+                            units='degrees')
     result.add_aux_coord(x_coord, (0, 1))
     result.add_aux_coord(y_coord, (0, 1))
     return result
 
 
-def true_vectors_from_grid_vectors(u_cube, v_cube,
-                                   grid_angles_cube=None,
-                                   grid_angles_kwargs=None):
+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

lib/iris/analysis/cartography.py

@@ -37,10 +37,7 @@
 import iris.coord_systems
 import iris.exceptions
 from iris.util import _meshgrid
-from ._grid_angles import (
-    gridcell_angles,
-    true_vectors_from_grid_vectors as rotate_grid_vectors)
-
+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.

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=''):

lib/iris/tests/stock/_stock_2d_latlons.py

@@ -233,18 +233,19 @@ def sample_cube(xargs, yargs):
         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:
-        # Extract small region.
+        # Make a small regional cube.
         cube = sample_cube(xargs=(150., 243.75, 6), yargs=(-10., 40., 5))
 
-        # Make contiguous bounds.
+        # Add contiguous bounds.
         for ax in ('x', 'y'):
             cube.coord(axis=ax).guess_bounds()
     else:
-        # Global data, but drastically reduced resolution.
+        # Global data, but at a drastically reduced resolution.
         cube = sample_cube(xargs=(37.5, 318.75, 6), yargs=(-85., 65., 5))
 
-        # Patch bounds to ensure it is still contiguous + global.
+        # Make contiguous bounds and adjust outer edges to ensure it is global.
         for name in ('longitude', 'latitude'):
             coord = cube.coord(name)
             coord.guess_bounds()
@@ -258,8 +259,11 @@ def sample_cube(xargs, yargs):
                 bds[-1, 1] = 90.0
             coord.bounds = bds
 
-    # Get 1d coordinate points + bounds + calculate 2d equivalents.
+    # 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.
@@ -271,7 +275,7 @@ def sample_cube(xargs, yargs):
         cube.add_aux_coord(coord, (0, 1))
 
     if transformed or rotated:
-        # Take the lats + lons as being in a rotated coord system.
+        # 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.