From ebfc4c10e6d2fdf8dc98d67c46dcdeae4fc3e729 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Fri, 27 Jul 2018 19:11:21 +0100 Subject: [PATCH 01/15] Small improvements; first sensible tests. --- lib/iris/analysis/_grid_angles.py | 99 ++-- lib/iris/analysis/cartography.py | 7 +- .../cartography/test_gridcell_angles.py | 468 ++++++++++++++++++ 3 files changed, 506 insertions(+), 68 deletions(-) create mode 100644 lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py index ffe6d51c7f..f847e0d00f 100644 --- a/lib/iris/analysis/_grid_angles.py +++ b/lib/iris/analysis/_grid_angles.py @@ -87,62 +87,20 @@ 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) + psi = np.arccos(cosine) * np.sign(r[1] - p[1]) + psi[index] = np.nan return psi @@ -223,6 +181,7 @@ 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') @@ -255,16 +214,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. + # Now have either 2 points arrays (ny, nx) or 2 bounds arrays (ny, nx, 4). # Construct (lhs, mid, rhs) where these represent 3 adjacent points with # increasing longitudes. + # 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 +240,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 points 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 +261,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,10 +272,11 @@ 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. + # '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. @@ -321,12 +289,13 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): 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 diff --git a/lib/iris/analysis/cartography.py b/lib/iris/analysis/cartography.py index dadcdedc86..af55a9de57 100644 --- a/lib/iris/analysis/cartography.py +++ b/lib/iris/analysis/cartography.py @@ -37,10 +37,11 @@ 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 +# Avoid pycodestyle warnings for unused imports. +gridcell_angles = gridcell_angles +rotate_grid_vectors = rotate_grid_vectors # List of contents to control Sphinx autodocs. # Unfortunately essential to get docs for the grid_angles functions. 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..26de728dad --- /dev/null +++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py @@ -0,0 +1,468 @@ +# (C) British Crown Copyright 2018, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +""" +Unit tests for the function +:func:`iris.analysis.cartography.gridcell_angles`. + +""" +from __future__ import (absolute_import, division, print_function) +from six.moves import (filter, input, map, range, zip) # noqa + +# Import iris.tests first so that some things can be initialised before +# importing anything else. +import iris.tests as tests + +import numpy as np +import numpy.ma as ma + +import cartopy.crs as ccrs +from iris.cube import Cube +from iris.coords import DimCoord, AuxCoord +import iris.coord_systems +from iris.analysis.cartography import unrotate_pole + +from iris.analysis.cartography import (gridcell_angles, + rotate_grid_vectors) + +import matplotlib.pyplot as plt +from orca_utils.plot_testing.blockplot_from_bounds import blockplot_2dll + + +def _rotated_grid_sample(pole_lat=15, pole_lon=-180, + lon_bounds=np.linspace(-30, 30, 6, endpoint=True), + lat_bounds=np.linspace(-30, 30, 6, endpoint=True)): + # Calculate *true* lat_bounds+lon_bounds for the rotated grid. + lon_bounds = np.array(lon_bounds, dtype=float) + lat_bounds = np.array(lat_bounds, dtype=float) + # Construct centrepoints. + lons = 0.5 * (lon_bounds[:-1] + lon_bounds[1:]) + lats = 0.5 * (lat_bounds[:-1] + lat_bounds[1:]) + # Convert all to full 2d arrays. + lon_bounds, lat_bounds = np.meshgrid(lon_bounds, lat_bounds) + lons, lats = np.meshgrid(lons, lats) + # Calculate true lats+lons for all points. + lons_true_bds, lats_true_bds = unrotate_pole(lon_bounds, lat_bounds, + pole_lon, pole_lat) + lons_true, lats_true = unrotate_pole(lons, lats, pole_lon, pole_lat) + # Make the 'unified' bounds into contiguous (ny, nx, 4) arrays. + def expand_unified_bds(bds): + ny, nx = bds.shape + bds_4 = np.zeros((ny - 1, nx - 1, 4)) + bds_4[:, :, 0] = bds[:-1, :-1] + bds_4[:, :, 1] = bds[:-1, 1:] + bds_4[:, :, 2] = bds[1:, 1:] + bds_4[:, :, 3] = bds[1:, :-1] + return bds_4 + + lon_true_bds4, lat_true_bds4 = (expand_unified_bds(bds) + for bds in (lons_true_bds, lats_true_bds)) + # Make these into a 2d-latlon grid for a cube + cube = Cube(np.zeros(lon_true_bds4.shape[:-1])) + co_x = AuxCoord(lons_true, bounds=lon_true_bds4, + standard_name='longitude', units='degrees') + co_y = AuxCoord(lats_true, bounds=lat_true_bds4, + standard_name='latitude', units='degrees') + cube.add_aux_coord(co_x, (0, 1)) + cube.add_aux_coord(co_y, (0, 1)) + return cube + + +class TestGridcellAngles(tests.IrisTest): + def _singlecell_30deg_cube(self, x0=90., y0=0., dx=20., dy=10.): + x_pts = np.array([[x0]]) + y_pts = np.array([[y0]]) + x_bds = x0 + dx * np.array([[[-1., 1, 0.5, -1.5]]]) +# self.assertArrayAllClose(x_bds, np.array([[[70., 110, 100, 60]]])) + y_bds = y0 + dy * np.array([[[-1., 1, 3, 1]]]) +# self.assertArrayAllClose(y_bds, np.array([[[-10., 10, 30, 10]]])) + co_x = AuxCoord(points=x_pts, bounds=x_bds, + standard_name='longitude', units='degrees') + co_y = AuxCoord(points=y_pts, bounds=y_bds, + standard_name='latitude', units='degrees') + cube = Cube(np.zeros((1, 1))) + cube.add_aux_coord(co_x, (0, 1)) + cube.add_aux_coord(co_y, (0, 1)) + return cube + + def _singlecell_diamond_cube(self, x0=90., y0=0., dy=10., dx_eq=None): + if dx_eq is None: + dx_eq = dy + x_pts = np.array([[x0]]) + y_pts = np.array([[y0]]) + dx = dx_eq / np.cos(np.deg2rad(y0)) + x_bds = np.array([[[x0, x0 + dx, x0, x0 - dx]]]) + y_bds = np.array([[[y0 - dy, y0, y0 + dy, y0]]]) + co_x = AuxCoord(points=x_pts, bounds=x_bds, + standard_name='longitude', units='degrees') + co_y = AuxCoord(points=y_pts, bounds=y_bds, + standard_name='latitude', units='degrees') + cube = Cube(np.zeros((1, 1))) + cube.add_aux_coord(co_x, (0, 1)) + cube.add_aux_coord(co_y, (0, 1)) + return cube + + def _check_different_orientations_and_latitudes( + self, + method='mid-lhs, mid-rhs', + atol_degrees=0.005): + ny, nx = 7, 9 + x0, x1 = -164, 164 + y0, y1 = -75, 75 + lats = np.linspace(y0, y1, ny, endpoint=True) + angles = np.linspace(x0, x1, nx, endpoint=True) + x_pts_2d, y_pts_2d = np.meshgrid(angles, lats) + + # Make gridcells rectangles surrounding these centrepoints, but also + # tilted at various angles (= same as x-point lons, as that's easy). + dx = 1.0 # half-width of gridcells, in degrees + dy = dx # half-height of gridcells, in degrees + # 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 = dx * np.cos(xangs) * dx_corners + x_ofs_2d -= dy * np.sin(xangs) * dy_corners + y_ofs_2d = dy * np.cos(xangs) * dy_corners + y_ofs_2d += dx * 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)) + + # 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: expand the original 1-d test angles into the full result shape, + # just to please 'np.testing.assert_allclose', which doesn't broadcast. + angles_expected = np.zeros(angles_cube.shape) + angles_expected[:] = angles + self.assertArrayAllClose(angles_calculated, angles_expected, + atol=atol_degrees) + return angles_calculated + + def test_different_orientations_and_latitudes(self): + self._check_different_orientations_and_latitudes() + + def test_different_methods(self): + r1 = self._check_different_orientations_and_latitudes( + 'mid-lhs, mid-rhs', atol_degrees=5.0) + r2 = self._check_different_orientations_and_latitudes( + 'lower-left, lower-right', atol_degrees=5.0) +# print(np.round(r1 - r2, 3)) +# for fn in (np.max, np.mean): +# print(fn.__name__, fn(np.abs(r1 - r2))) + atol = 1.0 # !!!!!!! + self.assertArrayAllClose(r1, r2, atol=atol) + + +# def test_single_cell_equatorial(self): +# plt.switch_backend('tkagg') +# plt.figure(figsize=(10,10)) +## ax = plt.axes(projection=ccrs.Mercator()) +## ax = plt.axes(projection=ccrs.NorthPolarStereo()) +# ax = plt.axes(projection=ccrs.Orthographic(central_longitude=90., +# central_latitude=30.)) +# +# lon0 = 90.0 +# dy = 1.0 +# dx = 3.0 +# y_0, y_n, ny = -80, 80, 9 +# angles = [] +# for lat in np.linspace(y_0, y_n, ny): +# cube = self._singlecell_diamond_cube(x0=lon0, y0=lat, +# dy=dy, dx_eq=dx) +# angles_cube = gridcell_angles(cube, +## cell_angle_boundpoints='mid-lhs, mid-rhs') +# cell_angle_boundpoints='lower-left, lower-right') +# tmp_cube = angles_cube.copy() +# tmp_cube.convert_units('degrees') +## print('') +## print(lat) +## co_x, co_y = (cube.coord(axis=ax) for ax in ('x', 'y')) +## print() +## print(' at : {}, {}'.format(co_x.points[0, 0], co_y.points[0, 0])) +## print(' x-bds:') +## print(co_x.bounds) +## print(' y-bds:') +## print(co_y.bounds) +# angle = tmp_cube.data[0, 0] +# angles.append(angle) +# print(lat, angle) +# blockplot_2dll(cube) +# +# ax.coastlines() +# ax.set_global() +# +# # Plot constant NEly (45deg) arrows. +# xx = np.array([lon0] * ny) +# yy = np.linspace(y_0, y_n, ny) - dy +# uu = np.array([1.0] * ny) +# plt.quiver(xx, yy, +# uu, np.cos(np.deg2rad(yy)), +# zorder=2, color='red', +## scale_units='xy', +# angles='uv', +# transform=ccrs.PlateCarree()) +# +# # Also plot returned angles. +# angles_arr_rad = np.deg2rad(angles) +# u_arr = uu * np.cos(angles_arr_rad) +# v_arr = uu * np.sin(angles_arr_rad) * np.cos(np.deg2rad(yy)) +# +# plt.quiver(xx, yy, +# u_arr, +# v_arr, +# zorder=2, color='magenta', +# scale_units='xy', +# width=0.005, +# scale=0.2e-6, +## width=0.5, +# transform=ccrs.PlateCarree()) +# +# plt.show() + + +# def test_values(self): +# # Construct a rotated-pole grid and check angle calculation. +# testcube = _rotated_grid_sample() +# +# cell_angle_boundpoints = 'mid-lhs, mid-rhs' +## cell_angle_boundpoints = 'lower-left, lower-right' +## cell_angle_boundpoints = 'garble' +# angles_cube = gridcell_angles( +# testcube, +# cell_angle_boundpoints=cell_angle_boundpoints) +# angles_cube.convert_units('radians') +# +# # testing phase... +# print(np.rad2deg(angles_cube.data)) +# +# import matplotlib.pyplot as plt +# plt.switch_backend('tkagg') +# +## plot_map = 'north_polar_stereographic' +## plot_map = 'plate_carree' +## plot_map = 'mercator' +# plot_map = 'north_polar_orthographic' +# if plot_map == 'plate_carree': +# scale = 0.1 +# map_proj = ccrs.PlateCarree() +# elif plot_map == 'mercator': +# scale = 3.0e-6 +# map_proj = ccrs.Mercator() +# map_proj._threshold *= 0.01 +# elif plot_map == 'north_polar_orthographic': +# scale = 3.0e-6 +# map_proj = ccrs.Orthographic(central_longitude=0.0, +# central_latitude=90.0,) +# map_proj._threshold *= 0.01 +# elif plot_map == 'north_polar_stereographic': +# scale = 3.0e-6 +# map_proj = ccrs.NorthPolarStereo() +# else: +# assert 0 +# +# ax = plt.axes(projection=map_proj) +# data_proj = ccrs.PlateCarree() +# +# deg_scale = 10.0 +# +## angles = 'uv' +# angles = 'xy' +# +# ax.coastlines() +# ax.gridlines() +# for i_bnd in range(4): +# color = ['black', 'red', 'blue', 'magenta'][i_bnd] +# plt.plot(testcube.coord('longitude').bounds[..., i_bnd], +# testcube.coord('latitude').bounds[..., i_bnd], +# '+', markersize=10., markeredgewidth=2., +# markerfacecolor=color, markeredgecolor=color, +# transform=data_proj) +# +# +# # Show plain 0,1 + 1,0 (PlateCarree) vectors unrotated at the given points. +# pts_shape = testcube.coord('longitude').shape +# ny, nx = pts_shape +# u0 = np.ones(pts_shape) +# v0 = np.zeros(pts_shape) +# u1 = v0.copy() +# v1 = u0.copy() +# +# x0s = testcube.coord('longitude').points +# y0s = testcube.coord('latitude').points +# yscale = np.cos(np.deg2rad(y0s)) +# plt.quiver(x0s, y0s, u0, v0 * yscale, +# color='blue', width=0.005, +# headwidth=2., # headlength=1.0, headaxislength=0.7, +# angles=angles, +# scale_units='xy', scale=scale, +# transform=data_proj) +# plt.quiver(x0s, y0s, u1, v1 * yscale, +# color='red', width=0.005, +# headwidth=2., # headlength=1.0, headaxislength=0.7, +# angles=angles, +# scale_units='xy', scale=scale, +# transform=data_proj) +# +# # Add 45deg arrows (NEly), still on a PlateCarree map. +# plt.quiver(x0s, y0s, v1, v1 * yscale, +# color='green', width=0.005, +# headwidth=2., # headlength=1.0, headaxislength=0.7, +# angles=angles, +# scale_units='xy', scale=scale, +# transform=data_proj) +# +# +# +# # +# # Repeat the above plotting short lines INSTEAD of quiver. +# # +# u0d = x0s + deg_scale * u0 +# v0d = y0s + deg_scale * v0 +# u1d = x0s + deg_scale * u1 +# v1d = y0s + deg_scale * v1 +# u2d = x0s + deg_scale * u0 +# v2d = y0s + deg_scale * v1 +# for iy in range(ny): +# for ix in range(nx): +# plt.plot([x0s[iy, ix], u0d[iy, ix]], +# [y0s[iy, ix], v0d[iy, ix]], +# ':', color='blue', linewidth=0.5, +# transform=data_proj) +# plt.plot([x0s[iy, ix], u1d[iy, ix]], +# [y0s[iy, ix], v1d[iy, ix]], +# ':', color='red', linewidth=0.5, +# transform=data_proj) +# plt.plot([x0s[iy, ix], u2d[iy, ix]], +# [y0s[iy, ix], v2d[iy, ix]], +# ':', color='green', linewidth=0.5, +# transform=data_proj) +# +# +# # Overplot BL-BR and BL-TL lines from the cell bounds. +# co_lon, co_lat = [testcube.coord(name).copy() +# for name in ('longitude', 'latitude')] +# for co in (co_lon, co_lat): +# co.convert_units('degrees') +# lon_bds, lat_bds = [co.bounds for co in (co_lon, co_lat)] +## ny, nx = lon_bds.shape[:-1] +# for iy in range(ny): +# for ix in range(nx): +# x0, y0 = lon_bds[iy, ix, 0], lat_bds[iy, ix, 0] +# x1, y1 = lon_bds[iy, ix, 1], lat_bds[iy, ix, 1] +# x2, y2 = lon_bds[iy, ix, 3], lat_bds[iy, ix, 3] +# plt.plot([x0, x1], [y0, y1], 'x-', +# color='orange', +# transform=data_proj) +# plt.plot([x0, x2], [y0, y2], 'x-', +# color='orange', linestyle='--', +# transform=data_proj) +# +# # Plot U0, rotated by cell angles, also at cell bottom-lefts. +# u0_cube, u1_cube, v0_cube, v1_cube = [testcube.copy(data=aa) +# for aa in (u0, v0, u1, v1)] +# u0r_cube, v0r_cube = rotate_grid_vectors( +# u0_cube, v0_cube, grid_angles_cube=angles_cube) +# u0r, v0r = [cube.data for cube in (u0r_cube, v0r_cube)] +# +# xbl, ybl = lon_bds[..., 0], lat_bds[..., 0] +# # +# # Replace quiver here with delta-based lineplot +# # +# urd = xbl + deg_scale * u0r +# vrd = ybl + deg_scale * v0r * yscale +# for iy in range(ny): +# for ix in range(nx): +# plt.plot([xbl[iy, ix], urd[iy, ix]], +# [ybl[iy, ix], vrd[iy, ix]], +# ':', color='magenta', linewidth=2.5, +# transform=data_proj) +# # Show this is the SAME as lineplot +# plt.quiver(xbl, ybl, u0r, v0r * yscale, +# color='magenta', width=0.01, +# headwidth=1.2, # headlength=1.0, headaxislength=0.7, +# angles=angles, +# scale_units='xy', scale=scale, +# transform=data_proj) +# +# plt.suptitle('angles from "{}"'.format(cell_angle_boundpoints)) +# +## # Also draw small lines pointing at the correct (TRUE, not ) angle. +## ny, nx = x0s.shape +## size_degrees = 1.0 +## angles = angles_cube.copy() +## angles.convert_units('radians') +## angles = angles.data +## lats = testcube.coord('latitude').copy() +## lats.convert_units('radians') +## lats = lats.points +## dxs = size_degrees * u0.copy() #* np.cos(angles) +## dys = size_degrees * u0.copy() # / np.sqrt(np.cos(lats)) +## x1s = x0s + dxs +## y1s = y0s + dys +### for iy in range(ny): +### for ix in range(nx): +### plt.plot([x0s[iy, ix], x1s[iy, ix]], +### [y0s[iy, ix], y1s[iy, ix]], +### 'o-', markersize=4., markeredgewidth=0., +### color='green', # scale_units='xy', scale=scale, +### transform=data_proj) +## plt.quiver(x0s, y0s, dxs, dys, +## color='green', linewidth=0.2, +## angles=angles, +## scale_units='xy', scale=scale * 0.6, +## transform=data_proj) +# +# +# +# ax.set_global() +# plt.show() +# +# angles_cube.convert_units('degrees') +# +# self.assertArrayAllClose( +# angles_cube.data, +# [[33.421, 17.928, 0., -17.928, -33.421], +# [41.981, 24.069, 0., -24.069, -41.981], +# [56.624, 37.809, 0., -37.809, -56.624], +# [79.940, 74.227, 0., -74.227, -79.940], +# [107.313, 126.361, -180., -126.361, -107.313]], +# atol=0.002) + + +if __name__ == "__main__": + tests.main() From 670982e3ef8c4232f9222829e5332dba775cc7c9 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Sun, 29 Jul 2018 17:19:01 +0100 Subject: [PATCH 02/15] Enhanced testing; better checking and crs awareness in grid_angles routine. --- lib/iris/analysis/_grid_angles.py | 46 ++++++++++---- .../cartography/test_gridcell_angles.py | 63 ++++++++++++------- 2 files changed, 76 insertions(+), 33 deletions(-) diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py index f847e0d00f..d371c8f1f9 100644 --- a/lib/iris/analysis/_grid_angles.py +++ b/lib/iris/analysis/_grid_angles.py @@ -24,6 +24,7 @@ import numpy as np +import cartopy.crs as ccrs import iris @@ -183,8 +184,6 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): 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') 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 {}.') @@ -193,19 +192,42 @@ 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 None: + # No coord system : assume x + y are lons + lats. + # Just make sure they are in degrees ! + # NOTE: raises an error if units are not angles. + x = x_coord.units.convert(x, 'degrees') + y = y_coord.units.convert(y, 'degrees') + else: + # 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') @@ -215,8 +237,8 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): raise ValueError(*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 adjacent points with - # increasing longitudes. + # 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: # Data is points arrays. diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py index 26de728dad..ce783e6b91 100644 --- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py +++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py @@ -115,10 +115,10 @@ def _singlecell_diamond_cube(self, x0=90., y0=0., dy=10., dx_eq=None): cube.add_aux_coord(co_y, (0, 1)) return cube - def _check_different_orientations_and_latitudes( - self, - method='mid-lhs, mid-rhs', - atol_degrees=0.005): + def _check_orientations_and_latitudes(self, + method='mid-lhs, mid-rhs', + atol_degrees=0.005, + cellsize_degrees=1.0): ny, nx = 7, 9 x0, x1 = -164, 164 y0, y1 = -75, 75 @@ -128,8 +128,9 @@ def _check_different_orientations_and_latitudes( # Make gridcells rectangles surrounding these centrepoints, but also # tilted at various angles (= same as x-point lons, as that's easy). - dx = 1.0 # half-width of gridcells, in degrees - dy = dx # half-height of gridcells, in degrees +# dx = cellsize_degrees # half-width of gridcells, in degrees +# dy = dx # half-height of gridcells, in degrees + # 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)] @@ -137,10 +138,10 @@ def _check_different_orientations_and_latitudes( 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 = dx * np.cos(xangs) * dx_corners - x_ofs_2d -= dy * np.sin(xangs) * dy_corners - y_ofs_2d = dy * np.cos(xangs) * dy_corners - y_ofs_2d += dx * np.sin(xangs) * dx_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. @@ -171,22 +172,42 @@ def _check_different_orientations_and_latitudes( angles_expected[:] = angles self.assertArrayAllClose(angles_calculated, angles_expected, atol=atol_degrees) - return angles_calculated + return angles_calculated, angles_expected - def test_different_orientations_and_latitudes(self): - self._check_different_orientations_and_latitudes() + def test_all_orientations_and_latitudes(self): + self._check_orientations_and_latitudes() def test_different_methods(self): - r1 = self._check_different_orientations_and_latitudes( - 'mid-lhs, mid-rhs', atol_degrees=5.0) - r2 = self._check_different_orientations_and_latitudes( - 'lower-left, lower-right', atol_degrees=5.0) -# print(np.round(r1 - r2, 3)) -# for fn in (np.max, np.mean): -# print(fn.__name__, fn(np.abs(r1 - r2))) - atol = 1.0 # !!!!!!! + # Get results with both calculation methods. + # A smallish cellsize should yield similar results in both cases. + r1, _ = self._check_orientations_and_latitudes( + method='mid-lhs, mid-rhs', + cellsize_degrees=0.1, atol_degrees=0.1) + r2, _ = self._check_orientations_and_latitudes( + method='lower-left, lower-right', + cellsize_degrees=0.1, atol_degrees=0.1) + + print(np.round(r1 - r2, 3)) + for fn in (np.max, np.mean): + print(fn.__name__, fn(np.abs(r1 - r2))) + + atol = 0.1 # A whole degree - significantly different at higher latitudes. self.assertArrayAllClose(r1, r2, atol=atol) + def test_methods_and_cellsizes(self): + for cellsize in (0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0): + r_mid, exp_mid = self._check_orientations_and_latitudes( + method='mid-lhs, mid-rhs', + cellsize_degrees=cellsize, atol_degrees=25) + r_btm, exp_btm = self._check_orientations_and_latitudes( + method='lower-left, lower-right', + cellsize_degrees=cellsize, atol_degrees=25) + wc_mid = np.max(np.abs(r_mid - exp_mid)) + wc_btm = np.max(np.abs(r_btm - exp_btm)) + msg = ('Cell size = {:5.2f} degrees, wc-abs-errors : ' + 'mid-lr={:7.3f} lower-lr={:7.3f}') + print(msg.format(cellsize, wc_mid, wc_btm)) + # def test_single_cell_equatorial(self): # plt.switch_backend('tkagg') From fdb693b7ce9fb69777c0452dd1daf5a5bf7f7adf Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Fri, 3 Aug 2018 17:03:55 +0100 Subject: [PATCH 03/15] Remove crud from test_gridcell_angles. --- .../cartography/test_gridcell_angles.py | 410 +----------------- 1 file changed, 21 insertions(+), 389 deletions(-) diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py index ce783e6b91..4023588dc4 100644 --- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py +++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py @@ -27,98 +27,19 @@ import iris.tests as tests import numpy as np -import numpy.ma as ma -import cartopy.crs as ccrs from iris.cube import Cube -from iris.coords import DimCoord, AuxCoord -import iris.coord_systems -from iris.analysis.cartography import unrotate_pole +from iris.coords import AuxCoord -from iris.analysis.cartography import (gridcell_angles, - rotate_grid_vectors) - -import matplotlib.pyplot as plt -from orca_utils.plot_testing.blockplot_from_bounds import blockplot_2dll - - -def _rotated_grid_sample(pole_lat=15, pole_lon=-180, - lon_bounds=np.linspace(-30, 30, 6, endpoint=True), - lat_bounds=np.linspace(-30, 30, 6, endpoint=True)): - # Calculate *true* lat_bounds+lon_bounds for the rotated grid. - lon_bounds = np.array(lon_bounds, dtype=float) - lat_bounds = np.array(lat_bounds, dtype=float) - # Construct centrepoints. - lons = 0.5 * (lon_bounds[:-1] + lon_bounds[1:]) - lats = 0.5 * (lat_bounds[:-1] + lat_bounds[1:]) - # Convert all to full 2d arrays. - lon_bounds, lat_bounds = np.meshgrid(lon_bounds, lat_bounds) - lons, lats = np.meshgrid(lons, lats) - # Calculate true lats+lons for all points. - lons_true_bds, lats_true_bds = unrotate_pole(lon_bounds, lat_bounds, - pole_lon, pole_lat) - lons_true, lats_true = unrotate_pole(lons, lats, pole_lon, pole_lat) - # Make the 'unified' bounds into contiguous (ny, nx, 4) arrays. - def expand_unified_bds(bds): - ny, nx = bds.shape - bds_4 = np.zeros((ny - 1, nx - 1, 4)) - bds_4[:, :, 0] = bds[:-1, :-1] - bds_4[:, :, 1] = bds[:-1, 1:] - bds_4[:, :, 2] = bds[1:, 1:] - bds_4[:, :, 3] = bds[1:, :-1] - return bds_4 - - lon_true_bds4, lat_true_bds4 = (expand_unified_bds(bds) - for bds in (lons_true_bds, lats_true_bds)) - # Make these into a 2d-latlon grid for a cube - cube = Cube(np.zeros(lon_true_bds4.shape[:-1])) - co_x = AuxCoord(lons_true, bounds=lon_true_bds4, - standard_name='longitude', units='degrees') - co_y = AuxCoord(lats_true, bounds=lat_true_bds4, - standard_name='latitude', units='degrees') - cube.add_aux_coord(co_x, (0, 1)) - cube.add_aux_coord(co_y, (0, 1)) - return cube +from iris.analysis.cartography import gridcell_angles class TestGridcellAngles(tests.IrisTest): - def _singlecell_30deg_cube(self, x0=90., y0=0., dx=20., dy=10.): - x_pts = np.array([[x0]]) - y_pts = np.array([[y0]]) - x_bds = x0 + dx * np.array([[[-1., 1, 0.5, -1.5]]]) -# self.assertArrayAllClose(x_bds, np.array([[[70., 110, 100, 60]]])) - y_bds = y0 + dy * np.array([[[-1., 1, 3, 1]]]) -# self.assertArrayAllClose(y_bds, np.array([[[-10., 10, 30, 10]]])) - co_x = AuxCoord(points=x_pts, bounds=x_bds, - standard_name='longitude', units='degrees') - co_y = AuxCoord(points=y_pts, bounds=y_bds, - standard_name='latitude', units='degrees') - cube = Cube(np.zeros((1, 1))) - cube.add_aux_coord(co_x, (0, 1)) - cube.add_aux_coord(co_y, (0, 1)) - return cube - - def _singlecell_diamond_cube(self, x0=90., y0=0., dy=10., dx_eq=None): - if dx_eq is None: - dx_eq = dy - x_pts = np.array([[x0]]) - y_pts = np.array([[y0]]) - dx = dx_eq / np.cos(np.deg2rad(y0)) - x_bds = np.array([[[x0, x0 + dx, x0, x0 - dx]]]) - y_bds = np.array([[[y0 - dy, y0, y0 + dy, y0]]]) - co_x = AuxCoord(points=x_pts, bounds=x_bds, - standard_name='longitude', units='degrees') - co_y = AuxCoord(points=y_pts, bounds=y_bds, - standard_name='latitude', units='degrees') - cube = Cube(np.zeros((1, 1))) - cube.add_aux_coord(co_x, (0, 1)) - cube.add_aux_coord(co_y, (0, 1)) - return cube - - def _check_orientations_and_latitudes(self, - method='mid-lhs, mid-rhs', - atol_degrees=0.005, - cellsize_degrees=1.0): + def _check_multiple_orientations_and_latitudes( + self, + method='mid-lhs, mid-rhs', + atol_degrees=0.005, + cellsize_degrees=1.0): ny, nx = 7, 9 x0, x1 = -164, 164 y0, y1 = -75, 75 @@ -166,324 +87,35 @@ def _check_orientations_and_latitudes(self, # itself are exact : Errors scale as the square of gridcell sizes. angles_cube.convert_units('degrees') angles_calculated = angles_cube.data + # Note: expand the original 1-d test angles into the full result shape, # just to please 'np.testing.assert_allclose', which doesn't broadcast. angles_expected = np.zeros(angles_cube.shape) angles_expected[:] = angles + + # Assert (toleranced) equality, and return results. self.assertArrayAllClose(angles_calculated, angles_expected, atol=atol_degrees) + return angles_calculated, angles_expected - def test_all_orientations_and_latitudes(self): - self._check_orientations_and_latitudes() + def test_various_orientations_and_locations(self): + self._check_multiple_orientations_and_latitudes() - def test_different_methods(self): - # Get results with both calculation methods. + 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_orientations_and_latitudes( - method='mid-lhs, mid-rhs', - cellsize_degrees=0.1, atol_degrees=0.1) - r2, _ = self._check_orientations_and_latitudes( + 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) - print(np.round(r1 - r2, 3)) - for fn in (np.max, np.mean): - print(fn.__name__, fn(np.abs(r1 - r2))) - - atol = 0.1 # A whole degree - significantly different at higher latitudes. + # They are not the same : checks we selected the 'other' method ! + self.assertFalse(np.allclose(r1, r2)) + # Note: results are rather different at higher latitudes. + atol = 0.1 self.assertArrayAllClose(r1, r2, atol=atol) - def test_methods_and_cellsizes(self): - for cellsize in (0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0): - r_mid, exp_mid = self._check_orientations_and_latitudes( - method='mid-lhs, mid-rhs', - cellsize_degrees=cellsize, atol_degrees=25) - r_btm, exp_btm = self._check_orientations_and_latitudes( - method='lower-left, lower-right', - cellsize_degrees=cellsize, atol_degrees=25) - wc_mid = np.max(np.abs(r_mid - exp_mid)) - wc_btm = np.max(np.abs(r_btm - exp_btm)) - msg = ('Cell size = {:5.2f} degrees, wc-abs-errors : ' - 'mid-lr={:7.3f} lower-lr={:7.3f}') - print(msg.format(cellsize, wc_mid, wc_btm)) - - -# def test_single_cell_equatorial(self): -# plt.switch_backend('tkagg') -# plt.figure(figsize=(10,10)) -## ax = plt.axes(projection=ccrs.Mercator()) -## ax = plt.axes(projection=ccrs.NorthPolarStereo()) -# ax = plt.axes(projection=ccrs.Orthographic(central_longitude=90., -# central_latitude=30.)) -# -# lon0 = 90.0 -# dy = 1.0 -# dx = 3.0 -# y_0, y_n, ny = -80, 80, 9 -# angles = [] -# for lat in np.linspace(y_0, y_n, ny): -# cube = self._singlecell_diamond_cube(x0=lon0, y0=lat, -# dy=dy, dx_eq=dx) -# angles_cube = gridcell_angles(cube, -## cell_angle_boundpoints='mid-lhs, mid-rhs') -# cell_angle_boundpoints='lower-left, lower-right') -# tmp_cube = angles_cube.copy() -# tmp_cube.convert_units('degrees') -## print('') -## print(lat) -## co_x, co_y = (cube.coord(axis=ax) for ax in ('x', 'y')) -## print() -## print(' at : {}, {}'.format(co_x.points[0, 0], co_y.points[0, 0])) -## print(' x-bds:') -## print(co_x.bounds) -## print(' y-bds:') -## print(co_y.bounds) -# angle = tmp_cube.data[0, 0] -# angles.append(angle) -# print(lat, angle) -# blockplot_2dll(cube) -# -# ax.coastlines() -# ax.set_global() -# -# # Plot constant NEly (45deg) arrows. -# xx = np.array([lon0] * ny) -# yy = np.linspace(y_0, y_n, ny) - dy -# uu = np.array([1.0] * ny) -# plt.quiver(xx, yy, -# uu, np.cos(np.deg2rad(yy)), -# zorder=2, color='red', -## scale_units='xy', -# angles='uv', -# transform=ccrs.PlateCarree()) -# -# # Also plot returned angles. -# angles_arr_rad = np.deg2rad(angles) -# u_arr = uu * np.cos(angles_arr_rad) -# v_arr = uu * np.sin(angles_arr_rad) * np.cos(np.deg2rad(yy)) -# -# plt.quiver(xx, yy, -# u_arr, -# v_arr, -# zorder=2, color='magenta', -# scale_units='xy', -# width=0.005, -# scale=0.2e-6, -## width=0.5, -# transform=ccrs.PlateCarree()) -# -# plt.show() - - -# def test_values(self): -# # Construct a rotated-pole grid and check angle calculation. -# testcube = _rotated_grid_sample() -# -# cell_angle_boundpoints = 'mid-lhs, mid-rhs' -## cell_angle_boundpoints = 'lower-left, lower-right' -## cell_angle_boundpoints = 'garble' -# angles_cube = gridcell_angles( -# testcube, -# cell_angle_boundpoints=cell_angle_boundpoints) -# angles_cube.convert_units('radians') -# -# # testing phase... -# print(np.rad2deg(angles_cube.data)) -# -# import matplotlib.pyplot as plt -# plt.switch_backend('tkagg') -# -## plot_map = 'north_polar_stereographic' -## plot_map = 'plate_carree' -## plot_map = 'mercator' -# plot_map = 'north_polar_orthographic' -# if plot_map == 'plate_carree': -# scale = 0.1 -# map_proj = ccrs.PlateCarree() -# elif plot_map == 'mercator': -# scale = 3.0e-6 -# map_proj = ccrs.Mercator() -# map_proj._threshold *= 0.01 -# elif plot_map == 'north_polar_orthographic': -# scale = 3.0e-6 -# map_proj = ccrs.Orthographic(central_longitude=0.0, -# central_latitude=90.0,) -# map_proj._threshold *= 0.01 -# elif plot_map == 'north_polar_stereographic': -# scale = 3.0e-6 -# map_proj = ccrs.NorthPolarStereo() -# else: -# assert 0 -# -# ax = plt.axes(projection=map_proj) -# data_proj = ccrs.PlateCarree() -# -# deg_scale = 10.0 -# -## angles = 'uv' -# angles = 'xy' -# -# ax.coastlines() -# ax.gridlines() -# for i_bnd in range(4): -# color = ['black', 'red', 'blue', 'magenta'][i_bnd] -# plt.plot(testcube.coord('longitude').bounds[..., i_bnd], -# testcube.coord('latitude').bounds[..., i_bnd], -# '+', markersize=10., markeredgewidth=2., -# markerfacecolor=color, markeredgecolor=color, -# transform=data_proj) -# -# -# # Show plain 0,1 + 1,0 (PlateCarree) vectors unrotated at the given points. -# pts_shape = testcube.coord('longitude').shape -# ny, nx = pts_shape -# u0 = np.ones(pts_shape) -# v0 = np.zeros(pts_shape) -# u1 = v0.copy() -# v1 = u0.copy() -# -# x0s = testcube.coord('longitude').points -# y0s = testcube.coord('latitude').points -# yscale = np.cos(np.deg2rad(y0s)) -# plt.quiver(x0s, y0s, u0, v0 * yscale, -# color='blue', width=0.005, -# headwidth=2., # headlength=1.0, headaxislength=0.7, -# angles=angles, -# scale_units='xy', scale=scale, -# transform=data_proj) -# plt.quiver(x0s, y0s, u1, v1 * yscale, -# color='red', width=0.005, -# headwidth=2., # headlength=1.0, headaxislength=0.7, -# angles=angles, -# scale_units='xy', scale=scale, -# transform=data_proj) -# -# # Add 45deg arrows (NEly), still on a PlateCarree map. -# plt.quiver(x0s, y0s, v1, v1 * yscale, -# color='green', width=0.005, -# headwidth=2., # headlength=1.0, headaxislength=0.7, -# angles=angles, -# scale_units='xy', scale=scale, -# transform=data_proj) -# -# -# -# # -# # Repeat the above plotting short lines INSTEAD of quiver. -# # -# u0d = x0s + deg_scale * u0 -# v0d = y0s + deg_scale * v0 -# u1d = x0s + deg_scale * u1 -# v1d = y0s + deg_scale * v1 -# u2d = x0s + deg_scale * u0 -# v2d = y0s + deg_scale * v1 -# for iy in range(ny): -# for ix in range(nx): -# plt.plot([x0s[iy, ix], u0d[iy, ix]], -# [y0s[iy, ix], v0d[iy, ix]], -# ':', color='blue', linewidth=0.5, -# transform=data_proj) -# plt.plot([x0s[iy, ix], u1d[iy, ix]], -# [y0s[iy, ix], v1d[iy, ix]], -# ':', color='red', linewidth=0.5, -# transform=data_proj) -# plt.plot([x0s[iy, ix], u2d[iy, ix]], -# [y0s[iy, ix], v2d[iy, ix]], -# ':', color='green', linewidth=0.5, -# transform=data_proj) -# -# -# # Overplot BL-BR and BL-TL lines from the cell bounds. -# co_lon, co_lat = [testcube.coord(name).copy() -# for name in ('longitude', 'latitude')] -# for co in (co_lon, co_lat): -# co.convert_units('degrees') -# lon_bds, lat_bds = [co.bounds for co in (co_lon, co_lat)] -## ny, nx = lon_bds.shape[:-1] -# for iy in range(ny): -# for ix in range(nx): -# x0, y0 = lon_bds[iy, ix, 0], lat_bds[iy, ix, 0] -# x1, y1 = lon_bds[iy, ix, 1], lat_bds[iy, ix, 1] -# x2, y2 = lon_bds[iy, ix, 3], lat_bds[iy, ix, 3] -# plt.plot([x0, x1], [y0, y1], 'x-', -# color='orange', -# transform=data_proj) -# plt.plot([x0, x2], [y0, y2], 'x-', -# color='orange', linestyle='--', -# transform=data_proj) -# -# # Plot U0, rotated by cell angles, also at cell bottom-lefts. -# u0_cube, u1_cube, v0_cube, v1_cube = [testcube.copy(data=aa) -# for aa in (u0, v0, u1, v1)] -# u0r_cube, v0r_cube = rotate_grid_vectors( -# u0_cube, v0_cube, grid_angles_cube=angles_cube) -# u0r, v0r = [cube.data for cube in (u0r_cube, v0r_cube)] -# -# xbl, ybl = lon_bds[..., 0], lat_bds[..., 0] -# # -# # Replace quiver here with delta-based lineplot -# # -# urd = xbl + deg_scale * u0r -# vrd = ybl + deg_scale * v0r * yscale -# for iy in range(ny): -# for ix in range(nx): -# plt.plot([xbl[iy, ix], urd[iy, ix]], -# [ybl[iy, ix], vrd[iy, ix]], -# ':', color='magenta', linewidth=2.5, -# transform=data_proj) -# # Show this is the SAME as lineplot -# plt.quiver(xbl, ybl, u0r, v0r * yscale, -# color='magenta', width=0.01, -# headwidth=1.2, # headlength=1.0, headaxislength=0.7, -# angles=angles, -# scale_units='xy', scale=scale, -# transform=data_proj) -# -# plt.suptitle('angles from "{}"'.format(cell_angle_boundpoints)) -# -## # Also draw small lines pointing at the correct (TRUE, not ) angle. -## ny, nx = x0s.shape -## size_degrees = 1.0 -## angles = angles_cube.copy() -## angles.convert_units('radians') -## angles = angles.data -## lats = testcube.coord('latitude').copy() -## lats.convert_units('radians') -## lats = lats.points -## dxs = size_degrees * u0.copy() #* np.cos(angles) -## dys = size_degrees * u0.copy() # / np.sqrt(np.cos(lats)) -## x1s = x0s + dxs -## y1s = y0s + dys -### for iy in range(ny): -### for ix in range(nx): -### plt.plot([x0s[iy, ix], x1s[iy, ix]], -### [y0s[iy, ix], y1s[iy, ix]], -### 'o-', markersize=4., markeredgewidth=0., -### color='green', # scale_units='xy', scale=scale, -### transform=data_proj) -## plt.quiver(x0s, y0s, dxs, dys, -## color='green', linewidth=0.2, -## angles=angles, -## scale_units='xy', scale=scale * 0.6, -## transform=data_proj) -# -# -# -# ax.set_global() -# plt.show() -# -# angles_cube.convert_units('degrees') -# -# self.assertArrayAllClose( -# angles_cube.data, -# [[33.421, 17.928, 0., -17.928, -33.421], -# [41.981, 24.069, 0., -24.069, -41.981], -# [56.624, 37.809, 0., -37.809, -56.624], -# [79.940, 74.227, 0., -74.227, -79.940], -# [107.313, 126.361, -180., -126.361, -107.313]], -# atol=0.002) - if __name__ == "__main__": tests.main() From e832fe387664886b95c9574313e0788ed157746d Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Fri, 3 Aug 2018 19:27:40 +0100 Subject: [PATCH 04/15] Use degree units for everything in _grid_angles. --- lib/iris/analysis/_grid_angles.py | 35 +++++++++++++++++-------------- 1 file changed, 19 insertions(+), 16 deletions(-) diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py index d371c8f1f9..e5179b05f9 100644 --- a/lib/iris/analysis/_grid_angles.py +++ b/lib/iris/analysis/_grid_angles.py @@ -66,13 +66,13 @@ def _latlon_from_xyz(xyz): Returns: lonlat : (array) - spherical angles, of dims (2, ), in radians. + spherical angles, of dims (2, ), 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]) @@ -103,7 +103,7 @@ def _angle(p, q, r): psi = np.arccos(cosine) * np.sign(r[1] - p[1]) psi[index] = np.nan - return psi + return np.rad2deg(psi) def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): @@ -160,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 @@ -184,6 +184,15 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): 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 {}.') @@ -211,13 +220,7 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): x, y = x_coord.points, y_coord.points # Make sure these arrays are ordinary lats+lons, in degrees. - if cs is None: - # No coord system : assume x + y are lons + lats. - # Just make sure they are in degrees ! - # NOTE: raises an error if units are not angles. - x = x_coord.units.convert(x, 'degrees') - y = y_coord.units.convert(y, 'degrees') - else: + if cs is not None: # Transform points into true lats + lons. crs_src = cs.as_cartopy_crs() crs_pc = ccrs.PlateCarree() @@ -262,7 +265,7 @@ 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: - # Data is points arrays. + # 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) @@ -296,8 +299,8 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): 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)) + 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]) @@ -305,7 +308,7 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): 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 From 56944dfe4c81e3bb2558693a90d5a900e904f753 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Fri, 3 Aug 2018 19:39:02 +0100 Subject: [PATCH 05/15] Make assertArrayAllClose print details when it fails. --- lib/iris/tests/__init__.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/lib/iris/tests/__init__.py b/lib/iris/tests/__init__.py index 710a3fbac6..c6fb748740 100644 --- a/lib/iris/tests/__init__.py +++ b/lib/iris/tests/__init__.py @@ -663,6 +663,22 @@ def assertArrayAllClose(self, a, b, rtol=1.0e-7, atol=0.0, **kwargs): For full details see underlying routine numpy.testing.assert_allclose. """ + ok = np.allclose(a, b, rtol=rtol, atol=atol, **kwargs) + if not ok: + # Print out the worst-case failure values. + # 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), errors.shape) + print('') + print('FAILED ARRAY CHECK "assertArrayAllClose" :') + msg = ' with shapes={} {}, atol={}, rtol={},' + print(msg.format(a.shape, b.shape, atol, rtol)) + aval, bval = a[worst_inds], b[worst_inds] + absdiff = np.abs(aval - bval) + print(' worst at element {} : abs({} - {}) = {} '.format( + worst_inds, aval, bval, absdiff)) + np.testing.assert_allclose(a, b, rtol=rtol, atol=atol, **kwargs) @contextlib.contextmanager From f1493411bcdcc35095a8b8a693552a89aa812e01 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Fri, 3 Aug 2018 19:47:23 +0100 Subject: [PATCH 06/15] Rework and extend testing for gridcell_angles. --- .../cartography/test_gridcell_angles.py | 218 ++++++++++++++---- 1 file changed, 168 insertions(+), 50 deletions(-) diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py index 4023588dc4..14b1fb0bf4 100644 --- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py +++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py @@ -32,51 +32,71 @@ from iris.coords import AuxCoord from iris.analysis.cartography import gridcell_angles +from iris.tests.stock import sample_2d_latlons + + +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 _check_multiple_orientations_and_latitudes( - self, - method='mid-lhs, mid-rhs', - atol_degrees=0.005, - cellsize_degrees=1.0): - ny, nx = 7, 9 - x0, x1 = -164, 164 - y0, y1 = -75, 75 - lats = np.linspace(y0, y1, ny, endpoint=True) - angles = np.linspace(x0, x1, nx, endpoint=True) - x_pts_2d, y_pts_2d = np.meshgrid(angles, lats) - - # Make gridcells rectangles surrounding these centrepoints, but also - # tilted at various angles (= same as x-point lons, as that's easy). -# dx = cellsize_degrees # half-width of gridcells, in degrees -# dy = dx # half-height of gridcells, in degrees - - # 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)) + 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_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) @@ -88,12 +108,25 @@ def _check_multiple_orientations_and_latitudes( angles_cube.convert_units('degrees') angles_calculated = angles_cube.data - # Note: expand the original 1-d test angles into the full result shape, - # just to please 'np.testing.assert_allclose', which doesn't broadcast. - angles_expected = np.zeros(angles_cube.shape) - angles_expected[:] = angles + # 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. + ok = np.allclose(angles_calculated, angles_expected, + atol=atol_degrees) + if not ok: + print('FAIL') + diffs = angles_calculated - angles_expected + worst_inds = np.unravel_index(np.argmax(np.abs(diffs)), + angles_calculated.shape) + print('max abs-error : got({}) - expected({}) at {}'.format( + angles_calculated[worst_inds], angles_expected[worst_inds], + worst_inds)) self.assertArrayAllClose(angles_calculated, angles_expected, atol=atol_degrees) @@ -110,11 +143,96 @@ def test_bottom_edge_method(self): method='lower-left, lower-right', cellsize_degrees=0.1, atol_degrees=0.1) - # They are not the same : checks we selected the 'other' method ! + # Not *exactly* the same : this checks we tested the 'other' method ! self.assertFalse(np.allclose(r1, r2)) - # Note: results are rather different at higher latitudes. - atol = 0.1 - self.assertArrayAllClose(r1, r2, atol=atol) + # 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 gets 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_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_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 values alone. + 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 > 100 degrees ! if __name__ == "__main__": From dff77502d2df403ee16ca566d76c75a37b28bca8 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Mon, 13 Aug 2018 10:59:33 +0100 Subject: [PATCH 07/15] Fix assertArrayAllClose; remove debug code from test_gridcell_angles. --- lib/iris/tests/__init__.py | 22 +++++++++---------- .../cartography/test_gridcell_angles.py | 8 ------- 2 files changed, 11 insertions(+), 19 deletions(-) diff --git a/lib/iris/tests/__init__.py b/lib/iris/tests/__init__.py index c6fb748740..7550fce1ca 100644 --- a/lib/iris/tests/__init__.py +++ b/lib/iris/tests/__init__.py @@ -663,23 +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. """ - ok = np.allclose(a, b, rtol=rtol, atol=atol, **kwargs) + ok = np.allclose(a, b, atol=atol, rtol=rtol, equal_nan=True) if not ok: - # Print out the worst-case failure values. # 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), errors.shape) - print('') - print('FAILED ARRAY CHECK "assertArrayAllClose" :') - msg = ' with shapes={} {}, atol={}, rtol={},' - print(msg.format(a.shape, b.shape, atol, rtol)) + 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) - print(' worst at element {} : abs({} - {}) = {} '.format( - worst_inds, aval, bval, absdiff)) - - np.testing.assert_allclose(a, b, rtol=rtol, atol=atol, **kwargs) + 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=''): diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py index 14b1fb0bf4..3fce8c29fa 100644 --- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py +++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py @@ -119,14 +119,6 @@ def _check_multiple_orientations_and_latitudes(self, # Assert (toleranced) equality, and return results. ok = np.allclose(angles_calculated, angles_expected, atol=atol_degrees) - if not ok: - print('FAIL') - diffs = angles_calculated - angles_expected - worst_inds = np.unravel_index(np.argmax(np.abs(diffs)), - angles_calculated.shape) - print('max abs-error : got({}) - expected({}) at {}'.format( - angles_calculated[worst_inds], angles_expected[worst_inds], - worst_inds)) self.assertArrayAllClose(angles_calculated, angles_expected, atol=atol_degrees) From fd14631316d2df8a30fba20f88f5447de6362282 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 14 Aug 2018 15:43:14 +0100 Subject: [PATCH 08/15] Remove obsolete assignments. --- lib/iris/analysis/cartography.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/lib/iris/analysis/cartography.py b/lib/iris/analysis/cartography.py index af55a9de57..a778f5f846 100644 --- a/lib/iris/analysis/cartography.py +++ b/lib/iris/analysis/cartography.py @@ -39,10 +39,6 @@ from iris.util import _meshgrid from ._grid_angles import gridcell_angles, rotate_grid_vectors -# Avoid pycodestyle warnings for unused imports. -gridcell_angles = gridcell_angles -rotate_grid_vectors = rotate_grid_vectors - # List of contents to control Sphinx autodocs. # Unfortunately essential to get docs for the grid_angles functions. __all__ = [ From d19f3d99414b71e374865873033b2f7a4b4c1482 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 14 Aug 2018 15:48:15 +0100 Subject: [PATCH 09/15] Remove obsolete code. --- .../tests/unit/analysis/cartography/test_gridcell_angles.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py index 3fce8c29fa..fcb22ae090 100644 --- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py +++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py @@ -117,8 +117,6 @@ def _check_multiple_orientations_and_latitudes(self, angles_expected = (angles_expected + 360.) % 360. # Assert (toleranced) equality, and return results. - ok = np.allclose(angles_calculated, angles_expected, - atol=atol_degrees) self.assertArrayAllClose(angles_calculated, angles_expected, atol=atol_degrees) From 330c71e3a7115f80c7b099a508165728589f8361 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 14 Aug 2018 16:03:44 +0100 Subject: [PATCH 10/15] Small comment improvements. --- lib/iris/tests/stock/_stock_2d_latlons.py | 16 ++++++++++------ .../analysis/cartography/test_gridcell_angles.py | 6 +++--- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/lib/iris/tests/stock/_stock_2d_latlons.py b/lib/iris/tests/stock/_stock_2d_latlons.py index 75b62f31ba..395a50d8ea 100644 --- a/lib/iris/tests/stock/_stock_2d_latlons.py +++ b/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. diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py index fcb22ae090..45ef52f607 100644 --- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py +++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py @@ -139,7 +139,7 @@ def test_bottom_edge_method(self): self.assertArrayAllClose(r1, r2, atol=0.1) def test_bounded_coord_args(self): - # Check that passing the coords gets the same result as the cube. + # 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) @@ -188,7 +188,7 @@ def test_unbounded_regional_coord_args(self): self.standard_small_cube_results[:, 1:-1]) def test_points_array_args(self): - # Check we can calculate from points values alone. + # 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. @@ -222,7 +222,7 @@ def test_unbounded_global(self): 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 > 100 degrees ! + # maximum errors there in the end columns are actually > 100 degrees ! if __name__ == "__main__": From 995864de1582f0aac5f8ccc5402cae140b1dc343 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 15 Aug 2018 17:00:43 +0100 Subject: [PATCH 11/15] Attempt to clarify docstrings of low-level routines. --- lib/iris/analysis/_grid_angles.py | 60 +++++++++++++++++++++++-------- 1 file changed, 45 insertions(+), 15 deletions(-) diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py index e5179b05f9..0492bf7143 100644 --- a/lib/iris/analysis/_grid_angles.py +++ b/lib/iris/analysis/_grid_angles.py @@ -34,12 +34,16 @@ def _3d_xyz_from_latlon(lon, lat): Args: - * lon, lat: (arrays in degrees) + * 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. Dimension 0 maps x,y,z. + * 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) @@ -61,13 +65,16 @@ def _latlon_from_xyz(xyz): Args: * xyz: (array) - positions array, of dims (3, ), where index 0 maps x/y/z. + Array of 3-D cartesian coordinates. + Shape (3, ). + x / y / z values are in xyz[0 / 1 / 2], Returns: - lonlat : (array) - spherical angles, of dims (2, ), in degrees. - Dim 0 maps longitude, latitude. + * 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])) @@ -78,14 +85,37 @@ def _latlon_from_xyz(xyz): def _angle(p, q, r): """ - Return angle (in _radians_) of grid wrt local east. - Anticlockwise +ve, as usual. - {P, Q, R} are consecutive points in the same row, - eg {v(i,j),f(i,j),v(i+1,j)}, or {T(i-1,j),T(i,j),T(i+1,j)} - Calculate dot product of PR with lambda_hat at Q. - This gives us cos(required angle). - Disciminate between +/- angles by comparing latitudes of P and R. - p, q, r, are all 2-element arrays [lon, lat] of angles in degrees. + 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. + + 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]) From 4820572e1037c9fa5e57da8deefba6b6d22a69cd Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 15 Aug 2018 19:03:20 +0100 Subject: [PATCH 12/15] More tests, and some functional fixes. --- lib/iris/analysis/_grid_angles.py | 48 +++++++--- .../cartography/test_gridcell_angles.py | 93 +++++++++++++++++-- 2 files changed, 117 insertions(+), 24 deletions(-) diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py index 0492bf7143..a9d31412c5 100644 --- a/lib/iris/analysis/_grid_angles.py +++ b/lib/iris/analysis/_grid_angles.py @@ -205,7 +205,7 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): x, y = cube.coord(axis='x'), cube.coord(axis='y') # Now should have either 2 coords or 2 arrays. - if not hasattr(x, 'shape') and hasattr(y, 'shape'): + 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))) @@ -224,11 +224,11 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): coord.convert_units('degrees') if x_coord.ndim != 2 or y_coord.ndim != 2: - msg = ('Coordinate inputs must have 2-dimensional shape. ', + 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. ', + 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: @@ -239,7 +239,7 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): 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. ', + msg = ('Coordinate inputs must have same coordinate system. ' 'Got x of {} and y of {}.') raise ValueError(msg.format(cs, y_coord.coord_system)) @@ -254,12 +254,30 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): # 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) + 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') + x_coord = iris.coords.AuxCoord( + points=xpts, bounds=xbds, + standard_name='longitude', units='degrees') elif hasattr(x, 'bounds') or hasattr(y, 'bounds'): # One was a Coord, and the other not ? @@ -267,15 +285,15 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): if hasattr(y, 'bounds'): is_and_not = reversed(is_and_not) msg = 'Input {!r} is a Coordinate, but {!r} is not.' - raise ValueError(*is_and_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 X - # indices. + # 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 longitude-axis as references. + # 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' @@ -296,7 +314,7 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): units='degrees') else: # Data is bounds arrays. - # Use different bounds points in the longitude-axis as references. + # 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. diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py index 45ef52f607..289ffab8c1 100644 --- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py +++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py @@ -32,7 +32,7 @@ from iris.coords import AuxCoord from iris.analysis.cartography import gridcell_angles -from iris.tests.stock import sample_2d_latlons +from iris.tests.stock import sample_2d_latlons, global_pp def _2d_multicells_testcube(cellsize_degrees=1.0): @@ -156,14 +156,6 @@ def test_coords_radians_args(self): self.assertArrayAllClose(result.data, self.standard_small_cube_results) - 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_bounds_array_args(self): # Check we can calculate from bounds values alone. co_x, co_y = (self.standard_regional_cube.coord(axis=ax) @@ -224,6 +216,89 @@ def test_unbounded_global(self): # 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) + + 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 = global_pp() + 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() From 605ef0cecc71ec3898788d50ae78c397978693e2 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 16 Aug 2018 16:14:27 +0100 Subject: [PATCH 13/15] Codestyle fixes. --- lib/iris/analysis/_grid_angles.py | 1 + lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py | 1 - 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py index a9d31412c5..2a172305ac 100644 --- a/lib/iris/analysis/_grid_angles.py +++ b/lib/iris/analysis/_grid_angles.py @@ -254,6 +254,7 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): # 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 diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py index 289ffab8c1..aa9b59c1b1 100644 --- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py +++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py @@ -299,6 +299,5 @@ def test_fail_bad_method(self): method='something_unknown') - if __name__ == "__main__": tests.main() From 1d373678fb7266c025decb3b04c214a691c62920 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 28 Aug 2018 17:24:12 +0100 Subject: [PATCH 14/15] Review changes + fixes. --- lib/iris/analysis/_grid_angles.py | 38 ++++++++++++++----- .../cartography/test_gridcell_angles.py | 24 ++++++++++++ 2 files changed, 52 insertions(+), 10 deletions(-) diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py index 2a172305ac..90876840af 100644 --- a/lib/iris/analysis/_grid_angles.py +++ b/lib/iris/analysis/_grid_angles.py @@ -78,8 +78,8 @@ def _latlon_from_xyz(xyz): """ lons = np.rad2deg(np.arctan2(xyz[1], xyz[0])) - axial_radii = np.sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]) - lats = np.rad2deg(np.arctan2(xyz[2], axial_radii)) + radii = np.sqrt(np.sum(xyz * xyz, axis=0)) + lats = np.rad2deg(np.arcsin(xyz[2] / radii)) return np.array([lons, lats]) @@ -102,6 +102,18 @@ def _angle(p, q, r): 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) @@ -191,11 +203,17 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): Cube of angles of grid-x vector from true Eastward direction for 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 - arrays, the output coords have bounds and centrepoints which are - the average of the 4 bounds. + 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 @@ -276,9 +294,9 @@ def transform_xy_arrays(x, y): x_coord = iris.coords.AuxCoord( points=xpts, bounds=xbds, standard_name='longitude', units='degrees') - 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 ? diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py index aa9b59c1b1..0a1ec0c03d 100644 --- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py +++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py @@ -28,6 +28,7 @@ import numpy as np +from cf_units import Unit from iris.cube import Cube from iris.coords import AuxCoord @@ -89,6 +90,7 @@ def setUp(self): # 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, @@ -125,6 +127,19 @@ def _check_multiple_orientations_and_latitudes(self, 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. @@ -222,6 +237,15 @@ def test_nonlatlon_coord_system(self): 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. From 5a047a637e750065ca06b8eddc2496dcabee7edb Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 28 Aug 2018 18:51:03 +0100 Subject: [PATCH 15/15] Avoid using sample data. --- .../tests/unit/analysis/cartography/test_gridcell_angles.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py index 0a1ec0c03d..4a1344e83d 100644 --- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py +++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py @@ -33,7 +33,7 @@ from iris.coords import AuxCoord from iris.analysis.cartography import gridcell_angles -from iris.tests.stock import sample_2d_latlons, global_pp +from iris.tests.stock import sample_2d_latlons, lat_lon_cube def _2d_multicells_testcube(cellsize_degrees=1.0): @@ -264,7 +264,7 @@ def test_fail_nonarraylike(self): def test_fail_non2d_coords(self): # Check error with bad args. - cube = global_pp() + cube = lat_lon_cube() with self.assertRaisesRegexp(ValueError, 'inputs must have 2-dimensional shape'): gridcell_angles(cube)