From 22dd241c1baff8ac6f6abb144d688b73f3b2ec55 Mon Sep 17 00:00:00 2001 From: Peter Killick Date: Fri, 13 Jul 2018 12:27:02 +0100 Subject: [PATCH 01/40] Quick fix to get tests passing by pinning dask version (#3086) --- requirements/core.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements/core.txt b/requirements/core.txt index c28cdb04a1..23a42dafe8 100644 --- a/requirements/core.txt +++ b/requirements/core.txt @@ -6,7 +6,7 @@ cartopy cf_units>=2 cftime -dask[array]>=0.17.1 #conda: dask>=0.17.1 +dask[array]>=0.17.1,<0.18.1 #conda: dask>=0.17.1,<0.18.1 matplotlib>=2 netcdf4 numpy>=1.14 From fcfea63277778fda98312127a07c5161190cdfc7 Mon Sep 17 00:00:00 2001 From: Carwyn Pelley Date: Fri, 13 Jul 2018 13:09:51 +0100 Subject: [PATCH 02/40] ENH: Working zonal mean linear regridding for circular sources or with use of extrapolation (#3085) * Working zonal mean from circular/extrapolated source to target * MAINT: Refactor of zonal mean testing * MAINT: Documentation changes from review --- lib/iris/analysis/_regrid.py | 5 +- lib/iris/tests/integration/test_regridding.py | 111 ++++++++++++++++++ 2 files changed, 114 insertions(+), 2 deletions(-) diff --git a/lib/iris/analysis/_regrid.py b/lib/iris/analysis/_regrid.py index eb7a16d075..71304a5dde 100644 --- a/lib/iris/analysis/_regrid.py +++ b/lib/iris/analysis/_regrid.py @@ -231,8 +231,9 @@ def _regrid(src_data, x_dim, y_dim, data = np.empty(shape, dtype=dtype) # The interpolation class requires monotonically increasing - # coordinates, so flip the coordinate(s) and data if the aren't. - reverse_x = src_x_coord.points[0] > src_x_coord.points[1] + # coordinates, so flip the coordinate(s) and data if they aren't. + reverse_x = (src_x_coord.points[0] > src_x_coord.points[1] if + src_x_coord.points.size > 1 else False) reverse_y = src_y_coord.points[0] > src_y_coord.points[1] flip_index = [slice(None)] * src_data.ndim if reverse_x: diff --git a/lib/iris/tests/integration/test_regridding.py b/lib/iris/tests/integration/test_regridding.py index ef521aa395..47d155f3f2 100644 --- a/lib/iris/tests/integration/test_regridding.py +++ b/lib/iris/tests/integration/test_regridding.py @@ -109,5 +109,116 @@ def test_nearest(self): self.assertArrayShapeStats(res, (1, 6, 3, 4), 315.890808, 11.000724) +class TestZonalMean_global(tests.IrisTest): + def setUp(self): + np.random.seed(0) + self.src = iris.cube.Cube(np.random.random_integers(0, 10, (140, 1))) + s_crs = iris.coord_systems.GeogCS(6371229.0) + sy_coord = iris.coords.DimCoord( + np.linspace(-90, 90, 140), standard_name='latitude', + units='degrees', coord_system=s_crs) + sx_coord = iris.coords.DimCoord( + -180, bounds=[-180, 180], standard_name='longitude', + units='degrees', circular=True, coord_system=s_crs) + self.src.add_dim_coord(sy_coord, 0) + self.src.add_dim_coord(sx_coord, 1) + + def test_linear_same_crs_global(self): + # Regrid the zonal mean onto an identical coordinate system target, but + # on a different set of longitudes - which should result in no change. + points = [-150, -90, -30, 30, 90, 150] + bounds = [[-180, -120], [-120, -60], [-60, 0], [0, 60], [60, 120], + [120, 180]] + sx_coord = self.src.coord(axis='x') + sy_coord = self.src.coord(axis='y') + x_coord = sx_coord.copy(points, bounds=bounds) + grid = iris.cube.Cube( + np.zeros([sy_coord.points.size, x_coord.points.size])) + grid.add_dim_coord(sy_coord, 0) + grid.add_dim_coord(x_coord, 1) + + res = self.src.regrid(grid, iris.analysis.Linear()) + + # Ensure data remains unchanged. + # (the same along each column) + self.assertTrue( + np.array([(res.data[:, 0]-res.data[:, i]).max() for i in + range(1, res.shape[1])]).max() < 1e-10) + self.assertArrayAlmostEqual(res.data[:, 0], self.src.data.reshape(-1)) + + +class TestZonalMean_regional(TestZonalMean_global, tests.IrisTest): + def setUp(self): + super(TestZonalMean_regional, self).setUp() + + # Define a target grid and a target result (what we expect the + # regridder to return). + sx_coord = self.src.coord(axis='x') + sy_coord = self.src.coord(axis='y') + grid_crs = iris.coord_systems.RotatedGeogCS( + 37.5, 177.5, ellipsoid=iris.coord_systems.GeogCS(6371229.0)) + grid_x = sx_coord.copy(np.linspace(350, 370, 100)) + grid_x.circular = False + grid_x.coord_system = grid_crs + grid_y = sy_coord.copy(np.linspace(-10, 10, 100)) + grid_y.coord_system = grid_crs + grid = iris.cube.Cube( + np.zeros([grid_y.points.size, grid_x.points.size])) + grid.add_dim_coord(grid_y, 0) + grid.add_dim_coord(grid_x, 1) + + # The target result is derived by regridding a multi-column version of + # the source to the target (i.e. turning a zonal mean regrid into a + # conventional regrid). + self.tar = self.zonal_mean_as_multi_column(self.src).regrid( + grid, iris.analysis.Linear()) + self.grid = grid + + def zonal_mean_as_multi_column(self, src_cube): + # Munge the source (duplicate source latitudes) so that we can + # utilise linear regridding as a conventional problem (that is, to + # duplicate columns so that it is no longer a zonal mean problem). + src_cube2 = src_cube.copy() + src_cube2.coord(axis='x').points = -90 + src_cube2.coord(axis='x').bounds = [-180, 0] + src_cube.coord(axis='x').points = 90 + src_cube.coord(axis='x').bounds = [0, 180] + src_cubes = iris.cube.CubeList([src_cube, src_cube2]) + return src_cubes.concatenate_cube() + + def test_linear_rotated_regional(self): + # Ensure that zonal mean source data is linearly interpolated onto a + # high resolution target. + regridder = iris.analysis.Linear() + res = self.src.regrid(self.grid, regridder) + self.assertArrayAlmostEqual(res.data, self.tar.data) + + def test_linear_rotated_regional_no_extrapolation(self): + # Capture the case where our source remains circular but we don't use + # extrapolation. + regridder = iris.analysis.Linear(extrapolation_mode='nan') + res = self.src.regrid(self.grid, regridder) + self.assertArrayAlmostEqual(res.data, self.tar.data) + + def test_linear_rotated_regional_not_circular(self): + # Capture the case where our source is not circular but we utilise + # extrapolation. + regridder = iris.analysis.Linear() + self.src.coord(axis='x').circular = False + res = self.src.regrid(self.grid, regridder) + self.assertArrayAlmostEqual(res.data, self.tar.data) + + def test_linear_rotated_regional_no_extrapolation_not_circular(self): + # Confirm how zonal mean actually works in so far as, that + # extrapolation and circular source handling is the means by which + # these usecases are supported. + # In the case where the source is neither using extrapolation and is + # not circular, then 'nan' values will result (as we would expect). + regridder = iris.analysis.Linear(extrapolation_mode='nan') + self.src.coord(axis='x').circular = False + res = self.src.regrid(self.grid, regridder) + self.assertTrue(np.isnan(res.data).all()) + + if __name__ == "__main__": tests.main() From 24c03299a8e401cae2dfa7ae8250c6d3a29fed5d Mon Sep 17 00:00:00 2001 From: hdyson Date: Fri, 13 Jul 2018 14:12:50 +0100 Subject: [PATCH 03/40] Avoid pandas deprecation warning. (#3079) Avoid pandas deprecation warning. This also simplifies code by being explicit about usage of "base" and "values" for numpy and pandas objects, and removes legacy code that may no longer be valid. --- lib/iris/pandas.py | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/lib/iris/pandas.py b/lib/iris/pandas.py index b4e1e9787a..244a7c011b 100644 --- a/lib/iris/pandas.py +++ b/lib/iris/pandas.py @@ -130,18 +130,7 @@ def _as_pandas_coord(coord): def _assert_shared(np_obj, pandas_obj): """Ensure the pandas object shares memory.""" - if hasattr(pandas_obj, 'base'): - base = pandas_obj.base - else: - base = pandas_obj[0].base - - # Prior to Pandas 0.17, when pandas_obj is a Series, pandas_obj.values - # returns a view of the underlying array, and pandas_obj.base, which calls - # pandas_obj.values.base, returns the underlying array. In 0.17 and 0.18 - # pandas_obj.values returns the underlying array, so base may be None even - # if the array is shared. - if base is None: - base = pandas_obj.values + values = pandas_obj.values def _get_base(array): # Chase the stack of NumPy `base` references back to the original array @@ -149,7 +138,7 @@ def _get_base(array): array = array.base return array - base = _get_base(base) + base = _get_base(values) np_base = _get_base(np_obj) if base is not np_base: msg = ('Pandas {} does not share memory' From 163e95854321e5ee82767693d17bd9e6e6adc869 Mon Sep 17 00:00:00 2001 From: Peter Killick Date: Mon, 16 Jul 2018 16:09:19 +0100 Subject: [PATCH 04/40] Workaround for dask array copy bug (#3088) --- lib/iris/tests/unit/cube/test_Cube.py | 4 +++- requirements/core.txt | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/lib/iris/tests/unit/cube/test_Cube.py b/lib/iris/tests/unit/cube/test_Cube.py index 9c96cf6b48..5633417c80 100644 --- a/lib/iris/tests/unit/cube/test_Cube.py +++ b/lib/iris/tests/unit/cube/test_Cube.py @@ -1430,7 +1430,9 @@ def test__masked_scalar_arraymask(self): self._check_copy(cube, cube.copy()) def test__lazy(self): - cube = Cube(as_lazy_data(np.array([1, 0]))) + # Note: multiple chunks added as a workaround suggested to dask#3751, + # which is fixed in dask#3754. + cube = Cube(as_lazy_data(np.array([1, 0]), chunks=(1, 1))) self._check_copy(cube, cube.copy()) diff --git a/requirements/core.txt b/requirements/core.txt index 23a42dafe8..c28cdb04a1 100644 --- a/requirements/core.txt +++ b/requirements/core.txt @@ -6,7 +6,7 @@ cartopy cf_units>=2 cftime -dask[array]>=0.17.1,<0.18.1 #conda: dask>=0.17.1,<0.18.1 +dask[array]>=0.17.1 #conda: dask>=0.17.1 matplotlib>=2 netcdf4 numpy>=1.14 From a0d18011b8b54f592fa6d417838aa0d79392f88f Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 2 Aug 2018 16:22:07 +0100 Subject: [PATCH 05/40] Pin Dask for avoid 0.18.2 bug with masked arrays. --- requirements/core.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements/core.txt b/requirements/core.txt index c28cdb04a1..997c91c665 100644 --- a/requirements/core.txt +++ b/requirements/core.txt @@ -6,7 +6,7 @@ cartopy cf_units>=2 cftime -dask[array]>=0.17.1 #conda: dask>=0.17.1 +dask[array]==0.18.1 #conda: dask==0.18.1 matplotlib>=2 netcdf4 numpy>=1.14 From c89f43e8a741b7126df2e92809062cca62eec146 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Fri, 3 Aug 2018 11:23:13 +0100 Subject: [PATCH 06/40] Tiny fix for dask, Python3 only? --- lib/iris/analysis/geometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/iris/analysis/geometry.py b/lib/iris/analysis/geometry.py index a415e10091..cca5d836ec 100644 --- a/lib/iris/analysis/geometry.py +++ b/lib/iris/analysis/geometry.py @@ -1,4 +1,4 @@ -# (C) British Crown Copyright 2010 - 2015, Met Office +# (C) British Crown Copyright 2010 - 2018, Met Office # # This file is part of Iris. # @@ -223,7 +223,7 @@ def geometry_area_weights(cube, geometry, normalize=False): else: slices.append(slice(None)) - weights[slices] = subweights + weights[tuple(slices)] = subweights # Fix for the limitation of iris.analysis.MEAN weights handling. # Broadcast the array to the full shape of the cube From a6de324d7240a3f9bca5226fe92526865f832c71 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Mon, 13 Aug 2018 16:46:34 +0100 Subject: [PATCH 07/40] Disable doctests for Python 2. --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index b92ecdbfb2..bd4c64b6f1 100644 --- a/.travis.yml +++ b/.travis.yml @@ -128,7 +128,7 @@ script: - INSTALL_DIR=$(pwd) - > - if [[ $TEST_TARGET == 'doctest' ]]; then + if [[ $TEST_TARGET == 'doctest' && ${TRAVIS_PYTHON_VERSION} == 3* ]]; then MPL_RC_DIR=$HOME/.config/matplotlib; mkdir -p $MPL_RC_DIR; echo 'backend : agg' > $MPL_RC_DIR/matplotlibrc; From 6e32e14263df7167f0fe5287d917754c88be2c8b Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Mon, 13 Aug 2018 16:46:56 +0100 Subject: [PATCH 08/40] Ignore warnings and update array printouts. --- .../src/userguide/interpolation_and_regridding.rst | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/iris/src/userguide/interpolation_and_regridding.rst b/docs/iris/src/userguide/interpolation_and_regridding.rst index 4f2b06de44..e3cd622541 100644 --- a/docs/iris/src/userguide/interpolation_and_regridding.rst +++ b/docs/iris/src/userguide/interpolation_and_regridding.rst @@ -5,6 +5,8 @@ import numpy as np import iris + import warnings + warnings.simplefilter('ignore') ================================= Cube interpolation and regridding @@ -137,10 +139,9 @@ This cube has a "hybrid-height" vertical coordinate system, meaning that the ver coordinate is unevenly spaced in altitude: >>> print(column.coord('altitude').points) - [418.6983642578125 434.57049560546875 456.79278564453125 485.3664855957031 - 520.2932739257812 561.5751953125 609.2144775390625 663.214111328125 - 723.5769653320312 790.306640625 863.4072265625 942.88232421875 - 1028.737060546875 1120.9764404296875 1219.6051025390625] + [ 418.69836 434.5705 456.7928 485.3665 520.2933 561.5752 + 609.2145 663.2141 723.57697 790.30664 863.4072 942.8823 + 1028.737 1120.9764 1219.6051 ] We could regularise the vertical coordinate by defining 10 equally spaced altitude sample points between 400 and 1250 and interpolating our vertical coordinate onto @@ -184,9 +185,8 @@ For example, to mask values that lie beyond the range of the original data: >>> scheme = iris.analysis.Linear(extrapolation_mode='mask') >>> new_column = column.interpolate(sample_points, scheme) >>> print(new_column.coord('altitude').points) - [nan 494.44451904296875 588.888916015625 683.333251953125 777.77783203125 - 872.2222290039062 966.666748046875 1061.111083984375 1155.555419921875 - nan] + [ nan 494.44452 588.8889 683.33325 777.77783 872.2222 + 966.66675 1061.1111 1155.5554 nan] .. _caching_an_interpolator: From d69222a795e147462d2f0044c946622703d4769f Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Mon, 13 Aug 2018 17:53:04 +0100 Subject: [PATCH 09/40] Better way to disable Python2 doctests. --- .travis.yml | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index bd4c64b6f1..ac5d341307 100644 --- a/.travis.yml +++ b/.travis.yml @@ -31,6 +31,15 @@ install: export IRIS_TEST_DATA_REF="2f3a6bcf25f81bd152b3d66223394074c9069a96"; export IRIS_TEST_DATA_SUFFIX=$(echo "${IRIS_TEST_DATA_REF}" | sed "s/^v//"); + # Cut short doctest phase under Python 2 : now only supports Python 3 + # SEE : https://github.com/SciTools/iris/pull/3134 + # ------------ + - > + if [[ $TEST_TARGET == 'doctest' && ${TRAVIS_PYTHON_VERSION} != 3* ]]; then + echo "DOCTEST phase only valid in Python 3 : ABORTING during 'install'." + exit 0 + fi + # Install miniconda # ----------------- - > @@ -128,7 +137,7 @@ script: - INSTALL_DIR=$(pwd) - > - if [[ $TEST_TARGET == 'doctest' && ${TRAVIS_PYTHON_VERSION} == 3* ]]; then + if [[ $TEST_TARGET == 'doctest' ]]; then MPL_RC_DIR=$HOME/.config/matplotlib; mkdir -p $MPL_RC_DIR; echo 'backend : agg' > $MPL_RC_DIR/matplotlibrc; From ff7ec92cd1294d188bebbe034b1e6b22c57e0ca6 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Mon, 23 Jul 2018 16:10:40 +0100 Subject: [PATCH 10/40] Add gridcell_angles and rotate_grid_vectors to iris.analysis.cartography, with tests: INCOMPLETE WIP --- lib/iris/analysis/_grid_angles.py | 261 ++++++++++++++++++ lib/iris/analysis/cartography.py | 3 + .../cartography/test_gridcell_angles.py | 171 ++++++++++++ 3 files changed, 435 insertions(+) create mode 100644 lib/iris/analysis/_grid_angles.py 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 new file mode 100644 index 0000000000..6f7dcf15db --- /dev/null +++ b/lib/iris/analysis/_grid_angles.py @@ -0,0 +1,261 @@ +import numpy as np + +import iris + + +def _3d_xyz_from_latlon(lon, lat): + """ + Return locations of (lon, lat) in 3D space. + + Args: + + * lon, lat: (arrays in degrees) + + Returns: + + x, y, z : (array, dtype=float64) + cartesian coordinates on a unit sphere. + + """ + lon1 = np.deg2rad(lon).astype(np.float64) + lat1 = np.deg2rad(lat).astype(np.float64) + + old_way = True + if old_way: + x3 = np.empty((3,) + lon.shape, dtype=float) + x3[0] = np.cos(lat1) * np.cos(lon1) + x3[1] = np.cos(lat1) * np.sin(lon1) + x3[2] = np.sin(lat1) + result = x3 + else: + x = np.cos(lat1) * np.cos(lon1) + y = np.cos(lat1) * np.sin(lon1) + z = np.sin(lat1) + + result = np.concatenate([array[np.newaxis] for array in (x, y, z)]) + + return result + + +def _latlon_from_xyz(xyz): + """ + Return arrays of lons+lats angles from xyz locations. + + Args: + + * xyz: (array) + positions array, of dims (3, ), where index 0 maps x/y/z. + + Returns: + + lonlat : (array) + spherical angles, of dims (2, ), in radians. + Dim 0 maps longitude, latitude. + + """ + lons = 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) + return np.array([lons, lats]) + + +def _angle(p, q, r): + """ + Return angle (in _radians_) of grid wrt local east. Anticlockwise +ve, as usual. + {P, Q, R} are consecutve 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. + + """ + q1 = np.deg2rad(q) + + 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(q1[0]) - pr[0] * np.sin(q1[0]) + + index = np.where(pr_norm == 0) + pr_norm[index] = 1 + + cosine = np.maximum(np.minimum(pr_top / pr_norm, 1), -1) + cosine[index] = 0 + + psi = np.arccos(cosine) * np.sign(r[1] - p[1]) + psi[index] = np.nan + + return psi + + +def gridcell_angles(x, y=None): + """ + Calculate gridcell orientation angles. + + The inputs (x [,y]) can be any of the folliwing : + + * x (:class:`~iris.cube.Cube`): + a grid cube with 2D longitude and latitude coordinates. + + * x, y (:class:`~iris.coords.Coord`): + longitude and latitude coordinates. + + * x, y (2-dimensional arrays of same shape (ny, nx)): + longitude and latitude cell center locations, in degrees. + + * x, y (3-dimensional arrays of same shape (ny, nx, 4)): + longitude and latitude cell bounds, in degrees. + The last index maps cell corners anticlockwise from bottom-left. + + Returns: + + angles : (2-dimensional cube) + + Cube of angles of grid-x vector from true Eastward direction for + each gridcell, in radians. + 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. + + """ + if hasattr(x, 'core_data'): + # N.B. only "true" lats + longs will do : Cannot handle rotated ! + x, y = x.coord('longitude'), x.coord('latitude') + + # Now should have either 2 coords or 2 arrays. + if not hasattr(x, 'shape') and hasattr(y, 'shape'): + msg = ('Inputs (x,y) must have array shape property.' + 'Got type(x)={} and type(y)={}.') + raise ValueError(msg.format(type(x), type(y))) + + x_coord, y_coord = None, None + if isinstance(x, iris.coords.Coord) and isinstance(y, iris.coords.Coord): + 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 {}.') + raise ValueError(msg.format(x_coord.shape, y_coord.shape)) + if x_coord.shape != y_coord.shape: + msg = ('Coordinate inputs must have same shape. ', + 'Got x-shape of {} and y-shape of {}.') + raise ValueError(msg.format(x_coord.shape, y_coord.shape)) +# 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 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 + + elif isinstance(x, iris.coords.Coord) or isinstance(y, iris.coords.Coord): + is_and_not = ('x', 'y') + if isinstance(y, iris.coords.Coord): + is_and_not = reversed(is_and_not) + 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. + 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 ? + + # 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. + mid = np.array([x, y]) + lhs = np.roll(mid, 1, 2) + rhs = np.roll(mid, -1, 2) + if not x_coord: + # Create coords for result cube : with no bounds. + y_coord = iris.coords.AuxCoord(x, standard_name='latitude', + units='degrees') + x_coord = iris.coords.AuxCoord(y, standard_name='longitude', + units='degrees') + else: + # Get lhs and rhs locations by averaging top+bottom each side. + # NOTE: so with bounds, we *don't* need full circular longitudes. + xyz = _3d_xyz_from_latlon(x, y) + lhs_xyz = 0.5 * (xyz[..., 0] + xyz[..., 3]) + rhs_xyz = 0.5 * (xyz[..., 1] + xyz[..., 2]) + 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. + mid_xyz = 0.5 * (lhs_xyz + rhs_xyz) + mid_latlons = _latlon_from_xyz(mid_xyz) + # Create coords with given bounds, and averaged centrepoints. + x_coord = iris.coords.AuxCoord( + points=mid_latlons[0], bounds=x, + standard_name='longitude', units='degrees') + y_coord = iris.coords.AuxCoord( + points=mid_latlons[1], bounds=y, + standard_name='latitude', units='degrees') + # Convert lhs and rhs points back to latlon form -- IN DEGREES ! + lhs = 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 = 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') + 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): + """ + Rotate distance vectors from grid-oriented to true-latlon-oriented. + + Args: + + * u_cube, v_cube : (cube) + Cubes of grid-u and grid-v vector components. + Units should be differentials of true-distance, e.g. 'm/s'. + + Optional args: + + * grid_angles_cube : (cube) + gridcell orientation angles. + Units must be angular, i.e. can be converted to 'radians'. + If not provided, grid angles are estimated from 'u_cube' using the + :func:`gridcell_angles` method. + + Returns: + + true_u, true_v : (cube) + Cubes of true-north oriented vector components. + Units are same as inputs. + + """ + u_out, v_out = (cube.copy() for cube in (u_cube, v_cube)) + if not grid_angles_cube: + grid_angles_cube = gridcell_angles(u_cube) + gridangles = grid_angles_cube.copy() + gridangles.convert_units('radians') + uu, vv, aa = (cube.data for cube in (u_out, v_out, gridangles)) + mags = np.sqrt(uu*uu + vv*vv) + angs = np.arctan2(vv, uu) + aa + uu, vv = mags * np.cos(angs), mags * np.sin(angs) + + # Promote all to masked arrays, and also apply mask at bad (NaN) angles. + mask = np.isnan(aa) + for cube in (u_out, v_out, aa): + if hasattr(cube.data, 'mask'): + mask |= cube.data.mask + u_out.data = np.ma.masked_array(uu, mask=mask) + v_out.data = np.ma.masked_array(vv, mask=mask) + + return u_out, v_out diff --git a/lib/iris/analysis/cartography.py b/lib/iris/analysis/cartography.py index 4b68dcc949..1c55d48fd3 100644 --- a/lib/iris/analysis/cartography.py +++ b/lib/iris/analysis/cartography.py @@ -37,6 +37,9 @@ 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) # This value is used as a fall-back if the cube does not define the earth 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..b3a10a783f --- /dev/null +++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py @@ -0,0 +1,171 @@ +# (C) British Crown Copyright 2015 - 2016, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +""" +Unit tests for the 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) + + +def _rotated_grid_sample(pole_lat=15, pole_lon=-180, + lon_bounds=np.linspace(-30, 30, 8, endpoint=True), + lat_bounds=np.linspace(-30, 30, 8, 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 test_values(self): + # Construct a rotated-pole grid and check angle calculation. + testcube = _rotated_grid_sample() + angles_cube = gridcell_angles(testcube) + angles_cube.convert_units('radians') + + # testing phase... + print(np.rad2deg(angles_cube.data)) + + import matplotlib.pyplot as plt + plt.switch_backend('tkagg') + ax = plt.axes(projection=ccrs.Orthographic(central_longitude=0.0, + central_latitude=90.0,)) + 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=ccrs.PlateCarree()) + + + # Show plain 0,1 + 1,0 vectors unrotated at the given points. + pts_shape = testcube.coord('longitude').shape + u0 = np.ones(pts_shape) + v0 = np.zeros(pts_shape) + u1 = v0.copy() + v1 = u0.copy() + +# 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) + + scale = 4.0e-6 + plt.quiver(testcube.coord('longitude').points, + testcube.coord('latitude').points, + u0, v0, color='blue', linewidth=0.5, scale_units='xy', scale=scale, + transform=ccrs.PlateCarree()) + plt.quiver(testcube.coord('longitude').points, + testcube.coord('latitude').points, + u1, v1, color='red', linewidth=0.5, scale_units='xy', scale=scale, + transform=ccrs.PlateCarree()) + +# plt.quiver(testcube.coord('longitude').points, +# testcube.coord('latitude').points, +# u0r_cube.data, v0r_cube.data, +# color='red', +# transform=ccrs.PlateCarree()) + + # Also draw small lines pointing at the correct angle. + x0s = testcube.coord('longitude').points + y0s = testcube.coord('latitude').points + ny, nx = x0s.shape + size_degrees = 5.0 + angles = angles_cube.copy() + angles.convert_units('radians') + angles = angles.data + lats = testcube.coord('latitude').copy() + lats.convert_units('radians') + lats = lats.points + dys = size_degrees * np.sin(angles) / np.cos(-lats) + dxs = size_degrees * np.cos(angles) + 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', + transform=ccrs.PlateCarree()) + + + + ax.set_global() + plt.show() + + + self.assertArrayAllClose( + angles_cube.data, + [[100.0, 100.0, 100.0], + [100.0, 100.0, 100.0], + [100.0, 100.0, 100.0]]) + + +if __name__ == "__main__": + tests.main() From 393489eb6f827a4675fb9126bc2d066254b28371 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 25 Jul 2018 20:22:48 +0100 Subject: [PATCH 11/40] Roughly working, snapshotted with complex test plot code, to be reduced. --- lib/iris/analysis/_grid_angles.py | 136 ++++++-- lib/iris/analysis/cartography.py | 2 +- .../cartography/test_gridcell_angles.py | 323 +++++++++++++++--- 3 files changed, 375 insertions(+), 86 deletions(-) diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py index 6f7dcf15db..f5fa94a256 100644 --- a/lib/iris/analysis/_grid_angles.py +++ b/lib/iris/analysis/_grid_angles.py @@ -13,26 +13,18 @@ def _3d_xyz_from_latlon(lon, lat): Returns: - x, y, z : (array, dtype=float64) - cartesian coordinates on a unit sphere. + xyz : (array, dtype=float64) + cartesian coordinates on a unit sphere. Dimension 0 maps x,y,z. """ lon1 = np.deg2rad(lon).astype(np.float64) lat1 = np.deg2rad(lat).astype(np.float64) - old_way = True - if old_way: - x3 = np.empty((3,) + lon.shape, dtype=float) - x3[0] = np.cos(lat1) * np.cos(lon1) - x3[1] = np.cos(lat1) * np.sin(lon1) - x3[2] = np.sin(lat1) - result = x3 - else: - x = np.cos(lat1) * np.cos(lon1) - y = np.cos(lat1) * np.sin(lon1) - z = np.sin(lat1) + x = np.cos(lat1) * np.cos(lon1) + y = np.cos(lat1) * np.sin(lon1) + z = np.sin(lat1) - result = np.concatenate([array[np.newaxis] for array in (x, y, z)]) + result = np.concatenate([array[np.newaxis] for array in (x, y, z)]) return result @@ -61,35 +53,82 @@ 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 consecutve 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. + 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. """ - q1 = np.deg2rad(q) +# old_style = True + old_style = False + if old_style: + mid_lons = np.deg2rad(q[0]) - pr = _3d_xyz_from_latlon(r[0], r[1]) - _3d_xyz_from_latlon(p[0], p[1]) - pr_norm = np.sqrt(np.sum(pr**2, axis=0)) - pr_top = pr[1] * np.cos(q1[0]) - pr[0] * np.sin(q1[0]) + pr = _3d_xyz_from_latlon(r[0], r[1]) - _3d_xyz_from_latlon(p[0], p[1]) + pr_norm = np.sqrt(np.sum(pr**2, axis=0)) + pr_top = pr[1] * np.cos(mid_lons) - pr[0] * np.sin(mid_lons) - index = np.where(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 + 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 -def gridcell_angles(x, y=None): +def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): """ Calculate gridcell orientation angles. + Args: + The inputs (x [,y]) can be any of the folliwing : * x (:class:`~iris.cube.Cube`): @@ -105,6 +144,16 @@ def gridcell_angles(x, y=None): longitude and latitude cell bounds, in degrees. The last index maps cell corners anticlockwise from bottom-left. + Optional Args: + + * cell_angle_boundpoints (string): + Controls which gridcell bounds locations are used to calculate angles, + if the inputs are bounds or bounded coordinates. + Valid values are 'lower-left, lower-right', which takes the angle from + the lower left to the lower right corner, and 'mid-lhs, mid-rhs' which + takes an angles between the average of the left-hand and right-hand + pairs of corners. The default is 'mid-lhs, mid-rhs'. + Returns: angles : (2-dimensional cube) @@ -184,8 +233,20 @@ def gridcell_angles(x, y=None): # Get lhs and rhs locations by averaging top+bottom each side. # NOTE: so with bounds, we *don't* need full circular longitudes. xyz = _3d_xyz_from_latlon(x, y) - lhs_xyz = 0.5 * (xyz[..., 0] + xyz[..., 3]) - rhs_xyz = 0.5 * (xyz[..., 1] + xyz[..., 2]) + angle_boundpoints_vals = {'mid-lhs, mid-rhs': '03_to_12', + 'lower-left, lower-right': '0_to_1'} + bounds_pos = angle_boundpoints_vals.get(cell_angle_boundpoints) + if bounds_pos == '0_to_1': + lhs_xyz = xyz[..., 0] + rhs_xyz = xyz[..., 1] + elif bounds_pos == '03_to_12': + lhs_xyz = 0.5 * (xyz[..., 0] + xyz[..., 3]) + rhs_xyz = 0.5 * (xyz[..., 1] + xyz[..., 2]) + else: + msg = ('unrecognised cell_angle_boundpoints of "{}", ' + 'must be one of {}') + raise ValueError(msg.format(cell_angle_boundpoints, + list(angle_boundpoints_vals.keys()))) if not x_coord: # Create bounded coords for result cube. # Use average lhs+rhs points in 3d to get 'mid' points, as coords @@ -215,7 +276,9 @@ def gridcell_angles(x, y=None): return result -def true_vectors_from_grid_vectors(u_cube, v_cube, grid_angles_cube=None): +def true_vectors_from_grid_vectors(u_cube, v_cube, + grid_angles_cube=None, + grid_angles_kwargs=None): """ Rotate distance vectors from grid-oriented to true-latlon-oriented. @@ -233,6 +296,10 @@ def true_vectors_from_grid_vectors(u_cube, v_cube, grid_angles_cube=None): If not provided, grid angles are estimated from 'u_cube' using the :func:`gridcell_angles` method. + * grid_angles_kwargs : (dict or None) + Additional keyword args to be passed to the :func:`gridcell_angles` + method, if it is used. + Returns: true_u, true_v : (cube) @@ -242,7 +309,8 @@ def true_vectors_from_grid_vectors(u_cube, v_cube, grid_angles_cube=None): """ u_out, v_out = (cube.copy() for cube in (u_cube, v_cube)) if not grid_angles_cube: - grid_angles_cube = gridcell_angles(u_cube) + grid_angles_kwargs = grid_angles_kwargs or {} + grid_angles_cube = gridcell_angles(u_cube, **grid_angles_kwargs) gridangles = grid_angles_cube.copy() gridangles.convert_units('radians') uu, vv, aa = (cube.data for cube in (u_out, v_out, gridangles)) diff --git a/lib/iris/analysis/cartography.py b/lib/iris/analysis/cartography.py index 1c55d48fd3..16c99c83a8 100644 --- a/lib/iris/analysis/cartography.py +++ b/lib/iris/analysis/cartography.py @@ -1,4 +1,4 @@ -# (C) British Crown Copyright 2010 - 2017, Met Office +# (C) British Crown Copyright 2010 - 2018, Met Office # # This file is part of Iris. # 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 b3a10a783f..a73988f668 100644 --- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py +++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py @@ -1,4 +1,4 @@ -# (C) British Crown Copyright 2015 - 2016, Met Office +# (C) British Crown Copyright 2018, Met Office # # This file is part of Iris. # @@ -38,10 +38,13 @@ 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, 8, endpoint=True), - lat_bounds=np.linspace(-30, 30, 8, endpoint=True)): + 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) @@ -79,19 +82,151 @@ def expand_unified_bds(bds): class TestGridcellAngles(tests.IrisTest): + def _singlecell_30deg_cube(self, x0=90., y0=0., dx=20., dy=10.): + x_pts = np.array([[x0]]) + y_pts = np.array([[y0]]) + x_bds = x0 + dx * np.array([[[-1., 1, 0.5, -1.5]]]) +# self.assertArrayAllClose(x_bds, np.array([[[70., 110, 100, 60]]])) + y_bds = y0 + dy * np.array([[[-1., 1, 3, 1]]]) +# self.assertArrayAllClose(y_bds, np.array([[[-10., 10, 30, 10]]])) + co_x = AuxCoord(points=x_pts, bounds=x_bds, + standard_name='longitude', units='degrees') + co_y = AuxCoord(points=y_pts, bounds=y_bds, + standard_name='latitude', units='degrees') + cube = Cube(np.zeros((1, 1))) + cube.add_aux_coord(co_x, (0, 1)) + cube.add_aux_coord(co_y, (0, 1)) + return cube + + def _singlecell_diamond_cube(self, x0=90., y0=0., dy=10., dx_eq=None): + if dx_eq is None: + dx_eq = dy + x_pts = np.array([[x0]]) + y_pts = np.array([[y0]]) + dx = dx_eq / np.cos(np.deg2rad(y0)) + x_bds = np.array([[[x0, x0 + dx, x0, x0 - dx]]]) + y_bds = np.array([[[y0 - dy, y0, y0 + dy, y0]]]) + co_x = AuxCoord(points=x_pts, bounds=x_bds, + standard_name='longitude', units='degrees') + co_y = AuxCoord(points=y_pts, bounds=y_bds, + standard_name='latitude', units='degrees') + cube = Cube(np.zeros((1, 1))) + cube.add_aux_coord(co_x, (0, 1)) + cube.add_aux_coord(co_y, (0, 1)) + return cube + + def test_single_cell_equatorial(self): + plt.switch_backend('tkagg') + plt.figure(figsize=(10,10)) + ax = plt.axes(projection=ccrs.Mercator()) +# ax = plt.axes(projection=ccrs.NorthPolarStereo()) +# ax = plt.axes(projection=ccrs.Orthographic(central_longitude=90., +# central_latitude=30.)) + + lon0 = 90.0 + dy = 2.0 + y_0, y_n, ny = -80, 80, 9 + angles = [] + for lat in np.linspace(y_0, y_n, ny): + cube = self._singlecell_diamond_cube(x0=lon0, y0=lat, dy=dy) + angles_cube = gridcell_angles(cube, +# cell_angle_boundpoints='mid-lhs, mid-rhs') + cell_angle_boundpoints='lower-left, lower-right') + tmp_cube = angles_cube.copy() + tmp_cube.convert_units('degrees') + print('') + print(lat) + co_x, co_y = (cube.coord(axis=ax) for ax in ('x', 'y')) + print() + print(' at : {}, {}'.format(co_x.points[0, 0], co_y.points[0, 0])) + print(' x-bds:') + print(co_x.bounds) + print(' y-bds:') + print(co_y.bounds) + angle = tmp_cube.data[0, 0] + angles.append(angle) + print(angle) + blockplot_2dll(cube) + + ax.coastlines() + ax.set_global() + + # Plot constant NEly (45deg) arrows. + xx = np.array([lon0] * ny) + yy = np.linspace(y_0, y_n, ny) - dy + uu = np.array([1.0] * ny) + plt.quiver(xx, yy, + uu, np.cos(np.deg2rad(yy)), + zorder=2, color='red', + scale_units='xy', + transform=ccrs.PlateCarree()) + + # Also plot returned angles. + angles_arr_rad = np.deg2rad(angles) + u_arr = uu * np.cos(angles_arr_rad) + v_arr = uu * np.sin(angles_arr_rad) * np.cos(np.deg2rad(yy)) + + plt.quiver(xx, yy, + u_arr, + v_arr, + zorder=2, color='magenta', + scale_units='xy', + width=0.005, + scale=0.2e-6, +# width=0.5, + transform=ccrs.PlateCarree()) + + plt.show() + + def test_values(self): # Construct a rotated-pole grid and check angle calculation. testcube = _rotated_grid_sample() - angles_cube = gridcell_angles(testcube) + + 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') - ax = plt.axes(projection=ccrs.Orthographic(central_longitude=0.0, - central_latitude=90.0,)) + +# 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): @@ -100,71 +235,157 @@ def test_values(self): testcube.coord('latitude').bounds[..., i_bnd], '+', markersize=10., markeredgewidth=2., markerfacecolor=color, markeredgecolor=color, - transform=ccrs.PlateCarree()) + transform=data_proj) - # Show plain 0,1 + 1,0 vectors unrotated at the given points. + # 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() -# 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) + 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) - scale = 4.0e-6 - plt.quiver(testcube.coord('longitude').points, - testcube.coord('latitude').points, - u0, v0, color='blue', linewidth=0.5, scale_units='xy', scale=scale, - transform=ccrs.PlateCarree()) - plt.quiver(testcube.coord('longitude').points, - testcube.coord('latitude').points, - u1, v1, color='red', linewidth=0.5, scale_units='xy', scale=scale, - transform=ccrs.PlateCarree()) + # 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) -# plt.quiver(testcube.coord('longitude').points, -# testcube.coord('latitude').points, -# u0r_cube.data, v0r_cube.data, -# color='red', -# transform=ccrs.PlateCarree()) - # Also draw small lines pointing at the correct angle. - x0s = testcube.coord('longitude').points - y0s = testcube.coord('latitude').points - ny, nx = x0s.shape - size_degrees = 5.0 - angles = angles_cube.copy() - angles.convert_units('radians') - angles = angles.data - lats = testcube.coord('latitude').copy() - lats.convert_units('radians') - lats = lats.points - dys = size_degrees * np.sin(angles) / np.cos(-lats) - dxs = size_degrees * np.cos(angles) - x1s = x0s + dxs - y1s = y0s + dys + + # + # 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], x1s[iy, ix]], - [y0s[iy, ix], y1s[iy, ix]], - 'o-', markersize=4., markeredgewidth=0., - color='green', - transform=ccrs.PlateCarree()) + 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() +# plt.show() + angles_cube.convert_units('degrees') self.assertArrayAllClose( angles_cube.data, - [[100.0, 100.0, 100.0], - [100.0, 100.0, 100.0], - [100.0, 100.0, 100.0]]) + [[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__": From d8bb331ca2774f34f538b9746e4cdc43e935dac8 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 26 Jul 2018 16:14:04 +0100 Subject: [PATCH 12/40] Small improvements. --- lib/iris/analysis/_grid_angles.py | 79 ++++++++++++++++--- .../cartography/test_gridcell_angles.py | 2 +- 2 files changed, 71 insertions(+), 10 deletions(-) diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py index f5fa94a256..1746834a65 100644 --- a/lib/iris/analysis/_grid_angles.py +++ b/lib/iris/analysis/_grid_angles.py @@ -1,3 +1,27 @@ +# (C) British Crown Copyright 2010 - 2018, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +""" +Code to implement vector rotation by angles, and inferring gridcell angles +from coordinate points and bounds. + +""" +from __future__ import (absolute_import, division, print_function) +from six.moves import (filter, input, map, range, zip) # noqa + import numpy as np import iris @@ -125,23 +149,41 @@ def _angle(p, q, r): def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): """ - Calculate gridcell orientation angles. + Calculate gridcell orientations for an arbitrary 2-dimensional grid. + + The input grid is defined by two 2-dimensional coordinate arrays with the + same dimensions (ny, nx), specifying the geolocations of a 2D mesh. + + Input values may be coordinate points (ny, nx) or bounds (ny, nx, 4). + However, if points, the edges in the X direction are assumed to be + connected by wraparound. + + Input can be either two arrays, two coordinates, or a single cube + containing two suitable coordinates identified with the 'x' and'y' axes. Args: The inputs (x [,y]) can be any of the folliwing : * x (:class:`~iris.cube.Cube`): - a grid cube with 2D longitude and latitude coordinates. + a grid cube with 2D X and Y coordinates, identified by 'axis'. + The coordinates must be 2-dimensional with the same shape. + The two dimensions represent grid dimensions in the order Y, then X. * x, y (:class:`~iris.coords.Coord`): - longitude and latitude coordinates. + X and Y coordinates, specifying grid locations on the globe. + The coordinates must be 2-dimensional with the same shape. + The two dimensions represent grid dimensions in the order Y, then X. + If there is no coordinate system, they are assumed to be true + longitudes and latitudes. Units must convertible to 'degrees'. * x, y (2-dimensional arrays of same shape (ny, nx)): longitude and latitude cell center locations, in degrees. + The two dimensions represent grid dimensions in the order Y, then X. * x, y (3-dimensional arrays of same shape (ny, nx, 4)): longitude and latitude cell bounds, in degrees. + The first two dimensions are grid dimensions in the order Y, then X. The last index maps cell corners anticlockwise from bottom-left. Optional Args: @@ -167,9 +209,11 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): the average of the 4 bounds. """ - if hasattr(x, 'core_data'): - # N.B. only "true" lats + longs will do : Cannot handle rotated ! - x, y = x.coord('longitude'), x.coord('latitude') + cube = None + if hasattr(x, 'add_aux_coord'): + # Passed a cube : extract 'x' and ;'y' axis coordinates. + cube = x # Save for later checking. + x, y = cube.coord(axis='x'), cube.coord(axis='y') # Now should have either 2 coords or 2 arrays. if not hasattr(x, 'shape') and hasattr(y, 'shape'): @@ -178,7 +222,7 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): raise ValueError(msg.format(type(x), type(y))) x_coord, y_coord = None, None - if isinstance(x, iris.coords.Coord) and isinstance(y, iris.coords.Coord): + if hasattr(x, 'bounds') and hasattr(y, 'bounds'): x_coord, y_coord = x.copy(), y.copy() x_coord.convert_units('degrees') y_coord.convert_units('degrees') @@ -203,9 +247,10 @@ def gridcell_angles(x, y=None, cell_angle_boundpoints='mid-lhs, mid-rhs'): else: x, y = x_coord.points, y_coord.points - elif isinstance(x, iris.coords.Coord) or isinstance(y, iris.coords.Coord): + elif hasattr(x, 'bounds') or hasattr(y, 'bounds'): + # One was a Coord, and the other not ? is_and_not = ('x', 'y') - if isinstance(y, iris.coords.Coord): + 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) @@ -282,6 +327,18 @@ def true_vectors_from_grid_vectors(u_cube, v_cube, """ Rotate distance vectors from grid-oriented to true-latlon-oriented. + .. Note:: + + This operation overlaps somewhat in function with + :func:`iris.analysis.cartography.rotate_winds`. + However, that routine only rotates vectors according to transformations + between coordinate systems. + This function, by contrast, can rotate vectors by arbitrary angles. + Most commonly, the angles are estimated solely from grid sampling + points, using :func:`gridcell_angles` : This allows operation on + complex meshes defined by two-dimensional coordinates, such as most + ocean grids. + Args: * u_cube, v_cube : (cube) @@ -306,6 +363,10 @@ def true_vectors_from_grid_vectors(u_cube, v_cube, Cubes of true-north oriented vector components. Units are same as inputs. + .. Note:: + + Vector magnitudes will always be the same as the inputs. + """ u_out, v_out = (cube.copy() for cube in (u_cube, v_cube)) if not grid_angles_cube: 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 a73988f668..de00d2bc76 100644 --- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py +++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py @@ -374,7 +374,7 @@ def test_values(self): ax.set_global() -# plt.show() + plt.show() angles_cube.convert_units('degrees') From c902ce3460279588e4d7cf0d9523372479f9d1d2 Mon Sep 17 00:00:00 2001 From: lbdreyer Date: Fri, 20 Jul 2018 10:19:20 +0100 Subject: [PATCH 13/40] Support plotting 2D bounded coords --- lib/iris/coords.py | 129 ++++++++++-- lib/iris/plot.py | 160 ++++++++++++--- lib/iris/tests/stock.py | 217 ++++++++++++--------- lib/iris/tests/unit/plot/test_2d_coords.py | 158 +++++++++++++++ 4 files changed, 529 insertions(+), 135 deletions(-) create mode 100644 lib/iris/tests/unit/plot/test_2d_coords.py diff --git a/lib/iris/coords.py b/lib/iris/coords.py index e4a8eb35f2..57ae723205 100644 --- a/lib/iris/coords.py +++ b/lib/iris/coords.py @@ -143,6 +143,85 @@ def __new__(cls, name_or_coord, minimum, maximum, 'groupby_point, groupby_slice') +def _discontinuity_in_2d_bounds(bds, abs_tol=1e-4): + """ + Check bounds of a 2-dimensional coordinate are contiguous + Args: + bds: Array of bounds of shape (X,Y,4) + abs_tol: tolerance + + Returns: + Bool, if there are no discontinuities + absolute difference along the x axis + absolute difference along the y axis + + """ + # Check form is (ny, nx, 4) + if not bds.ndim == 3 and bds.shape[2] == 4: + raise ValueError('2D coordinates must have 4 bounds per point ' + 'for 2D coordinate plotting') + + # Check ordering: + # i i+1 + # j @0 @1 + # j+1 @3 @2 + def mod360_diff(x1, x2): + diff = x1 - x2 + diff = (diff + 360.0 + 180.0) % 360.0 - 180.0 + return diff + + # Compare cell with the cell next to it (i+1) + diffs_along_x = mod360_diff(bds[:, :-1, 1], bds[:, 1:, 0]) + # Compare cell with the cell above it (j+1) + diffs_along_y = mod360_diff(bds[:-1, :, 3], bds[1:, :, 0]) + + def eq_diffs(x1): + return np.all(np.abs(x1) < abs_tol) + + match_y0_x1 = eq_diffs(diffs_along_x) + match_y1_x0 = eq_diffs(diffs_along_y) + + all_eq = match_y0_x1 and match_y1_x0 + + return all_eq, np.abs(diffs_along_x), np.abs(diffs_along_y) + + +def _get_2d_coord_bound_grid(bds): + """ + Function used that takes a bounds array for a 2-D coordinate variable with + 4 sides and returns the bounds grid. + + Cf standards requires the four vertices of the cell to be traversed + anti-clockwise if the coordinates are defined in a right handed coordinate + system. + + selects the zeroth vertex of each cell and then adds the column the first + vertex at the end. For the top row it uses the thirs vertex, with the + second added on to the end. + e.g. + # 0-0-0-0-1 + # 0-0-0-0-1 + # 3-3-3-3-2 + + + Args: + bounds: array of shape (X,Y,4) + + Returns: + array of shape (X+1, Y+1) + + """ + bds_shape = bds.shape + result = np.zeros((bds_shape[0] + 1, bds_shape[1] + 1)) + + result[:-1, :-1] = bds[:, :, 0] + result[:-1, -1] = bds[:, -1, 1] + result[-1, :-1] = bds[-1, :, 3] + result[-1, -1] = bds[-1, -1, 2] + + return result + + class Cell(collections.namedtuple('Cell', ['point', 'bound'])): """ An immutable representation of a single cell of a coordinate, including the @@ -951,15 +1030,22 @@ def cells(self): return _CellIterator(self) def _sanity_check_contiguous(self): - if self.ndim != 1: - raise iris.exceptions.CoordinateMultiDimError( - 'Invalid operation for {!r}. Contiguous bounds are not defined' - ' for multi-dimensional coordinates.'.format(self.name())) - if self.nbounds != 2: - raise ValueError( - 'Invalid operation for {!r}, with {} bounds. Contiguous bounds' - ' are only defined for coordinates with 2 bounds.'.format( - self.name(), self.nbounds)) + if self.ndim == 1: + if self.nbounds != 2: + raise ValueError('Invalid operation for {!r}, with {} bounds. ' + 'Contiguous bounds are only defined for 1D ' + 'coordinates with 2 bounds.'.format + (self.name(), self.nbounds)) + elif self.ndim == 2: + if self.nbounds != 4: + raise ValueError('Invalid operation for {!r}, with {} bounds. ' + 'Contiguous bounds are only defined for 2D ' + 'coordinates with 4 bounds.'.format + (self.name(), self.nbounds)) + else: + raise ValueError('Invalid operation for {!r}. Contiguous bounds ' + 'are not defined for coordinates with more than ' + '2 dimensions.'.format(self.name())) def is_contiguous(self, rtol=1e-05, atol=1e-08): """ @@ -979,19 +1065,29 @@ def is_contiguous(self, rtol=1e-05, atol=1e-08): """ if self.has_bounds(): self._sanity_check_contiguous() - return np.allclose(self.bounds[1:, 0], self.bounds[:-1, 1], - rtol=rtol, atol=atol) + if self.ndim == 1: + return np.allclose(self.bounds[1:, 0], self.bounds[:-1, 1], + rtol=rtol, atol=atol) + elif self.ndim == 2: + allclose, _, _ = _discontinuity_in_2d_bounds(self.bounds, + abs_tol=atol) + return allclose else: return False def contiguous_bounds(self): """ - Returns the N+1 bound values for a contiguous bounded coordinate + Returns the N+1 bound values for a contiguous bounded 1D coordinate of length N. + Returns the (N+1, M+1) bound values for a contiguous bounded 2D + coordinate of shape (N, M). + + Assumes input is contiguous. + .. note:: - If the coordinate is does not have bounds, this method will + If the coordinate does not have bounds, this method will return bounds positioned halfway between the coordinate's points. """ @@ -1003,8 +1099,11 @@ def contiguous_bounds(self): self._sanity_check_contiguous() bounds = self.bounds - c_bounds = np.resize(bounds[:, 0], bounds.shape[0] + 1) - c_bounds[-1] = bounds[-1, 1] + if self.ndim == 1: + c_bounds = np.resize(bounds[:, 0], bounds.shape[0] + 1) + c_bounds[-1] = bounds[-1, 1] + elif self.ndim == 2: + c_bounds = _get_2d_coord_bound_grid(bounds) return c_bounds def is_monotonic(self): diff --git a/lib/iris/plot.py b/lib/iris/plot.py index 64ce54c04c..5e00d5f713 100644 --- a/lib/iris/plot.py +++ b/lib/iris/plot.py @@ -31,16 +31,17 @@ import cartopy.crs as ccrs import cartopy.mpl.geoaxes from cartopy.geodesic import Geodesic +import cftime import matplotlib.axes import matplotlib.collections as mpl_collections import matplotlib.dates as mpl_dates import matplotlib.pyplot as plt from matplotlib.offsetbox import AnchoredText import matplotlib.ticker as mpl_ticker -import cftime import matplotlib.transforms as mpl_transforms import numpy as np import numpy.ma as ma +import warnings import iris.cube import iris.analysis.cartography as cartography @@ -51,7 +52,6 @@ import iris.palette from iris.util import _meshgrid - # Cynthia Brewer citation text. BREWER_CITE = 'Colours based on ColorBrewer.org' @@ -65,7 +65,7 @@ def names(coords): if isinstance(coord, int): result.append('dim={}'.format(coord)) else: - result.append(coord.name()) + result.append(cube.coord.name()) return ', '.join(result) def as_coord(coord): @@ -95,15 +95,17 @@ def get_span(coord): else: span = set(cube.coord_dims(coord)) return span + spans = list(map(get_span, coords)) for span, coord in zip(spans, coords): if not span: msg = 'The coordinate {!r} doesn\'t span a data dimension.' raise ValueError(msg.format(coord.name())) - if mode == iris.coords.BOUND_MODE and len(span) != 1: - raise ValueError('The coordinate {!r} is multi-dimensional and' - ' cannot be used in a cell-based plot.' - .format(coord.name())) + if mode == iris.coords.BOUND_MODE and len(span) not in [1, 2]: + raise ValueError('The coordinate {!r} has {} dimensions.' + 'Cell-based plotting is only supporting for' + 'coordinates with one or two dimensions.' + .format(coord.name()), len(span)) # Check the combination of coordinates spans enough (ndims) data # dimensions. @@ -125,11 +127,13 @@ def get_span(coord): # Note the use of `reversed` to convert from the X-then-Y # convention of the end-user API to the V-then-U convention used by # the plotting routines. + coords = names(coords) + plot_coords = list(reversed(coords)) return PlotDefn(plot_coords, transpose) -def _valid_bound_coord(coord): +def _valid_bound_dim_coord(coord): result = None if coord and coord.ndim == 1 and coord.nbounds: result = coord @@ -154,7 +158,7 @@ def _get_plot_defn(cube, mode, ndims=2): # When appropriate, restrict to 1D with bounds. if mode == iris.coords.BOUND_MODE: - coords = list(map(_valid_bound_coord, coords)) + coords = list(map(_valid_bound_dim_coord, coords)) def guess_axis(coord): axis = None @@ -172,6 +176,19 @@ def guess_axis(coord): aux_coords.sort(key=lambda coord: coord._as_defn()) coords[dim] = aux_coords[0] + # If plotting a 2 dimensional plot, check for 2d coordinates + if ndims == 2: + missing_dims = [dim for dim, coord in enumerate(coords) if + coord is None] + if missing_dims: + # Note that this only picks up coordinates that span the dims + two_dim_coords = cube.coords(dimensions=missing_dims) + two_dim_coords = [coord for coord in two_dim_coords + if coord.ndim == 2] + if len(two_dim_coords) >= 2: + two_dim_coords.sort(key=lambda coord: coord._as_defn()) + coords = two_dim_coords[:2] + if mode == iris.coords.POINT_MODE: # Allow multi-dimensional aux_coords to override the dim_coords # along the Z axis. This results in a preference for using the @@ -195,6 +212,7 @@ def sort_key(coord): return (order.get(axis, 0), coords.index(coord), coord and coord.name()) + sorted_coords = sorted(coords, key=sort_key) transpose = (sorted_coords != coords) @@ -256,6 +274,61 @@ def _invert_yaxis(v_coord, axes=None): axes.invert_yaxis() +def _check_contiguity_and_bounds(coord, data, abs_tol=1e-4, transpose=False): + """ + Check that the discontinuous bounds occur where the data is masked. + + If discontinuity occurs but data is masked, raise warning + If discontinuity occurs and data is NOT masked, raise error + + Args: + coords: + Array of the bounds of a 2D coord, of shape (X,Y,4) + data: + data of the the cube we are plotting + abs_tol: + tolerance when checking the contiguity + + """ + if transpose: + bounds = coord.bounds.T + else: + bounds = coord.bounds + + both_dirs_contiguous, diffs_along_x, diffs_along_y = \ + iris.coords._discontinuity_in_2d_bounds(bounds, abs_tol=abs_tol) + + if not both_dirs_contiguous: + + # True where data exists + mask_invert = np.logical_not(data.mask) + + # Where a discontinuity occurs the grid will not be created correctly. + # This does not matter if the data is masked as this is not plotted. + # So for places where data exists (opposite of the mask) AND a + # diff exists. If any exist, raise a warning + not_masked_at_discontinuity_along_x = np.any( + np.logical_and(mask_invert[:, :-1], diffs_along_x)) + + not_masked_at_discontinuity_along_y = np.any( + np.logical_and(mask_invert[:-1, ], diffs_along_y)) + + # If discontinuity occurs but not masked, any grid will be created + # incorrectly, so raise a warning + if not_masked_at_discontinuity_along_x or \ + not_masked_at_discontinuity_along_y: + raise ValueError('The bounds of the {} coordinate are not' + ' contiguous and data is not masked where the' + ' discontiguity occurs. Not able to create a' + ' suitable mesh to give to' + ' Matplotlib'.format(coord.name())) + else: + msg = 'The bounds of the {} coordinate are not contiguous. ' \ + 'However, data is masked where the discontiguity occurs ' \ + 'so plotting anyway.'.format(coord.name()) + warnings.warn(msg) + + def _draw_2d_from_bounds(draw_method_name, cube, *args, **kwargs): # NB. In the interests of clarity we use "u" and "v" to refer to the # horizontal and vertical axes on the matplotlib plot. @@ -263,10 +336,25 @@ def _draw_2d_from_bounds(draw_method_name, cube, *args, **kwargs): # Get & remove the coords entry from kwargs. coords = kwargs.pop('coords', None) if coords is not None: - plot_defn = _get_plot_defn_custom_coords_picked(cube, coords, mode) + plot_defn = _get_plot_defn_custom_coords_picked(cube, coords, mode, + ndims=2) else: plot_defn = _get_plot_defn(cube, mode, ndims=2) + two_dim_contig_atol = kwargs.pop('two_dim_coord_contiguity_atol', + 1e-4) + for coord in plot_defn.coords: + if coord.ndim == 2 and coord.has_bounds(): + try: + _check_contiguity_and_bounds(coord, data=cube.data, + abs_tol=two_dim_contig_atol) + except ValueError: + if _check_contiguity_and_bounds(coord, data=cube.data, + abs_tol=two_dim_contig_atol, + transpose=True)\ + is True: + plot_defn.transpose = True + if _can_draw_map(plot_defn.coords): result = _map_common(draw_method_name, None, iris.coords.BOUND_MODE, cube, plot_defn, *args, **kwargs) @@ -286,7 +374,7 @@ def _draw_2d_from_bounds(draw_method_name, cube, *args, **kwargs): string_axes = {} for coord, axis_name, data_dim in zip([u_coord, v_coord], - ['xaxis', 'yaxis'], + ['plotxaxis', 'yaxis'], [1, 0]): if coord is None: values = np.arange(data.shape[data_dim] + 1) @@ -309,7 +397,16 @@ def _draw_2d_from_bounds(draw_method_name, cube, *args, **kwargs): plot_arrays.append(values) u, v = plot_arrays - u, v = _broadcast_2d(u, v) + + # If the data is tranposed, 2D coordinates will also need to be + # tranposed. + if u.ndim == v.ndim == 2 and plot_defn.transpose is True: + u = u.T + v = v.T + + if u.ndim == v.ndim == 1: + u, v = _broadcast_2d(u, v) + axes = kwargs.pop('axes', None) draw_method = getattr(axes if axes else plt, draw_method_name) result = draw_method(u, v, data, *args, **kwargs) @@ -434,8 +531,8 @@ def _fixup_dates(coord, values): raise IrisError(msg) r = [nc_time_axis.CalendarDateTime( - cftime.datetime(*date), coord.units.calendar) - for date in dates] + cftime.datetime(*date), coord.units.calendar) + for date in dates] values = np.empty(len(r), dtype=object) values[:] = r return values @@ -675,8 +772,8 @@ def _ensure_cartopy_axes_and_determine_kwargs(x_coord, y_coord, kwargs): axes = kwargs.get('axes') if axes is None: if (isinstance(cs, iris.coord_systems.RotatedGeogCS) and - x_coord.points.max() > 180 and x_coord.points.max() < 360 and - x_coord.points.min() > 0): + x_coord.points.max() > 180 and x_coord.points.max() < 360 and + x_coord.points.min() > 0): # The RotatedGeogCS has 0 - 360 extent, different from the # assumptions made by Cartopy: rebase longitudes for the map axes # to set the datum longitude to the International Date Line. @@ -721,17 +818,22 @@ def _map_common(draw_method_name, arg_func, mode, cube, plot_defn, else: raise ValueError("Expected 1D or 2D XY coords") else: - try: - x, y = _meshgrid(x_coord.contiguous_bounds(), - y_coord.contiguous_bounds()) - # Exception translation. - except iris.exceptions.CoordinateMultiDimError: - raise ValueError("Could not get XY grid from bounds. " - "X or Y coordinate not 1D.") - except ValueError: - raise ValueError("Could not get XY grid from bounds. " - "X or Y coordinate doesn't have 2 bounds " - "per point.") + if not x_coord.ndim == y_coord.ndim == 2: + try: + x, y = _meshgrid(x_coord.contiguous_bounds(), + y_coord.contiguous_bounds()) + # Exception translation. + except iris.exceptions.CoordinateMultiDimError: + raise ValueError("Expected two 1D coords. Could not get XY" + " grid from bounds. X or Y coordinate not" + " 1D.") + except ValueError: + raise ValueError("Could not get XY grid from bounds. " + "X or Y coordinate doesn't have 2 bounds " + "per point.") + else: + x = x_coord.contiguous_bounds() + y = y_coord.contiguous_bounds() # Obtain the data array. data = cube.data @@ -1045,6 +1147,9 @@ def pcolormesh(cube, *args, **kwargs): * axes: the :class:`matplotlib.axes.Axes` to use for drawing. Defaults to the current axes if none provided. + * two_dim_coord_contiguity_atol: absolute tolerance when checking for + contiguity between cells in a two dimensional coordinate. + See :func:`matplotlib.pyplot.pcolormesh` for details of other valid keyword arguments. @@ -1073,6 +1178,7 @@ def points(cube, *args, **kwargs): keyword arguments. """ + def _scatter_args(u, v, data, *args, **kwargs): return ((u, v) + args, kwargs) diff --git a/lib/iris/tests/stock.py b/lib/iris/tests/stock.py index 25df8eda33..1b55510d2b 100644 --- a/lib/iris/tests/stock.py +++ b/lib/iris/tests/stock.py @@ -1,4 +1,4 @@ -# (C) British Crown Copyright 2010 - 2017, Met Office +# (C) British Crown Copyright 2010 - 2018, Met Office # # This file is part of Iris. # @@ -46,12 +46,12 @@ def lat_lon_cube(): """ cube = Cube(np.arange(12, dtype=np.int32).reshape((3, 4))) cs = GeogCS(6371229) - coord = iris.coords.DimCoord(points=np.array([-1, 0, 1], dtype=np.int32), + coord = DimCoord(points=np.array([-1, 0, 1], dtype=np.int32), standard_name='latitude', units='degrees', coord_system=cs) cube.add_dim_coord(coord, 0) - coord = iris.coords.DimCoord(points=np.array([-1, 0, 1, 2], dtype=np.int32), + coord = DimCoord(points=np.array([-1, 0, 1, 2], dtype=np.int32), standard_name='longitude', units='degrees', coord_system=cs) @@ -99,7 +99,7 @@ def simple_1d(with_bounds=True): cube.units = '1' points = np.arange(11, dtype=np.int32) + 1 bounds = np.column_stack([np.arange(11, dtype=np.int32), np.arange(11, dtype=np.int32) + 1]) - coord = iris.coords.DimCoord(points, long_name='foo', units='1', bounds=bounds) + coord = DimCoord(points, long_name='foo', units='1', bounds=bounds) cube.add_dim_coord(coord, 0) return cube @@ -124,12 +124,15 @@ def simple_2d(with_bounds=True): cube = Cube(np.arange(12, dtype=np.int32).reshape((3, 4))) cube.long_name = 'thingness' cube.units = '1' - y_points = np.array([ 2.5, 7.5, 12.5]) + y_points = np.array([2.5, 7.5, 12.5]) y_bounds = np.array([[0, 5], [5, 10], [10, 15]], dtype=np.int32) - y_coord = iris.coords.DimCoord(y_points, long_name='bar', units='1', bounds=y_bounds if with_bounds else None) + y_coord = DimCoord(y_points, long_name='bar', units='1', + bounds=y_bounds if with_bounds else None) x_points = np.array([ -7.5, 7.5, 22.5, 37.5]) - x_bounds = np.array([[-15, 0], [0, 15], [15, 30], [30, 45]], dtype=np.int32) - x_coord = iris.coords.DimCoord(x_points, long_name='foo', units='1', bounds=x_bounds if with_bounds else None) + x_bounds = np.array([[-15, 0], [0, 15], [15, 30], [30, 45]], + dtype=np.int32) + x_coord = DimCoord(x_points, long_name='foo', units='1', + bounds=x_bounds if with_bounds else None) cube.add_dim_coord(y_coord, 0) cube.add_dim_coord(x_coord, 1) @@ -191,9 +194,8 @@ def simple_3d_w_multidim_coords(with_bounds=True): [[5, 15], [15, 20], [20, 35], [35, 50]], [[10, 20], [20, 25], [25, 40], [40, 60]]], dtype=np.int32) - y_coord = iris.coords.AuxCoord(points=y_points, long_name='bar', - units='1', - bounds=y_bounds if with_bounds else None) + y_coord = AuxCoord(points=y_points, long_name='bar', units='1', + bounds=y_bounds if with_bounds else None) x_points = np.array([[-7.5, 7.5, 22.5, 37.5], [-12.5, 4., 26.5, 47.5], [2.5, 14., 36.5, 44.]]) @@ -201,12 +203,10 @@ def simple_3d_w_multidim_coords(with_bounds=True): [[-25, 0], [0, 8], [8, 45], [45, 50]], [[-5, 10], [10, 18], [18, 55], [18, 70]]], dtype=np.int32) - x_coord = iris.coords.AuxCoord(points=x_points, long_name='foo', - units='1', - bounds=x_bounds if with_bounds else None) - wibble_coord = iris.coords.DimCoord(np.array([10., 30.], - dtype=np.float32), - long_name='wibble', units='1') + x_coord = AuxCoord(points=x_points, long_name='foo', units='1', + bounds=x_bounds if with_bounds else None) + wibble_coord = DimCoord(np.array([10., 30.], dtype=np.float32), + long_name='wibble', units='1') cube.add_dim_coord(wibble_coord, [0]) cube.add_aux_coord(y_coord, [1, 2]) @@ -238,13 +238,13 @@ def simple_3d(): cube = Cube(np.arange(24, dtype=np.int32).reshape((2, 3, 4))) cube.long_name = 'thingness' cube.units = '1' - wibble_coord = iris.coords.DimCoord(np.array([10., 30.], + wibble_coord = DimCoord(np.array([10., 30.], dtype=np.float32), long_name='wibble', units='1') - lon = iris.coords.DimCoord([-180, -90, 0, 90], + lon = DimCoord([-180, -90, 0, 90], standard_name='longitude', units='degrees', circular=True) - lat = iris.coords.DimCoord([90, 0, -90], + lat = DimCoord([90, 0, -90], standard_name='latitude', units='degrees') cube.add_dim_coord(wibble_coord, [0]) cube.add_dim_coord(lat, [1]) @@ -294,39 +294,46 @@ def track_1d(duplicate_x=False): array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) """ - cube = Cube(np.arange(11, dtype=np.int32), standard_name='air_temperature', units='K') - bounds = np.column_stack([np.arange(11, dtype=np.int32), np.arange(11, dtype=np.int32) + 1]) + cube = Cube(np.arange(11, dtype=np.int32), standard_name='air_temperature', + units='K') + bounds = np.column_stack([np.arange(11, dtype=np.int32), + np.arange(11, dtype=np.int32) + 1]) pts = bounds[:, 1] - coord = iris.coords.AuxCoord(pts, 'projection_x_coordinate', units='1', bounds=bounds) + coord = AuxCoord(pts, 'projection_x_coordinate', units='1', bounds=bounds) cube.add_aux_coord(coord, [0]) if duplicate_x: - coord = iris.coords.AuxCoord(pts, 'projection_x_coordinate', units='1', bounds=bounds) + coord = AuxCoord(pts, 'projection_x_coordinate', units='1', + bounds=bounds) cube.add_aux_coord(coord, [0]) - coord = iris.coords.AuxCoord(pts * 2, 'projection_y_coordinate', units='1', bounds=bounds * 2) + coord = AuxCoord(pts * 2, 'projection_y_coordinate', units='1', + bounds=bounds * 2) cube.add_aux_coord(coord, 0) return cube def simple_2d_w_multidim_and_scalars(): data = np.arange(50, dtype=np.int32).reshape((5, 10)) - cube = iris.cube.Cube(data, long_name='test 2d dimensional cube', units='meters') + cube = iris.cube.Cube(data, long_name='test 2d dimensional cube', + units='meters') # DimCoords - dim1 = iris.coords.DimCoord(np.arange(5, dtype=np.float32) * 5.1 + 3.0, long_name='dim1', units='meters') - dim2 = iris.coords.DimCoord(np.arange(10, dtype=np.int32), long_name='dim2', units='meters', - bounds=np.arange(20, dtype=np.int32).reshape(10, 2)) + dim1 = DimCoord(np.arange(5, dtype=np.float32) * 5.1 + 3.0, + long_name='dim1', units='meters') + dim2 = DimCoord(np.arange(10, dtype=np.int32), + long_name='dim2', units='meters', + bounds=np.arange(20, dtype=np.int32).reshape(10, 2)) # Scalars - an_other = iris.coords.AuxCoord(3.0, long_name='an_other', units='meters') - yet_an_other = iris.coords.DimCoord(23.3, standard_name='air_temperature', - long_name='custom long name', - var_name='custom_var_name', - units='K') + an_other = AuxCoord(3.0, long_name='an_other', units='meters') + yet_an_other = DimCoord(23.3, standard_name='air_temperature', + long_name='custom long name', + var_name='custom_var_name', units='K') # Multidim - my_multi_dim_coord = iris.coords.AuxCoord(np.arange(50, dtype=np.int32).reshape(5, 10), - long_name='my_multi_dim_coord', units='1', - bounds=np.arange(200, dtype=np.int32).reshape(5, 10, 4)) + my_multi_dim_coord = AuxCoord(np.arange(50, dtype=np.int32).reshape(5, 10), + long_name='my_multi_dim_coord', units='1', + bounds=np.arange(200, dtype=np.int32). + reshape(5, 10, 4)) cube.add_dim_coord(dim1, 0) cube.add_dim_coord(dim2, 1) @@ -361,18 +368,21 @@ def hybrid_height(): """ data = np.arange(12, dtype='i8').reshape((3, 4)) - orography = icoords.AuxCoord([10, 25, 50, 5], standard_name='surface_altitude', units='m') - model_level = icoords.AuxCoord([2, 1, 0], standard_name='model_level_number') - level_height = icoords.DimCoord([100, 50, 10], long_name='level_height', - units='m', attributes={'positive': 'up'}, - bounds=[[150, 75], [75, 20], [20, 0]]) - sigma = icoords.AuxCoord([0.8, 0.9, 0.95], long_name='sigma', - bounds=[[0.7, 0.85], [0.85, 0.97], [0.97, 1.0]]) - hybrid_height = iris.aux_factory.HybridHeightFactory(level_height, sigma, orography) + orography = AuxCoord([10, 25, 50, 5], standard_name='surface_altitude', + units='m') + model_level = AuxCoord([2, 1, 0], standard_name='model_level_number') + level_height = DimCoord([100, 50, 10], long_name='level_height', + units='m', attributes={'positive': 'up'}, + bounds=[[150, 75], [75, 20], [20, 0]]) + sigma = AuxCoord([0.8, 0.9, 0.95], long_name='sigma', + bounds=[[0.7, 0.85], [0.85, 0.97], [0.97, 1.0]]) + hybrid_height = iris.aux_factory.HybridHeightFactory(level_height, + sigma, orography) cube = iris.cube.Cube(data, standard_name='air_temperature', units='K', dim_coords_and_dims=[(level_height, 0)], - aux_coords_and_dims=[(orography, 1), (model_level, 0), (sigma, 0)], + aux_coords_and_dims=[(orography, 1), + (model_level, 0), (sigma, 0)], aux_factories=[hybrid_height]) return cube @@ -381,25 +391,25 @@ def simple_4d_with_hybrid_height(): cube = iris.cube.Cube(np.arange(3*4*5*6, dtype='i8').reshape(3, 4, 5, 6), "air_temperature", units="K") - cube.add_dim_coord(iris.coords.DimCoord(np.arange(3, dtype='i8'), "time", - units="hours since epoch"), 0) - cube.add_dim_coord(iris.coords.DimCoord(np.arange(4, dtype='i8')+10, - "model_level_number", units="1"), 1) - cube.add_dim_coord(iris.coords.DimCoord(np.arange(5, dtype='i8')+20, - "grid_latitude", - units="degrees"), 2) - cube.add_dim_coord(iris.coords.DimCoord(np.arange(6, dtype='i8')+30, - "grid_longitude", - units="degrees"), 3) - - cube.add_aux_coord(iris.coords.AuxCoord(np.arange(4, dtype='i8')+40, - long_name="level_height", - units="m"), 1) - cube.add_aux_coord(iris.coords.AuxCoord(np.arange(4, dtype='i8')+50, - long_name="sigma", units="1"), 1) - cube.add_aux_coord(iris.coords.AuxCoord(np.arange(5*6, dtype='i8').reshape(5, 6)+100, - long_name="surface_altitude", - units="m"), [2, 3]) + cube.add_dim_coord(DimCoord(np.arange(3, dtype='i8'), "time", + units="hours since epoch"), 0) + cube.add_dim_coord(DimCoord(np.arange(4, dtype='i8')+10, + "model_level_number", units="1"), 1) + cube.add_dim_coord(DimCoord(np.arange(5, dtype='i8')+20, + "grid_latitude", + units="degrees"), 2) + cube.add_dim_coord(DimCoord(np.arange(6, dtype='i8')+30, + "grid_longitude", + units="degrees"), 3) + + cube.add_aux_coord(AuxCoord(np.arange(4, dtype='i8')+40, + long_name="level_height", + units="m"), 1) + cube.add_aux_coord(AuxCoord(np.arange(4, dtype='i8')+50, + long_name="sigma", units="1"), 1) + cube.add_aux_coord(AuxCoord(np.arange(5*6, dtype='i8').reshape(5, 6)+100, + long_name="surface_altitude", + units="m"), [2, 3]) cube.add_aux_factory(iris.aux_factory.HybridHeightFactory( delta=cube.coord("level_height"), @@ -449,7 +459,8 @@ def realistic_4d(): Returns a realistic 4d cube. >>> print(repr(realistic_4d())) - + """ # the stock arrays were created in Iris 0.8 with: @@ -477,7 +488,8 @@ def realistic_4d(): if not os.path.isfile(data_path): raise IOError('Test data is not available at {}.'.format(data_path)) r = np.load(data_path) - # sort the arrays based on the order they were originally given. The names given are of the form 'arr_1' or 'arr_10' + # sort the arrays based on the order they were originally given. + # The names given are of the form 'arr_1' or 'arr_10' _, arrays = zip(*sorted(six.iteritems(r), key=lambda item: int(item[0][4:]))) lat_pts, lat_bnds, lon_pts, lon_bnds, level_height_pts, \ @@ -486,25 +498,37 @@ def realistic_4d(): ll_cs = RotatedGeogCS(37.5, 177.5, ellipsoid=GeogCS(6371229.0)) - lat = icoords.DimCoord(lat_pts, standard_name='grid_latitude', units='degrees', - bounds=lat_bnds, coord_system=ll_cs) - lon = icoords.DimCoord(lon_pts, standard_name='grid_longitude', units='degrees', - bounds=lon_bnds, coord_system=ll_cs) + lat = icoords.DimCoord(lat_pts, standard_name='grid_latitude', + units='degrees', bounds=lat_bnds, + coord_system=ll_cs) + lon = icoords.DimCoord(lon_pts, standard_name='grid_longitude', + units='degrees', bounds=lon_bnds, + coord_system=ll_cs) level_height = icoords.DimCoord(level_height_pts, long_name='level_height', units='m', bounds=level_height_bnds, attributes={'positive': 'up'}) - model_level = icoords.DimCoord(model_level_pts, standard_name='model_level_number', + model_level = icoords.DimCoord(model_level_pts, + standard_name='model_level_number', units='1', attributes={'positive': 'up'}) - sigma = icoords.AuxCoord(sigma_pts, long_name='sigma', units='1', bounds=sigma_bnds) - orography = icoords.AuxCoord(orography, standard_name='surface_altitude', units='m') - time = icoords.DimCoord(time_pts, standard_name='time', units='hours since 1970-01-01 00:00:00') - forecast_period = icoords.DimCoord(forecast_period_pts, standard_name='forecast_period', units='hours') + sigma = icoords.AuxCoord(sigma_pts, long_name='sigma', units='1', + bounds=sigma_bnds) + orography = icoords.AuxCoord(orography, standard_name='surface_altitude', + units='m') + time = icoords.DimCoord(time_pts, standard_name='time', + units='hours since 1970-01-01 00:00:00') + forecast_period = icoords.DimCoord(forecast_period_pts, + standard_name='forecast_period', + units='hours') - hybrid_height = iris.aux_factory.HybridHeightFactory(level_height, sigma, orography) + hybrid_height = iris.aux_factory.HybridHeightFactory(level_height, + sigma, orography) - cube = iris.cube.Cube(data, standard_name='air_potential_temperature', units='K', - dim_coords_and_dims=[(time, 0), (model_level, 1), (lat, 2), (lon, 3)], - aux_coords_and_dims=[(orography, (2, 3)), (level_height, 1), (sigma, 1), + cube = iris.cube.Cube(data, standard_name='air_potential_temperature', + units='K', + dim_coords_and_dims=[(time, 0), (model_level, 1), + (lat, 2), (lon, 3)], + aux_coords_and_dims=[(orography, (2, 3)), + (level_height, 1), (sigma, 1), (forecast_period, None)], attributes={'source': 'Iris test case'}, aux_factories=[hybrid_height]) @@ -516,7 +540,8 @@ def realistic_4d_no_derived(): Returns a realistic 4d cube without hybrid height >>> print(repr(realistic_4d())) - + """ cube = realistic_4d() @@ -534,24 +559,30 @@ def realistic_4d_w_missing_data(): data_archive = np.load(data_path) data = ma.masked_array(data_archive['arr_0'], mask=data_archive['arr_1']) - # sort the arrays based on the order they were originally given. The names given are of the form 'arr_1' or 'arr_10' + # sort the arrays based on the order they were originally given. + # The names given are of the form 'arr_1' or 'arr_10' ll_cs = GeogCS(6371229) - lat = iris.coords.DimCoord(np.arange(20, dtype=np.float32), standard_name='grid_latitude', - units='degrees', coord_system=ll_cs) - lon = iris.coords.DimCoord(np.arange(20, dtype=np.float32), standard_name='grid_longitude', - units='degrees', coord_system=ll_cs) - time = iris.coords.DimCoord([1000., 1003., 1006.], standard_name='time', - units='hours since 1970-01-01 00:00:00') - forecast_period = iris.coords.DimCoord([0.0, 3.0, 6.0], standard_name='forecast_period', units='hours') - pressure = iris.coords.DimCoord(np.array([ 800., 900., 1000.], dtype=np.float32), - long_name='pressure', units='hPa') + lat = DimCoord(np.arange(20, dtype=np.float32), + standard_name='grid_latitude', + units='degrees', coord_system=ll_cs) + lon = DimCoord(np.arange(20, dtype=np.float32), + standard_name='grid_longitude', + units='degrees', coord_system=ll_cs) + time = DimCoord([1000., 1003., 1006.], standard_name='time', + units='hours since 1970-01-01 00:00:00') + forecast_period = DimCoord([0.0, 3.0, 6.0], + standard_name='forecast_period', + units='hours') + pressure = DimCoord(np.array([800., 900., 1000.], dtype=np.float32), + long_name='pressure', units='hPa') cube = iris.cube.Cube(data, long_name='missing data test data', units='K', - dim_coords_and_dims=[(time, 0), (pressure, 1), (lat, 2), (lon, 3)], + dim_coords_and_dims=[(time, 0), (pressure, 1), + (lat, 2), (lon, 3)], aux_coords_and_dims=[(forecast_period, 0)], - attributes={'source':'Iris test case'}) + attributes={'source': 'Iris test case'}) return cube diff --git a/lib/iris/tests/unit/plot/test_2d_coords.py b/lib/iris/tests/unit/plot/test_2d_coords.py new file mode 100644 index 0000000000..3b8ca91652 --- /dev/null +++ b/lib/iris/tests/unit/plot/test_2d_coords.py @@ -0,0 +1,158 @@ +# (C) British Crown Copyright 2018, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +"""Unit tests for handling and plotting of 2-dimensional coordinates""" + +from __future__ import (absolute_import, division, print_function) +from six.moves import (filter, input, map, range, zip) # noqa + +# Import iris.tests first so that some things can be initialised before +# importing anything else. +import iris.tests as tests +import iris.coords as coords + +import numpy.ma as ma + +from iris.tests.stock import simple_2d_w_multidim_coords as cube_2dcoords +from iris.tests.stock import simple_3d_w_multidim_coords as cube3d_2dcoords +from iris.tests.stock import lat_lon_cube + +from orca_utils.plot_testing import testdata_2d_coords + +if tests.MPL_AVAILABLE: + import iris.plot as iplt + + +class Test_2d_coords_plot_defn_bound_mode(tests.IrisTest): + def setUp(self): + self.multidim_cube = cube_2dcoords() + self.overspan_cube = cube3d_2dcoords() + + # latlon_2d is a cube with 2d coords, 4 bounds per point, + # discontiguities in the bounds but masked data at the discontiguities. + self.latlon_2d = testdata_2d_coords.full2d_global() + testdata_2d_coords.make_bounds_discontiguous_at_point\ + (self.latlon_2d, 2, 2) + + # Take a latlon cube with 1D coords, broadcast the coords into 2D + # ones, then add ONE of them back into the cube in place of original: + single_dims = lat_lon_cube() + lon = single_dims.coord('longitude') + lat = single_dims.coord('latitude') + big_lon, big_lat = testdata_2d_coords.grid_coords_2d_from_1d(lon, lat) + mixed_dims = single_dims.copy() + mixed_dims.remove_coord(lon) + # TODO Fix this coord addition: + # When adding an aux_coord, the function '_check_multidim_metadata' + # throws an error as it requires coord.shape to be (1, ) instead of + # (3, 4) or whatever. + mixed_dims.add_aux_coord(big_lon) + + # mixed_dims is now a cube with 2 1D dim coords and an additional + # 2D aux coord. + self.mixed_dims = mixed_dims + + self.mode = coords.BOUND_MODE + + def test_2d_coords_identified(self): + # Test that 2d coords are identified in the plot definition without + # having to be specified in the args without any errors. + cube = self.multidim_cube + defn = iplt._get_plot_defn(cube, mode=self.mode) + self.assertEqual([coord.name() for coord in defn.coords], + ['bar', 'foo']) + + def test_2d_coords_custom_picked(self): + # Test that 2d coords which are specified in the args will be + # accepted without any errors. + cube = self.multidim_cube + defn = iplt._get_plot_defn_custom_coords_picked(cube, ('foo', 'bar'), + self.mode) + self.assertEqual([coord.name() for coord in defn.coords], + ['bar', 'foo']) + + def test_2d_coords_as_integers(self): + # Test that if you pass in 2d coords as args in the form of integers, + # they will still be correctly identified without any errors. + cube = self.multidim_cube + defn = iplt._get_plot_defn_custom_coords_picked(cube, (0, 1), + self.mode) + self.assertEqual([coord for coord in defn.coords], + [1, 0]) + + def test_total_span_check(self): + # Test that an error is raised if a user tries to plot a 2d coord + # against a different coord, making total number of dimensions 3. + cube = self.overspan_cube + with self.assertRaises(ValueError): + iplt._get_plot_defn_custom_coords_picked(cube, ('wibble', 'foo'), + self.mode) + + def test_2dcoord_with_1dcoord(self): + # TODO Generate a cube with one 2d coord and one 1d coord + # TODO Try and plot them against each other + # TODO Find out where I can put a catch for this (if necessary) + cube = self.mixed_dims + with self.assertRaises(ValueError): + iplt._get_plot_defn_custom_coords_picked(cube, + ('latitude', 'longitude'), + self.mode) + + def test_map_common_not_enough_bounds(self): + # Test that a lat-lon cube with 2d coords and 2 bounds per point + # throws an error in contiguity checks. + cube = self.multidim_cube + cube.coord('foo').rename('longitude') + cube.coord('bar').rename('latitude') + with self.assertRaises(ValueError): + plot_defn = iplt._get_plot_defn(cube, self.mode) + iplt._map_common('pcolor', None, self.mode, cube, plot_defn) + + def test_map_common_2d(self): + # Test that a cube with 2d coords can be plotted as a map. + cube = self.latlon_2d + # Get necessary variables from _get_plot_defn to check that the test + # case will be accepted by _map_common. + plot_defn = iplt._get_plot_defn(cube, self.mode) + result = iplt._map_common('pcolor', None, self.mode, cube, plot_defn) + self.assertTrue(result) + + def test_discontiguous_masked(self): + # Test that a contiguity check will raise a warning (not an error) for + # discontiguous bounds but appropriately masked data. + cube = self.latlon_2d + coord = cube.coord('longitude') + msg = 'The bounds of the longitude coordinate are not contiguous. ' \ + 'However, data is masked where the discontiguity occurs so ' \ + 'plotting anyway.' + with self.assertWarnsRegexp(msg): + iplt._check_contiguity_and_bounds(coord, cube.data) + + def test_discontiguous_unmasked(self): + # Check that an error occurs when the contiguity check finds + # discontiguous bounds but unmasked data. + cube = self.latlon_2d + cube.data.mask = ma.nomask + coord = cube.coord('longitude') + with self.assertRaises(ValueError): + iplt._check_contiguity_and_bounds(coord, cube.data) + + def test_draw_2d_from_bounds(self): + # Test this function will not raise an error even with our most + # awkward but supported cube. + cube = self.latlon_2d + result = iplt._draw_2d_from_bounds('pcolormesh', cube) + self.assertTrue(result) From 723d45b1fe40afabcae54ff1d2868805e77da58a Mon Sep 17 00:00:00 2001 From: Corinne Bosley Date: Tue, 31 Jul 2018 09:44:38 +0100 Subject: [PATCH 14/40] some words for docs, WIP probably --- docs/iris/src/whatsnew/2.2.rst | 26 ++++++++++++++++++++++++++ lib/iris/plot.py | 15 +++++++++++++-- 2 files changed, 39 insertions(+), 2 deletions(-) create mode 100644 docs/iris/src/whatsnew/2.2.rst diff --git a/docs/iris/src/whatsnew/2.2.rst b/docs/iris/src/whatsnew/2.2.rst new file mode 100644 index 0000000000..4165ebf43d --- /dev/null +++ b/docs/iris/src/whatsnew/2.2.rst @@ -0,0 +1,26 @@ +What's New in Iris 2.2.0 +************************ + +:Release: 2.2.0a0 +:Date: + + +This document explains the new/changed features of Iris in the alpha release + of version 2.2.0 +(:doc:`View all changes `). + + +Iris 2.2.0 Features +=================== +.. _showcase: + +.. admonition:: 2-Dimensional Coordinate Plotting + + The iris plot functions :func:`~iris.plot.pcolor` and + :func:`~iris.plot.pcolormesh` now accommodate the plotting of 2-dimensional + coordinates as well as 1-dimensional coordinates. + + To enable this feature, each coordinate passed in for plotting will be + automatically checked for contiguity. Coordinate bounds must either be + contiguous, or the cube's data must be masked at the discontiguities in + order to avoid plotting errors. diff --git a/lib/iris/plot.py b/lib/iris/plot.py index 5e00d5f713..a464f5393b 100644 --- a/lib/iris/plot.py +++ b/lib/iris/plot.py @@ -1107,7 +1107,10 @@ def outline(cube, coords=None, color='k', linewidth=None, axes=None): def pcolor(cube, *args, **kwargs): """ - Draws a pseudocolor plot based on the given Cube. + Draws a pseudocolor plot based on the given 2-dimensional Cube. + + The cube must have either two 1-dimensional coordinates or two + 2-dimensional coordinates with contiguous bounds to plot against each other. Kwargs: @@ -1121,6 +1124,11 @@ def pcolor(cube, *args, **kwargs): * axes: the :class:`matplotlib.axes.Axes` to use for drawing. Defaults to the current axes if none provided. + * two_dim_coord_contiguity_atol: absolute tolerance when checking for + contiguity between cells in a two dimensional coordinate. + + + See :func:`matplotlib.pyplot.pcolor` for details of other valid keyword arguments. @@ -1133,7 +1141,10 @@ def pcolor(cube, *args, **kwargs): def pcolormesh(cube, *args, **kwargs): """ - Draws a pseudocolor plot based on the given Cube. + Draws a pseudocolor plot based on the given 2-dimensional Cube. + + The cube must have either two 1-dimensional coordinates or two + 2-dimensional coordinates with contiguous bounds to plot against each other. Kwargs: From 0b8baf55437f2dd7882cddc1c96e93e415ac1a23 Mon Sep 17 00:00:00 2001 From: lbdreyer Date: Mon, 30 Jul 2018 18:54:27 +0100 Subject: [PATCH 15/40] Handle custom coords correctly --- lib/iris/coords.py | 15 ++++++++------- lib/iris/plot.py | 33 +++++++++++++-------------------- 2 files changed, 21 insertions(+), 27 deletions(-) diff --git a/lib/iris/coords.py b/lib/iris/coords.py index 57ae723205..8f1e48fd18 100644 --- a/lib/iris/coords.py +++ b/lib/iris/coords.py @@ -143,7 +143,7 @@ def __new__(cls, name_or_coord, minimum, maximum, 'groupby_point, groupby_slice') -def _discontinuity_in_2d_bounds(bds, abs_tol=1e-4): +def _discontiguity_in_2d_bounds(bds, abs_tol=1e-4): """ Check bounds of a 2-dimensional coordinate are contiguous Args: @@ -156,7 +156,7 @@ def _discontinuity_in_2d_bounds(bds, abs_tol=1e-4): absolute difference along the y axis """ - # Check form is (ny, nx, 4) + # Check bds has the shape (ny, nx, 4) if not bds.ndim == 3 and bds.shape[2] == 4: raise ValueError('2D coordinates must have 4 bounds per point ' 'for 2D coordinate plotting') @@ -1066,14 +1066,15 @@ def is_contiguous(self, rtol=1e-05, atol=1e-08): if self.has_bounds(): self._sanity_check_contiguous() if self.ndim == 1: - return np.allclose(self.bounds[1:, 0], self.bounds[:-1, 1], - rtol=rtol, atol=atol) + contiguous = np.allclose(self.bounds[1:, 0], + self.bounds[:-1, 1], + rtol=rtol, atol=atol) elif self.ndim == 2: - allclose, _, _ = _discontinuity_in_2d_bounds(self.bounds, + contiguous, _, _ = _discontiguity_in_2d_bounds(self.bounds, abs_tol=atol) - return allclose else: - return False + contiguous = False + return contiguous def contiguous_bounds(self): """ diff --git a/lib/iris/plot.py b/lib/iris/plot.py index a464f5393b..6ede9e1735 100644 --- a/lib/iris/plot.py +++ b/lib/iris/plot.py @@ -41,7 +41,6 @@ import matplotlib.transforms as mpl_transforms import numpy as np import numpy.ma as ma -import warnings import iris.cube import iris.analysis.cartography as cartography @@ -65,7 +64,7 @@ def names(coords): if isinstance(coord, int): result.append('dim={}'.format(coord)) else: - result.append(cube.coord.name()) + result.append(coord.name()) return ', '.join(result) def as_coord(coord): @@ -127,8 +126,6 @@ def get_span(coord): # Note the use of `reversed` to convert from the X-then-Y # convention of the end-user API to the V-then-U convention used by # the plotting routines. - coords = names(coords) - plot_coords = list(reversed(coords)) return PlotDefn(plot_coords, transpose) @@ -296,7 +293,7 @@ def _check_contiguity_and_bounds(coord, data, abs_tol=1e-4, transpose=False): bounds = coord.bounds both_dirs_contiguous, diffs_along_x, diffs_along_y = \ - iris.coords._discontinuity_in_2d_bounds(bounds, abs_tol=abs_tol) + iris.coords._discontiguity_in_2d_bounds(bounds, abs_tol=abs_tol) if not both_dirs_contiguous: @@ -322,11 +319,6 @@ def _check_contiguity_and_bounds(coord, data, abs_tol=1e-4, transpose=False): ' discontiguity occurs. Not able to create a' ' suitable mesh to give to' ' Matplotlib'.format(coord.name())) - else: - msg = 'The bounds of the {} coordinate are not contiguous. ' \ - 'However, data is masked where the discontiguity occurs ' \ - 'so plotting anyway.'.format(coord.name()) - warnings.warn(msg) def _draw_2d_from_bounds(draw_method_name, cube, *args, **kwargs): @@ -344,16 +336,17 @@ def _draw_2d_from_bounds(draw_method_name, cube, *args, **kwargs): two_dim_contig_atol = kwargs.pop('two_dim_coord_contiguity_atol', 1e-4) for coord in plot_defn.coords: - if coord.ndim == 2 and coord.has_bounds(): - try: - _check_contiguity_and_bounds(coord, data=cube.data, - abs_tol=two_dim_contig_atol) - except ValueError: - if _check_contiguity_and_bounds(coord, data=cube.data, - abs_tol=two_dim_contig_atol, - transpose=True)\ + if hasattr(coord, 'has_bounds'): + if coord.ndim == 2 and coord.has_bounds(): + try: + _check_contiguity_and_bounds(coord, data=cube.data, + abs_tol=two_dim_contig_atol) + except ValueError: + if _check_contiguity_and_bounds( + coord, data=cube.data, + abs_tol=two_dim_contig_atol, transpose=True) \ is True: - plot_defn.transpose = True + plot_defn.transpose = True if _can_draw_map(plot_defn.coords): result = _map_common(draw_method_name, None, iris.coords.BOUND_MODE, @@ -374,7 +367,7 @@ def _draw_2d_from_bounds(draw_method_name, cube, *args, **kwargs): string_axes = {} for coord, axis_name, data_dim in zip([u_coord, v_coord], - ['plotxaxis', 'yaxis'], + ['xaxis', 'yaxis'], [1, 0]): if coord is None: values = np.arange(data.shape[data_dim] + 1) From 759e357134d980c314a5da45e53d861cfbb6f2bc Mon Sep 17 00:00:00 2001 From: Corinne Bosley Date: Tue, 31 Jul 2018 11:15:12 +0100 Subject: [PATCH 16/40] suggested tweaks and corrections --- lib/iris/plot.py | 21 ++++++++++----------- lib/iris/tests/unit/plot/test_2d_coords.py | 9 ++++----- 2 files changed, 14 insertions(+), 16 deletions(-) diff --git a/lib/iris/plot.py b/lib/iris/plot.py index 6ede9e1735..d81354fe4c 100644 --- a/lib/iris/plot.py +++ b/lib/iris/plot.py @@ -288,9 +288,7 @@ def _check_contiguity_and_bounds(coord, data, abs_tol=1e-4, transpose=False): """ if transpose: - bounds = coord.bounds.T - else: - bounds = coord.bounds + data = data.T both_dirs_contiguous, diffs_along_x, diffs_along_y = \ iris.coords._discontiguity_in_2d_bounds(bounds, abs_tol=abs_tol) @@ -333,19 +331,18 @@ def _draw_2d_from_bounds(draw_method_name, cube, *args, **kwargs): else: plot_defn = _get_plot_defn(cube, mode, ndims=2) - two_dim_contig_atol = kwargs.pop('two_dim_coord_contiguity_atol', + twodim_contig_atol = kwargs.pop('two_dim_coord_contiguity_atol', 1e-4) for coord in plot_defn.coords: if hasattr(coord, 'has_bounds'): if coord.ndim == 2 and coord.has_bounds(): try: _check_contiguity_and_bounds(coord, data=cube.data, - abs_tol=two_dim_contig_atol) + abs_tol=twodim_contig_atol) except ValueError: - if _check_contiguity_and_bounds( - coord, data=cube.data, - abs_tol=two_dim_contig_atol, transpose=True) \ - is True: + if _check_contiguity_and_bounds(coord, data=cube.data, + abs_tol=twodim_contig_atol, + transpose=True) is True: plot_defn.transpose = True if _can_draw_map(plot_defn.coords): @@ -1103,7 +1100,8 @@ def pcolor(cube, *args, **kwargs): Draws a pseudocolor plot based on the given 2-dimensional Cube. The cube must have either two 1-dimensional coordinates or two - 2-dimensional coordinates with contiguous bounds to plot against each other. + 2-dimensional coordinates with contiguous bounds to plot against each + other. Kwargs: @@ -1137,7 +1135,8 @@ def pcolormesh(cube, *args, **kwargs): Draws a pseudocolor plot based on the given 2-dimensional Cube. The cube must have either two 1-dimensional coordinates or two - 2-dimensional coordinates with contiguous bounds to plot against each other. + 2-dimensional coordinates with contiguous bounds to plot against each + other. Kwargs: diff --git a/lib/iris/tests/unit/plot/test_2d_coords.py b/lib/iris/tests/unit/plot/test_2d_coords.py index 3b8ca91652..ce96157bf4 100644 --- a/lib/iris/tests/unit/plot/test_2d_coords.py +++ b/lib/iris/tests/unit/plot/test_2d_coords.py @@ -30,7 +30,7 @@ from iris.tests.stock import simple_3d_w_multidim_coords as cube3d_2dcoords from iris.tests.stock import lat_lon_cube -from orca_utils.plot_testing import testdata_2d_coords +from orca_utils.plot_testing import testdata_2d_coords as testdata if tests.MPL_AVAILABLE: import iris.plot as iplt @@ -43,16 +43,15 @@ def setUp(self): # latlon_2d is a cube with 2d coords, 4 bounds per point, # discontiguities in the bounds but masked data at the discontiguities. - self.latlon_2d = testdata_2d_coords.full2d_global() - testdata_2d_coords.make_bounds_discontiguous_at_point\ - (self.latlon_2d, 2, 2) + self.latlon_2d = testdata.full2d_global() + testdata.make_bounds_discontiguous_at_point(self.latlon_2d, 2, 2) # Take a latlon cube with 1D coords, broadcast the coords into 2D # ones, then add ONE of them back into the cube in place of original: single_dims = lat_lon_cube() lon = single_dims.coord('longitude') lat = single_dims.coord('latitude') - big_lon, big_lat = testdata_2d_coords.grid_coords_2d_from_1d(lon, lat) + big_lon, big_lat = testdata.grid_coords_2d_from_1d(lon, lat) mixed_dims = single_dims.copy() mixed_dims.remove_coord(lon) # TODO Fix this coord addition: From ade18ad70cbd9fba6fd21802136e5456b143f49d Mon Sep 17 00:00:00 2001 From: lbdreyer Date: Tue, 31 Jul 2018 11:24:30 +0100 Subject: [PATCH 17/40] Update docs to included 2.2 versions (#3110) --- docs/iris/src/_templates/index.html | 2 +- docs/iris/src/whatsnew/index.rst | 1 + lib/iris/__init__.py | 2 +- 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/iris/src/_templates/index.html b/docs/iris/src/_templates/index.html index 0c4b5d958f..31acded447 100644 --- a/docs/iris/src/_templates/index.html +++ b/docs/iris/src/_templates/index.html @@ -134,7 +134,7 @@ extra information on specific technical issues

  • -

    What's new in Iris 2.1?
    +

    What's new in Iris 2.2?
    recent changes in Iris's capabilities

    • diff --git a/docs/iris/src/whatsnew/index.rst b/docs/iris/src/whatsnew/index.rst index 104c5074ca..c3a34303f0 100644 --- a/docs/iris/src/whatsnew/index.rst +++ b/docs/iris/src/whatsnew/index.rst @@ -9,6 +9,7 @@ Iris versions. .. toctree:: :maxdepth: 2 + 2.2.rst 2.1.rst 2.0.rst 1.13.rst diff --git a/lib/iris/__init__.py b/lib/iris/__init__.py index 2e7eba5892..882b52aaea 100644 --- a/lib/iris/__init__.py +++ b/lib/iris/__init__.py @@ -121,7 +121,7 @@ def callback(cube, field, filename): # Iris revision. -__version__ = '2.2.0dev0' +__version__ = '2.2.0a0' # Restrict the names imported when using "from iris import *" __all__ = ['load', 'load_cube', 'load_cubes', 'load_raw', From 91d69fa660d7fd1b340e5a2dbdad7ce69082dc97 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 31 Jul 2018 11:36:09 +0100 Subject: [PATCH 18/40] Provide a test skipper for 2d coords WIP. (#3099) --- lib/iris/tests/__init__.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/lib/iris/tests/__init__.py b/lib/iris/tests/__init__.py index 8c64dc15b1..710a3fbac6 100644 --- a/lib/iris/tests/__init__.py +++ b/lib/iris/tests/__init__.py @@ -1196,6 +1196,12 @@ class MyPlotTests(test.GraphicsTest): 'Test(s) require "python-stratify", which is not available.') +SKIP_2D_TESTS = True +skip_2d = unittest.skipIf( + SKIP_2D_TESTS, + 'Test(s) broken by WIP on 2d coords support -- temporarily disabled.') + + def no_warnings(func): """ Provides a decorator to ensure that there are no warnings raised From b6b305f5a17a5a8d658e71835b0a91aba2bdf36f Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Mon, 30 Jul 2018 14:43:17 +0100 Subject: [PATCH 19/40] First working quiver+streamplot. --- lib/iris/plot.py | 133 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 133 insertions(+) diff --git a/lib/iris/plot.py b/lib/iris/plot.py index d81354fe4c..ac8da822cc 100644 --- a/lib/iris/plot.py +++ b/lib/iris/plot.py @@ -1189,6 +1189,139 @@ def _scatter_args(u, v, data, *args, **kwargs): *args, **kwargs) +def _vector_component_args(x_points, y_points, u_data, *args, **kwargs): + """ + Callback from _draw_2d_from_points for 'quiver' and 'streamlines'. + + Returns arguments (x, y, u, v), to be passed to the underlying matplotlib + call. + + "u_data" will always be "u_cube.data". + The matching "v_cube.data" component is stored in kwargs['_v_data']. + + """ + v_data = kwargs.pop('_v_data') + + # Rescale u+v values for plot distortion. + crs = kwargs.get('transform', None) + if crs: + if not isinstance(crs, (ccrs.PlateCarree, ccrs.RotatedPole)): + msg = ('Can only plot vectors provided in a lat-lon ' + 'projection, i.e. "cartopy.crs.PlateCarree" or ' + '"cartopy.crs.RotatedPole". This ' + "cubes coordinate system is {}.") + raise ValueError(msg.format(crs)) + # Given the above check, the Y points must be latitudes. + # We therefore **assume** they are in degrees : I'm not sure this + # is wise, but all the rest of this plot code does that, e.g. in + # _map_common. + # TODO: investigate degree units assumptions, here + elsewhere. + + # Implement a latitude scaling, but preserve the given magnitudes. + v_data = v_data.copy() + mags = np.sqrt(u_data * u_data + v_data * v_data) + v_data *= np.cos(np.deg2rad(y_points)) + scales = mags / np.sqrt(u_data * u_data + v_data * v_data) + u_data *= scales + v_data *= scales + + return ((x_points, y_points, u_data, v_data), kwargs) + + +def quiver(u_cube, v_cube, *args, **kwargs): + """ + Draws an arrow plot from two vector component cubes. + + Args: + + * u_cube, v_cube : (:class:`~iris.cube.Cube`) + u and v vector components. Must have same shape and units. + If the cubes have geographic coordinates, the values are treated as + true distance differentials, e.g. windspeeds, and *not* map coordinate + vectors. The components are aligned with the North and East of the + cube coordinate system. + + .. Note: + + At present, if u_cube and v_cube have geographic coordinates, then they + must be in a lat-lon coordinate system, though it may be a rotated one. + To transform wind values between coordinate systems, use + :func:`iris.analysis.cartography.rotate_vectors`. + To transform coordinate grid points, you will need to create + 2-dimensional arrays of x and y values. These can be transformed with + :meth:`cartopy.crs.CRS.transform_points`. + + Kwargs: + + * coords: (list of :class:`~iris.coords.Coord` or string) + Coordinates or coordinate names. Use the given coordinates as the axes + for the plot. The order of the given coordinates indicates which axis + to use for each, where the first element is the horizontal + axis of the plot and the second element is the vertical axis + of the plot. + + * axes: the :class:`matplotlib.axes.Axes` to use for drawing. + Defaults to the current axes if none provided. + + See :func:`matplotlib.pyplot.quiver` for details of other valid + keyword arguments. + + """ + # + # TODO: check u + v cubes for compatibility. + # + kwargs['_v_data'] = v_cube.data + return _draw_2d_from_points('quiver', _vector_component_args, u_cube, + *args, **kwargs) + + +def streamplot(u_cube, v_cube, *args, **kwargs): + """ + Draws a streamline plot from two vector component cubes. + + Args: + + * u_cube, v_cube : (:class:`~iris.cube.Cube`) + u and v vector components. Must have same shape and units. + If the cubes have geographic coordinates, the values are treated as + true distance differentials, e.g. windspeeds, and *not* map coordinate + vectors. The components are aligned with the North and East of the + cube coordinate system. + + .. Note: + + At present, if u_cube and v_cube have geographic coordinates, then they + must be in a lat-lon coordinate system, though it may be a rotated one. + To transform wind values between coordinate systems, use + :func:`iris.analysis.cartography.rotate_vectors`. + To transform coordinate grid points, you will need to create + 2-dimensional arrays of x and y values. These can be transformed with + :meth:`cartopy.crs.CRS.transform_points`. + + Kwargs: + + * coords: (list of :class:`~iris.coords.Coord` or string) + Coordinates or coordinate names. Use the given coordinates as the axes + for the plot. The order of the given coordinates indicates which axis + to use for each, where the first element is the horizontal + axis of the plot and the second element is the vertical axis + of the plot. + + * axes: the :class:`matplotlib.axes.Axes` to use for drawing. + Defaults to the current axes if none provided. + + See :func:`matplotlib.pyplot.quiver` for details of other valid + keyword arguments. + + """ + # + # TODO: check u + v cubes for compatibility. + # + kwargs['_v_data'] = v_cube.data + return _draw_2d_from_points('streamplot', _vector_component_args, u_cube, + *args, **kwargs) + + def plot(*args, **kwargs): """ Draws a line plot based on the given cube(s) or coordinate(s). From aee7f561d9e6dcc0cf5a8829d2f6db4bc8c17b25 Mon Sep 17 00:00:00 2001 From: lbdreyer Date: Tue, 31 Jul 2018 14:05:40 +0100 Subject: [PATCH 20/40] Change version number (#3118) --- docs/iris/src/_templates/layout.html | 2 +- docs/iris/src/whatsnew/2.2.rst | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/iris/src/_templates/layout.html b/docs/iris/src/_templates/layout.html index 53ae5105cb..8ecc35bade 100644 --- a/docs/iris/src/_templates/layout.html +++ b/docs/iris/src/_templates/layout.html @@ -37,7 +37,7 @@

    - Iris v2.1 + Iris v2.2

    A powerful, format-agnostic, community-driven Python library for analysing and diff --git a/docs/iris/src/whatsnew/2.2.rst b/docs/iris/src/whatsnew/2.2.rst index 4165ebf43d..97eba3dee6 100644 --- a/docs/iris/src/whatsnew/2.2.rst +++ b/docs/iris/src/whatsnew/2.2.rst @@ -6,7 +6,7 @@ What's New in Iris 2.2.0 This document explains the new/changed features of Iris in the alpha release - of version 2.2.0 +of version 2.2.0 (:doc:`View all changes `). From 297becf2ddb6b74128f700cffd65876cdeac002b Mon Sep 17 00:00:00 2001 From: Corinne Bosley Date: Tue, 31 Jul 2018 14:22:55 +0100 Subject: [PATCH 21/40] re-added bounds definition for contiguity check, removed unnecessary test --- lib/iris/plot.py | 2 + lib/iris/tests/unit/plot/test_2d_coords.py | 72 +++++++++++----------- 2 files changed, 38 insertions(+), 36 deletions(-) diff --git a/lib/iris/plot.py b/lib/iris/plot.py index ac8da822cc..8f11c8d116 100644 --- a/lib/iris/plot.py +++ b/lib/iris/plot.py @@ -290,6 +290,8 @@ def _check_contiguity_and_bounds(coord, data, abs_tol=1e-4, transpose=False): if transpose: data = data.T + bounds = coord.bounds + both_dirs_contiguous, diffs_along_x, diffs_along_y = \ iris.coords._discontiguity_in_2d_bounds(bounds, abs_tol=abs_tol) diff --git a/lib/iris/tests/unit/plot/test_2d_coords.py b/lib/iris/tests/unit/plot/test_2d_coords.py index ce96157bf4..bb27341baa 100644 --- a/lib/iris/tests/unit/plot/test_2d_coords.py +++ b/lib/iris/tests/unit/plot/test_2d_coords.py @@ -46,23 +46,23 @@ def setUp(self): self.latlon_2d = testdata.full2d_global() testdata.make_bounds_discontiguous_at_point(self.latlon_2d, 2, 2) - # Take a latlon cube with 1D coords, broadcast the coords into 2D - # ones, then add ONE of them back into the cube in place of original: - single_dims = lat_lon_cube() - lon = single_dims.coord('longitude') - lat = single_dims.coord('latitude') - big_lon, big_lat = testdata.grid_coords_2d_from_1d(lon, lat) - mixed_dims = single_dims.copy() - mixed_dims.remove_coord(lon) - # TODO Fix this coord addition: - # When adding an aux_coord, the function '_check_multidim_metadata' - # throws an error as it requires coord.shape to be (1, ) instead of - # (3, 4) or whatever. - mixed_dims.add_aux_coord(big_lon) - - # mixed_dims is now a cube with 2 1D dim coords and an additional - # 2D aux coord. - self.mixed_dims = mixed_dims + # # Take a latlon cube with 1D coords, broadcast the coords into 2D + # # ones, then add ONE of them back into the cube in place of original: + # single_dims = lat_lon_cube() + # lon = single_dims.coord('longitude') + # lat = single_dims.coord('latitude') + # big_lon, big_lat = testdata.grid_coords_2d_from_1d(lon, lat) + # mixed_dims = single_dims.copy() + # mixed_dims.remove_coord(lon) + # # TODO Fix this coord addition: + # # When adding an aux_coord, the function '_check_multidim_metadata' + # # throws an error as it requires coord.shape to be (1, ) instead of + # # (3, 4) or whatever. + # mixed_dims.add_aux_coord(big_lon) + # + # # mixed_dims is now a cube with 2 1D dim coords and an additional + # # 2D aux coord. + # self.mixed_dims = mixed_dims self.mode = coords.BOUND_MODE @@ -100,15 +100,15 @@ def test_total_span_check(self): iplt._get_plot_defn_custom_coords_picked(cube, ('wibble', 'foo'), self.mode) - def test_2dcoord_with_1dcoord(self): - # TODO Generate a cube with one 2d coord and one 1d coord - # TODO Try and plot them against each other - # TODO Find out where I can put a catch for this (if necessary) - cube = self.mixed_dims - with self.assertRaises(ValueError): - iplt._get_plot_defn_custom_coords_picked(cube, - ('latitude', 'longitude'), - self.mode) + # def test_2dcoord_with_1dcoord(self): + # # TODO Generate a cube with one 2d coord and one 1d coord + # # TODO Try and plot them against each other + # # TODO Find out where I can put a catch for this (if necessary) + # cube = self.mixed_dims + # with self.assertRaises(ValueError): + # iplt._get_plot_defn_custom_coords_picked(cube, + # ('latitude', 'longitude'), + # self.mode) def test_map_common_not_enough_bounds(self): # Test that a lat-lon cube with 2d coords and 2 bounds per point @@ -129,16 +129,16 @@ def test_map_common_2d(self): result = iplt._map_common('pcolor', None, self.mode, cube, plot_defn) self.assertTrue(result) - def test_discontiguous_masked(self): - # Test that a contiguity check will raise a warning (not an error) for - # discontiguous bounds but appropriately masked data. - cube = self.latlon_2d - coord = cube.coord('longitude') - msg = 'The bounds of the longitude coordinate are not contiguous. ' \ - 'However, data is masked where the discontiguity occurs so ' \ - 'plotting anyway.' - with self.assertWarnsRegexp(msg): - iplt._check_contiguity_and_bounds(coord, cube.data) + # def test_discontiguous_masked(self): + # # Test that a contiguity check will raise a warning (not an error) for + # # discontiguous bounds but appropriately masked data. + # cube = self.latlon_2d + # coord = cube.coord('longitude') + # msg = 'The bounds of the longitude coordinate are not contiguous. ' \ + # 'However, data is masked where the discontiguity occurs so ' \ + # 'plotting anyway.' + # with self.assertWarnsRegexp(msg): + # iplt._check_contiguity_and_bounds(coord, cube.data) def test_discontiguous_unmasked(self): # Check that an error occurs when the contiguity check finds From a281805785e212f2f23cdb8a1ecc372978feca39 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 2 Aug 2018 12:32:50 +0100 Subject: [PATCH 22/40] Ensure Sphinx autodocs for grid_angles routines. --- lib/iris/analysis/cartography.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/lib/iris/analysis/cartography.py b/lib/iris/analysis/cartography.py index 16c99c83a8..0c19b84815 100644 --- a/lib/iris/analysis/cartography.py +++ b/lib/iris/analysis/cartography.py @@ -42,6 +42,23 @@ true_vectors_from_grid_vectors as rotate_grid_vectors) +# List of contents to control Sphinx autodocs. +# Unfortunately essential to get docs for the grid_angles functions. +__all__ = [ + 'area_weights', + 'cosine_latitude_weights', + 'get_xy_contiguous_bounded_grids', + 'get_xy_grids', + 'gridcell_angles', + 'project', + 'rotate_grid_vectors', + 'rotate_pole', + 'rotate_winds', + 'unrotate_pole', + 'wrap_lons', + 'DistanceDifferential', + 'PartialDifferential'] + # This value is used as a fall-back if the cube does not define the earth DEFAULT_SPHERICAL_EARTH_RADIUS = 6367470 # TODO: This should not be necessary, as CF is always in meters From e024f5019702013c228844da410af8bc014f30ea Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 2 Aug 2018 12:34:54 +0100 Subject: [PATCH 23/40] Codestyle fixes. --- lib/iris/analysis/cartography.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/iris/analysis/cartography.py b/lib/iris/analysis/cartography.py index 0c19b84815..dadcdedc86 100644 --- a/lib/iris/analysis/cartography.py +++ b/lib/iris/analysis/cartography.py @@ -352,7 +352,7 @@ def _quadrant_area(radian_lat_bounds, radian_lon_bounds, radius_of_earth): def area_weights(cube, normalize=False): - """ + r""" Returns an array of area weights, with the same dimensions as the cube. This is a 2D lat/lon area weights array, repeated over the non lat/lon @@ -460,7 +460,7 @@ def area_weights(cube, normalize=False): def cosine_latitude_weights(cube): - """ + r""" Returns an array of latitude weights, with the same dimensions as the cube. The weights are the cosine of latitude. @@ -970,7 +970,7 @@ def _transform_distance_vectors_tolerance_mask(src_crs, x, y, tgt_crs, def rotate_winds(u_cube, v_cube, target_cs): - """ + r""" Transform wind vectors to a different coordinate system. The input cubes contain U and V components parallel to the local X and Y From e3e82e458f68ebbc008f43eb5ce91bced038d445 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 2 Aug 2018 12:53:08 +0100 Subject: [PATCH 24/40] Whatsnew entries for 2d vector support. --- docs/iris/src/whatsnew/2.2.rst | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/docs/iris/src/whatsnew/2.2.rst b/docs/iris/src/whatsnew/2.2.rst index 97eba3dee6..4edc305025 100644 --- a/docs/iris/src/whatsnew/2.2.rst +++ b/docs/iris/src/whatsnew/2.2.rst @@ -24,3 +24,13 @@ Iris 2.2.0 Features automatically checked for contiguity. Coordinate bounds must either be contiguous, or the cube's data must be masked at the discontiguities in order to avoid plotting errors. + + The iris plot functions :func:`iris.plot.quiver` and + :func:`iris.plot.streamplot` have been added, and these also work with + 2-dimensional plot coordinates. + +.. admonition:: 2-Dimensional Grid Vectors + + The iris functions :func:`iris.analysis.cartography.gridcell_angles` and + :func:`iris.analysis.cartography.rotate_grid_vectors` have been added, + allowing you to convert gridcell-oriented vectors to true-North/East ones. From b041c50ca371dacbbfc14489b35aa9a9dca44b50 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 14 Aug 2018 15:55:29 +0100 Subject: [PATCH 25/40] Pin Dask to avoid 0.18.2 bug with masked arrays. (#3127) * Put orca_util routines in subpackage of iris.test.stock, to get existing iris.tests.unit.plot.test_2d_coords working. * Disable broken misused testcode in test_gridcell_angles. * Ditch test_gridcell_angles, none of it is functional. * Further style fixes. * Skip tests using iris-test-data, for Travis TEST_MINIMAL phases. * Codestyle fix (though this code obsolete anyway). * Fix unused imports. * Made-up test cube replaces use of iris-test-data. * Renamed keyword; improved docstring. * Change 'co' to 'coord' for clarity. * Review changes. --- lib/iris/analysis/_grid_angles.py | 6 +- lib/iris/coords.py | 2 +- lib/iris/plot.py | 2 +- .../tests/{stock.py => stock/__init__.py} | 3 +- lib/iris/tests/stock/_stock_2d_latlons.py | 346 ++++++++++++++++ lib/iris/tests/test_coding_standards.py | 2 +- .../cartography/test_gridcell_angles.py | 392 ------------------ lib/iris/tests/unit/plot/test_2d_coords.py | 46 +- 8 files changed, 382 insertions(+), 417 deletions(-) rename lib/iris/tests/{stock.py => stock/__init__.py} (99%) create mode 100644 lib/iris/tests/stock/_stock_2d_latlons.py delete 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 1746834a65..ffe6d51c7f 100644 --- a/lib/iris/analysis/_grid_angles.py +++ b/lib/iris/analysis/_grid_angles.py @@ -92,7 +92,7 @@ def _angle(p, q, r): if old_style: mid_lons = np.deg2rad(q[0]) - pr = _3d_xyz_from_latlon(r[0], r[1]) - _3d_xyz_from_latlon(p[0], p[1]) + pr = _3d_xyz_from_latlon(r[0], r[1]) - _3d_xyz_from_latlon(p[0], p[1]) pr_norm = np.sqrt(np.sum(pr**2, axis=0)) pr_top = pr[1] * np.cos(mid_lons) - pr[0] * np.sin(mid_lons) @@ -124,7 +124,7 @@ def _angle(p, q, r): lmb_hatvec_y, lmb_hatvec_z)]) - pr = _3d_xyz_from_latlon(r[0], r[1]) - _3d_xyz_from_latlon(p[0], p[1]) + pr = _3d_xyz_from_latlon(r[0], r[1]) - _3d_xyz_from_latlon(p[0], p[1]) # Dot products to form true-northward / true-eastward projections. pr_cmpt_e = np.sum(pr * lmb_hatvec, axis=0) @@ -141,7 +141,7 @@ def _angle(p, q, r): rtol = 1.e-3 check = np.allclose(mag_rot, mag_orig, rtol=rtol) if not check: - print (mag_rot, mag_orig) + print(mag_rot, mag_orig) assert np.allclose(mag_rot, mag_orig, rtol=rtol) return psi diff --git a/lib/iris/coords.py b/lib/iris/coords.py index 8f1e48fd18..ff2f170a2f 100644 --- a/lib/iris/coords.py +++ b/lib/iris/coords.py @@ -1071,7 +1071,7 @@ def is_contiguous(self, rtol=1e-05, atol=1e-08): rtol=rtol, atol=atol) elif self.ndim == 2: contiguous, _, _ = _discontiguity_in_2d_bounds(self.bounds, - abs_tol=atol) + abs_tol=atol) else: contiguous = False return contiguous diff --git a/lib/iris/plot.py b/lib/iris/plot.py index 8f11c8d116..7bf8ea07f3 100644 --- a/lib/iris/plot.py +++ b/lib/iris/plot.py @@ -334,7 +334,7 @@ def _draw_2d_from_bounds(draw_method_name, cube, *args, **kwargs): plot_defn = _get_plot_defn(cube, mode, ndims=2) twodim_contig_atol = kwargs.pop('two_dim_coord_contiguity_atol', - 1e-4) + 1e-4) for coord in plot_defn.coords: if hasattr(coord, 'has_bounds'): if coord.ndim == 2 and coord.has_bounds(): diff --git a/lib/iris/tests/stock.py b/lib/iris/tests/stock/__init__.py similarity index 99% rename from lib/iris/tests/stock.py rename to lib/iris/tests/stock/__init__.py index 1b55510d2b..5965d5a208 100644 --- a/lib/iris/tests/stock.py +++ b/lib/iris/tests/stock/__init__.py @@ -36,7 +36,8 @@ from iris.coords import DimCoord, AuxCoord import iris.tests as tests from iris.coord_systems import GeogCS, RotatedGeogCS - +from ._stock_2d_latlons import (sample_2d_latlons, + make_bounds_discontiguous_at_point) def lat_lon_cube(): """ diff --git a/lib/iris/tests/stock/_stock_2d_latlons.py b/lib/iris/tests/stock/_stock_2d_latlons.py new file mode 100644 index 0000000000..75b62f31ba --- /dev/null +++ b/lib/iris/tests/stock/_stock_2d_latlons.py @@ -0,0 +1,346 @@ +# (C) British Crown Copyright 2018, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +""" +Extra stock routines for making and manipulating cubes with 2d coordinates, +to mimic ocean grid data. + +""" +from __future__ import (absolute_import, division, print_function) +from six.moves import (filter, input, map, range, zip) # noqa + +import numpy as np +import numpy.ma as ma + +from iris.analysis.cartography import unrotate_pole +from iris.cube import Cube +from iris.coords import AuxCoord, DimCoord +from iris.coord_systems import RotatedGeogCS + + +def expand_1d_x_and_y_bounds_to_2d_xy(x_bounds_1d, y_bounds_1d): + """ + Convert bounds for separate 1-D X and Y coords into bounds for the + equivalent 2D coordinates. + + The output arrays have 4 points per cell, for 4 'corners' of a gridcell, + in the usual anticlockwise order + (bottom-left, bottom-right, top-right, top-left). + + If 1-dimensional X and Y coords have shapes nx and ny, then their bounds + have shapes (nx, 2) and (ny, 2). + The equivalent 2d coordinates would have values which are a "meshgrid" of + the original 1-D points, both having the shape (ny, ny). + The outputs are 2d bounds arrays suitable for such 2d coordinates. + + Args: + + * x_bounds_1d, y_bounds_1d : (array) + Coordinate bounds arrays, with shapes (nx, 2) and (ny, 2). + + Result: + + * x_bounds_2d, y_bounds_2d : (array) + Expanded 2d bounds arrays, both of shape (ny, nx, 4). + + """ + shapes = [bds.shape for bds in (x_bounds_1d, y_bounds_1d)] + for shape in shapes: + if len(shape) != 2 or shape[1] != 2: + msg = ('One-dimensional bounds arrays must have shapes (ny, 2) ' + 'and (nx, 2). Got {} and {}.') + raise ValueError(msg.format(*shapes)) + + # Construct output arrays, which are both (ny, nx, 4). + nx, ny = [shape[0] for shape in shapes] + bds_2d_x = np.zeros((ny, nx, 4)) + bds_2d_y = bds_2d_x.copy() + + # Expand x bounds to 2D array (ny, nx, 4) : the same over all 'Y'. + # Bottom left+right corners are the original 1-D x bounds. + bds_2d_x[:, :, 0] = x_bounds_1d[:, 0].reshape((1, nx)) + bds_2d_x[:, :, 1] = x_bounds_1d[:, 1].reshape((1, nx)) + # Top left+right corners are the same as bottom left+right. + bds_2d_x[:, :, 2] = bds_2d_x[:, :, 1].copy() + bds_2d_x[:, :, 3] = bds_2d_x[:, :, 0].copy() + + # Expand y bounds to 2D array (ny, nx, 4) : the same over all 'X'. + # Left-hand lower+upper corners are the original 1-D y bounds. + bds_2d_y[:, :, 0] = y_bounds_1d[:, 0].reshape((ny, 1)) + bds_2d_y[:, :, 3] = y_bounds_1d[:, 1].reshape((ny, 1)) + # Right-hand lower+upper corners are the same as the left-hand ones. + bds_2d_y[:, :, 1] = bds_2d_y[:, :, 0].copy() + bds_2d_y[:, :, 2] = bds_2d_y[:, :, 3].copy() + + return bds_2d_x, bds_2d_y + + +def grid_coords_2d_from_1d(x_coord_1d, y_coord_1d): + """ + Calculate a pair of 2d X+Y coordinates from 1d ones, in a "meshgrid" style. + If the inputs are bounded, the outputs have 4 points per bounds in the + usual way, i.e. points 0,1,2,3 are the gridcell corners anticlockwise from + the bottom left. + + """ + for coord in (x_coord_1d, y_coord_1d): + if coord.ndim != 1: + msg = ('Input coords must be one-dimensional. ' + 'Coordinate "{}" has shape {}.') + raise ValueError(msg.format(coord.name(), coord.shape)) + + # Calculate centre-points as a mesh of the 2 inputs. + pts_2d_x, pts_2d_y = np.meshgrid(x_coord_1d.points, y_coord_1d.points) + if not x_coord_1d.has_bounds() or not y_coord_1d.has_bounds(): + bds_2d_x = None + bds_2d_y = None + else: + bds_2d_x, bds_2d_y = expand_1d_x_and_y_bounds_to_2d_xy( + x_coord_1d.bounds, y_coord_1d.bounds) + + # Make two new coords + return them. + result = [] + for pts, bds, name in zip((pts_2d_x, pts_2d_y), + (bds_2d_x, bds_2d_y), + ('longitude', 'latitude')): + coord = AuxCoord(pts, bounds=bds, standard_name=name, units='degrees') + result.append(coord) + + return result + + +def sample_2d_latlons(regional=False, rotated=False, transformed=False): + """ + Construct small 2d cubes with 2d X and Y coordinates. + + This makes cubes with 'expanded' coordinates (4 bounds per cell), analagous + to ORCA data. + The coordinates are always geographical, so either it has a coord system + or they are "true" lats + lons. + ( At present, they are always latitudes and longitudes, but maybe in a + rotated system. ) + The results always have fully contiguous bounds. + + Kwargs: + * regional (bool): + If False (default), results cover the whole globe, and there is + implicit connectivity between rhs + lhs of the array. + If True, coverage is regional and edges do not connect. + * rotated (bool): + If False, X and Y coordinates are true-latitudes and longitudes, with + an implicit coordinate system (i.e. None). + If True, the X and Y coordinates are lats+lons in a selected + rotated-latlon coordinate system. + * transformed (bool): + Build coords from rotated coords as for 'rotated', but then replace + their values with the equivalent "true" lats + lons, and no + coord-system (defaults to true-latlon). + In this case, the X and Y coords are no longer 'meshgrid' style, + i.e. the points + bounds values vary in *both* dimensions. + + .. note:: + + 'transformed' is an alternative to 'rotated' : when 'transformed' is + set, then 'rotated' has no effect. + + .. Some sample results printouts :: + + >>> print(sample_2d_latlons()) + test_data / (unknown) (-- : 5; -- : 6) + Auxiliary coordinates: + latitude x x + longitude x x + >>> + >>> print(sample_2d_latlons().coord(axis='x')[0, :2]) + AuxCoord(array([ 37.5 , 93.75]), + bounds=array([[ 0. , 65.625, 65.625, 0. ], + [ 65.625, 121.875, 121.875, 65.625]]), + standard_name='longitude', units=Unit('degrees')) + >>> print(np.round(sample_2d_latlons().coord(axis='x').points, 3)) + [[ 37.5 93.75 150. 206.25 262.5 318.75] + [ 37.5 93.75 150. 206.25 262.5 318.75] + [ 37.5 93.75 150. 206.25 262.5 318.75] + [ 37.5 93.75 150. 206.25 262.5 318.75] + [ 37.5 93.75 150. 206.25 262.5 318.75]] + >>> print(np.round(sample_2d_latlons().coord(axis='y').points, 3)) + [[-85. -85. -85. -85. -85. -85. ] + [-47.5 -47.5 -47.5 -47.5 -47.5 -47.5] + [-10. -10. -10. -10. -10. -10. ] + [ 27.5 27.5 27.5 27.5 27.5 27.5] + [ 65. 65. 65. 65. 65. 65. ]] + + + >>> print(np.round( + sample_2d_latlons(rotated=True).coord(axis='x').points, 3)) + [[ 37.5 93.75 150. 206.25 262.5 318.75] + [ 37.5 93.75 150. 206.25 262.5 318.75] + [ 37.5 93.75 150. 206.25 262.5 318.75] + [ 37.5 93.75 150. 206.25 262.5 318.75] + [ 37.5 93.75 150. 206.25 262.5 318.75]] + >>> print(sample_2d_latlons(rotated=True).coord(axis='y').coord_system) + RotatedGeogCS(75.0, 120.0) + + + >>> print( + sample_2d_latlons(transformed=True).coord(axis='y').coord_system) + None + >>> print(np.round( + sample_2d_latlons(transformed=True).coord(axis='x').points, 3)) + [[ -50.718 -40.983 -46.74 -71.938 -79.293 -70.146] + [ -29.867 17.606 77.936 157.145 -141.037 -93.172] + [ -23.139 31.007 87.699 148.322 -154.639 -100.505] + [ -16.054 41.218 92.761 143.837 -164.738 -108.105] + [ 10.86 61.78 100.236 137.285 175.511 -135.446]] + >>> print(np.round( + sample_2d_latlons(transformed=True).coord(axis='y').points, 3)) + [[-70.796 -74.52 -79.048 -79.26 -74.839 -70.96 ] + [-34.99 -46.352 -59.721 -60.34 -47.305 -35.499] + [ 1.976 -10.626 -22.859 -23.349 -11.595 1.37 ] + [ 38.914 25.531 14.312 13.893 24.585 38.215] + [ 74.197 60.258 51.325 51.016 59.446 73.268]] + >>> + + """ + def sample_cube(xargs, yargs): + # Make a test cube with given latitude + longitude coordinates. + # xargs/yargs are args for np.linspace (start, stop, N), to make the X + # and Y coordinate points. + x0, x1, nx = xargs + y0, y1, ny = yargs + # Data has cycling values, staggered a bit in successive rows. + data = np.zeros((ny, nx)) + data.flat[:] = np.arange(ny * nx) % (nx + 2) + # Build a 2d cube with longitude + latitude coordinates. + cube = Cube(data, long_name='test_data') + x_pts = np.linspace(x0, x1, nx, endpoint=True) + y_pts = np.linspace(y0, y1, ny, endpoint=True) + co_x = DimCoord(x_pts, standard_name='longitude', units='degrees') + co_y = DimCoord(y_pts, standard_name='latitude', units='degrees') + cube.add_dim_coord(co_y, 0) + cube.add_dim_coord(co_x, 1) + return cube + + if regional: + # Extract small region. + cube = sample_cube(xargs=(150., 243.75, 6), yargs=(-10., 40., 5)) + + # Make contiguous bounds. + for ax in ('x', 'y'): + cube.coord(axis=ax).guess_bounds() + else: + # Global data, but drastically reduced resolution. + cube = sample_cube(xargs=(37.5, 318.75, 6), yargs=(-85., 65., 5)) + + # Patch bounds to ensure it is still contiguous + global. + for name in ('longitude', 'latitude'): + coord = cube.coord(name) + coord.guess_bounds() + bds = coord.bounds.copy() + # Make bounds global, by fixing lowest and uppermost values. + if name == 'longitude': + bds[0, 0] = 0.0 + bds[-1, 1] = 360.0 + else: + bds[0, 0] = -90.0 + bds[-1, 1] = 90.0 + coord.bounds = bds + + # Get 1d coordinate points + bounds + calculate 2d equivalents. + co_1d_x, co_1d_y = [cube.coord(axis=ax).copy() for ax in ('x', 'y')] + co_2d_x, co_2d_y = grid_coords_2d_from_1d(co_1d_x, co_1d_y) + + # Remove the old grid coords. + for coord in (co_1d_x, co_1d_y): + cube.remove_coord(coord) + + # Add the new grid coords. + for coord in (co_2d_x, co_2d_y): + cube.add_aux_coord(coord, (0, 1)) + + if transformed or rotated: + # Take the lats + lons as being in a rotated coord system. + pole_lat, pole_lon = 75.0, 120.0 + if transformed: + # Reproject coordinate values from rotated to true lat-lons. + co_x, co_y = [cube.coord(axis=ax) for ax in ('x', 'y')] + # Unrotate points. + lons, lats = co_x.points, co_y.points + lons, lats = unrotate_pole(lons, lats, pole_lon, pole_lat) + co_x.points, co_y.points = lons, lats + # Unrotate bounds. + lons, lats = co_x.bounds, co_y.bounds + # Note: save the shape, flatten + then re-apply the shape, because + # "unrotate_pole" uses "cartopy.crs.CRS.transform_points", which + # only works on arrays of 1 or 2 dimensions. + shape = lons.shape + lons, lats = unrotate_pole(lons.flatten(), lats.flatten(), + pole_lon, pole_lat) + co_x.bounds, co_y.bounds = lons.reshape(shape), lats.reshape(shape) + else: + # "Just" rotate operation : add a coord-system to each coord. + cs = RotatedGeogCS(pole_lat, pole_lon) + for coord in cube.coords(): + coord.coord_system = cs + + return cube + + +def make_bounds_discontiguous_at_point(cube, at_iy, at_ix): + """ + Meddle with the XY grid bounds of a cube to make the grid discontiguous. + + Changes the points and bounds of a single gridcell, so that it becomes + discontinuous with the next gridcell to its right. + Also masks the cube data at the given point. + + The cube must have bounded 2d 'x' and 'y' coordinates. + + TODO: add a switch to make a discontiguity in the *y*-direction instead ? + + """ + x_coord = cube.coord(axis='x') + y_coord = cube.coord(axis='y') + assert x_coord.shape == y_coord.shape + assert (coord.bounds.ndim == 3 and coord.shape[-1] == 4 + for coord in (x_coord, y_coord)) + + # For both X and Y coord, move points + bounds to create a discontinuity. + def adjust_coord(coord): + pts, bds = coord.points, coord.bounds + # Fetch the 4 bounds (bottom-left, bottom-right, top-right, top-left) + bds_bl, bds_br, bds_tr, bds_tl = bds[at_iy, at_ix] + # Make a discontinuity "at" (iy, ix), by moving the right-hand edge of + # the cell to the midpoint of the existing left+right bounds. + new_bds_br = 0.5 * (bds_bl + bds_br) + new_bds_tr = 0.5 * (bds_tl + bds_tr) + bds_br, bds_tr = new_bds_br, new_bds_tr + bds[at_iy, at_ix] = [bds_bl, bds_br, bds_tr, bds_tl] + # Also reset the cell midpoint to the middle of the 4 new corners, + # in case having a midpoint outside the corners might cause a problem. + new_pt = 0.25 * sum([bds_bl, bds_br, bds_tr, bds_tl]) + pts[at_iy, at_ix] = new_pt + # Write back the coord points+bounds (can only assign whole arrays). + coord.points, coord.bounds = pts, bds + + adjust_coord(x_coord) + adjust_coord(y_coord) + # Also mask the relevant data point. + data = cube.data # N.B. fetch all the data. + if not ma.isMaskedArray(data): + # Promote to masked array, to avoid converting mask to NaN. + data = ma.masked_array(data) + data[at_iy, at_ix] = ma.masked + cube.data = data diff --git a/lib/iris/tests/test_coding_standards.py b/lib/iris/tests/test_coding_standards.py index 39c73e8285..c24ac46cb3 100644 --- a/lib/iris/tests/test_coding_standards.py +++ b/lib/iris/tests/test_coding_standards.py @@ -93,7 +93,7 @@ class StandardReportWithExclusions(pep8.StandardReport): '*/iris/io/format_picker.py', '*/iris/tests/__init__.py', '*/iris/tests/pp.py', - '*/iris/tests/stock.py', + '*/iris/tests/stock/__init__.py', '*/iris/tests/system_test.py', '*/iris/tests/test_analysis.py', '*/iris/tests/test_analysis_calculus.py', diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py deleted file mode 100644 index de00d2bc76..0000000000 --- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py +++ /dev/null @@ -1,392 +0,0 @@ -# (C) British Crown Copyright 2018, Met Office -# -# This file is part of Iris. -# -# Iris is free software: you can redistribute it and/or modify it under -# the terms of the GNU Lesser General Public License as published by the -# Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# Iris is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Lesser General Public License for more details. -# -# You should have received a copy of the GNU Lesser General Public License -# along with Iris. If not, see . -""" -Unit tests for the function -:func:`iris.analysis.cartography.gridcell_angles`. - -""" -from __future__ import (absolute_import, division, print_function) -from six.moves import (filter, input, map, range, zip) # noqa - -# Import iris.tests first so that some things can be initialised before -# importing anything else. -import iris.tests as tests - -import numpy as np -import numpy.ma as ma - -import cartopy.crs as ccrs -from iris.cube import Cube -from iris.coords import DimCoord, AuxCoord -import iris.coord_systems -from iris.analysis.cartography import unrotate_pole - -from iris.analysis.cartography import (gridcell_angles, - rotate_grid_vectors) - -import matplotlib.pyplot as plt -from orca_utils.plot_testing.blockplot_from_bounds import blockplot_2dll - - -def _rotated_grid_sample(pole_lat=15, pole_lon=-180, - lon_bounds=np.linspace(-30, 30, 6, endpoint=True), - lat_bounds=np.linspace(-30, 30, 6, endpoint=True)): - # Calculate *true* lat_bounds+lon_bounds for the rotated grid. - lon_bounds = np.array(lon_bounds, dtype=float) - lat_bounds = np.array(lat_bounds, dtype=float) - # Construct centrepoints. - lons = 0.5 * (lon_bounds[:-1] + lon_bounds[1:]) - lats = 0.5 * (lat_bounds[:-1] + lat_bounds[1:]) - # Convert all to full 2d arrays. - lon_bounds, lat_bounds = np.meshgrid(lon_bounds, lat_bounds) - lons, lats = np.meshgrid(lons, lats) - # Calculate true lats+lons for all points. - lons_true_bds, lats_true_bds = unrotate_pole(lon_bounds, lat_bounds, - pole_lon, pole_lat) - lons_true, lats_true = unrotate_pole(lons, lats, pole_lon, pole_lat) - # Make the 'unified' bounds into contiguous (ny, nx, 4) arrays. - def expand_unified_bds(bds): - ny, nx = bds.shape - bds_4 = np.zeros((ny - 1, nx - 1, 4)) - bds_4[:, :, 0] = bds[:-1, :-1] - bds_4[:, :, 1] = bds[:-1, 1:] - bds_4[:, :, 2] = bds[1:, 1:] - bds_4[:, :, 3] = bds[1:, :-1] - return bds_4 - - lon_true_bds4, lat_true_bds4 = (expand_unified_bds(bds) - for bds in (lons_true_bds, lats_true_bds)) - # Make these into a 2d-latlon grid for a cube - cube = Cube(np.zeros(lon_true_bds4.shape[:-1])) - co_x = AuxCoord(lons_true, bounds=lon_true_bds4, - standard_name='longitude', units='degrees') - co_y = AuxCoord(lats_true, bounds=lat_true_bds4, - standard_name='latitude', units='degrees') - cube.add_aux_coord(co_x, (0, 1)) - cube.add_aux_coord(co_y, (0, 1)) - return cube - - -class TestGridcellAngles(tests.IrisTest): - def _singlecell_30deg_cube(self, x0=90., y0=0., dx=20., dy=10.): - x_pts = np.array([[x0]]) - y_pts = np.array([[y0]]) - x_bds = x0 + dx * np.array([[[-1., 1, 0.5, -1.5]]]) -# self.assertArrayAllClose(x_bds, np.array([[[70., 110, 100, 60]]])) - y_bds = y0 + dy * np.array([[[-1., 1, 3, 1]]]) -# self.assertArrayAllClose(y_bds, np.array([[[-10., 10, 30, 10]]])) - co_x = AuxCoord(points=x_pts, bounds=x_bds, - standard_name='longitude', units='degrees') - co_y = AuxCoord(points=y_pts, bounds=y_bds, - standard_name='latitude', units='degrees') - cube = Cube(np.zeros((1, 1))) - cube.add_aux_coord(co_x, (0, 1)) - cube.add_aux_coord(co_y, (0, 1)) - return cube - - def _singlecell_diamond_cube(self, x0=90., y0=0., dy=10., dx_eq=None): - if dx_eq is None: - dx_eq = dy - x_pts = np.array([[x0]]) - y_pts = np.array([[y0]]) - dx = dx_eq / np.cos(np.deg2rad(y0)) - x_bds = np.array([[[x0, x0 + dx, x0, x0 - dx]]]) - y_bds = np.array([[[y0 - dy, y0, y0 + dy, y0]]]) - co_x = AuxCoord(points=x_pts, bounds=x_bds, - standard_name='longitude', units='degrees') - co_y = AuxCoord(points=y_pts, bounds=y_bds, - standard_name='latitude', units='degrees') - cube = Cube(np.zeros((1, 1))) - cube.add_aux_coord(co_x, (0, 1)) - cube.add_aux_coord(co_y, (0, 1)) - return cube - - def test_single_cell_equatorial(self): - plt.switch_backend('tkagg') - plt.figure(figsize=(10,10)) - ax = plt.axes(projection=ccrs.Mercator()) -# ax = plt.axes(projection=ccrs.NorthPolarStereo()) -# ax = plt.axes(projection=ccrs.Orthographic(central_longitude=90., -# central_latitude=30.)) - - lon0 = 90.0 - dy = 2.0 - y_0, y_n, ny = -80, 80, 9 - angles = [] - for lat in np.linspace(y_0, y_n, ny): - cube = self._singlecell_diamond_cube(x0=lon0, y0=lat, dy=dy) - angles_cube = gridcell_angles(cube, -# cell_angle_boundpoints='mid-lhs, mid-rhs') - cell_angle_boundpoints='lower-left, lower-right') - tmp_cube = angles_cube.copy() - tmp_cube.convert_units('degrees') - print('') - print(lat) - co_x, co_y = (cube.coord(axis=ax) for ax in ('x', 'y')) - print() - print(' at : {}, {}'.format(co_x.points[0, 0], co_y.points[0, 0])) - print(' x-bds:') - print(co_x.bounds) - print(' y-bds:') - print(co_y.bounds) - angle = tmp_cube.data[0, 0] - angles.append(angle) - print(angle) - blockplot_2dll(cube) - - ax.coastlines() - ax.set_global() - - # Plot constant NEly (45deg) arrows. - xx = np.array([lon0] * ny) - yy = np.linspace(y_0, y_n, ny) - dy - uu = np.array([1.0] * ny) - plt.quiver(xx, yy, - uu, np.cos(np.deg2rad(yy)), - zorder=2, color='red', - scale_units='xy', - transform=ccrs.PlateCarree()) - - # Also plot returned angles. - angles_arr_rad = np.deg2rad(angles) - u_arr = uu * np.cos(angles_arr_rad) - v_arr = uu * np.sin(angles_arr_rad) * np.cos(np.deg2rad(yy)) - - plt.quiver(xx, yy, - u_arr, - v_arr, - zorder=2, color='magenta', - scale_units='xy', - width=0.005, - scale=0.2e-6, -# width=0.5, - transform=ccrs.PlateCarree()) - - plt.show() - - - def test_values(self): - # Construct a rotated-pole grid and check angle calculation. - testcube = _rotated_grid_sample() - - cell_angle_boundpoints = 'mid-lhs, mid-rhs' -# cell_angle_boundpoints = 'lower-left, lower-right' -# cell_angle_boundpoints = 'garble' - angles_cube = gridcell_angles( - testcube, - cell_angle_boundpoints=cell_angle_boundpoints) - angles_cube.convert_units('radians') - - # testing phase... - print(np.rad2deg(angles_cube.data)) - - import matplotlib.pyplot as plt - plt.switch_backend('tkagg') - -# plot_map = 'north_polar_stereographic' -# plot_map = 'plate_carree' -# plot_map = 'mercator' - plot_map = 'north_polar_orthographic' - if plot_map == 'plate_carree': - scale = 0.1 - map_proj = ccrs.PlateCarree() - elif plot_map == 'mercator': - scale = 3.0e-6 - map_proj = ccrs.Mercator() - map_proj._threshold *= 0.01 - elif plot_map == 'north_polar_orthographic': - scale = 3.0e-6 - map_proj = ccrs.Orthographic(central_longitude=0.0, - central_latitude=90.0,) - map_proj._threshold *= 0.01 - elif plot_map == 'north_polar_stereographic': - scale = 3.0e-6 - map_proj = ccrs.NorthPolarStereo() - else: - assert 0 - - ax = plt.axes(projection=map_proj) - data_proj = ccrs.PlateCarree() - - deg_scale = 10.0 - -# angles = 'uv' - angles = 'xy' - - ax.coastlines() - ax.gridlines() - for i_bnd in range(4): - color = ['black', 'red', 'blue', 'magenta'][i_bnd] - plt.plot(testcube.coord('longitude').bounds[..., i_bnd], - testcube.coord('latitude').bounds[..., i_bnd], - '+', markersize=10., markeredgewidth=2., - markerfacecolor=color, markeredgecolor=color, - transform=data_proj) - - - # Show plain 0,1 + 1,0 (PlateCarree) vectors unrotated at the given points. - pts_shape = testcube.coord('longitude').shape - ny, nx = pts_shape - u0 = np.ones(pts_shape) - v0 = np.zeros(pts_shape) - u1 = v0.copy() - v1 = u0.copy() - - x0s = testcube.coord('longitude').points - y0s = testcube.coord('latitude').points - yscale = np.cos(np.deg2rad(y0s)) - plt.quiver(x0s, y0s, u0, v0 * yscale, - color='blue', width=0.005, - headwidth=2., # headlength=1.0, headaxislength=0.7, - angles=angles, - scale_units='xy', scale=scale, - transform=data_proj) - plt.quiver(x0s, y0s, u1, v1 * yscale, - color='red', width=0.005, - headwidth=2., # headlength=1.0, headaxislength=0.7, - angles=angles, - scale_units='xy', scale=scale, - transform=data_proj) - - # Add 45deg arrows (NEly), still on a PlateCarree map. - plt.quiver(x0s, y0s, v1, v1 * yscale, - color='green', width=0.005, - headwidth=2., # headlength=1.0, headaxislength=0.7, - angles=angles, - scale_units='xy', scale=scale, - transform=data_proj) - - - - # - # Repeat the above plotting short lines INSTEAD of quiver. - # - u0d = x0s + deg_scale * u0 - v0d = y0s + deg_scale * v0 - u1d = x0s + deg_scale * u1 - v1d = y0s + deg_scale * v1 - u2d = x0s + deg_scale * u0 - v2d = y0s + deg_scale * v1 - for iy in range(ny): - for ix in range(nx): - plt.plot([x0s[iy, ix], u0d[iy, ix]], - [y0s[iy, ix], v0d[iy, ix]], - ':', color='blue', linewidth=0.5, - transform=data_proj) - plt.plot([x0s[iy, ix], u1d[iy, ix]], - [y0s[iy, ix], v1d[iy, ix]], - ':', color='red', linewidth=0.5, - transform=data_proj) - plt.plot([x0s[iy, ix], u2d[iy, ix]], - [y0s[iy, ix], v2d[iy, ix]], - ':', color='green', linewidth=0.5, - transform=data_proj) - - - # Overplot BL-BR and BL-TL lines from the cell bounds. - co_lon, co_lat = [testcube.coord(name).copy() - for name in ('longitude', 'latitude')] - for co in (co_lon, co_lat): - co.convert_units('degrees') - lon_bds, lat_bds = [co.bounds for co in (co_lon, co_lat)] -# ny, nx = lon_bds.shape[:-1] - for iy in range(ny): - for ix in range(nx): - x0, y0 = lon_bds[iy, ix, 0], lat_bds[iy, ix, 0] - x1, y1 = lon_bds[iy, ix, 1], lat_bds[iy, ix, 1] - x2, y2 = lon_bds[iy, ix, 3], lat_bds[iy, ix, 3] - plt.plot([x0, x1], [y0, y1], 'x-', - color='orange', - transform=data_proj) - plt.plot([x0, x2], [y0, y2], 'x-', - color='orange', linestyle='--', - transform=data_proj) - - # Plot U0, rotated by cell angles, also at cell bottom-lefts. - u0_cube, u1_cube, v0_cube, v1_cube = [testcube.copy(data=aa) - for aa in (u0, v0, u1, v1)] - u0r_cube, v0r_cube = rotate_grid_vectors( - u0_cube, v0_cube, grid_angles_cube=angles_cube) - u0r, v0r = [cube.data for cube in (u0r_cube, v0r_cube)] - - xbl, ybl = lon_bds[..., 0], lat_bds[..., 0] - # - # Replace quiver here with delta-based lineplot - # - urd = xbl + deg_scale * u0r - vrd = ybl + deg_scale * v0r * yscale - for iy in range(ny): - for ix in range(nx): - plt.plot([xbl[iy, ix], urd[iy, ix]], - [ybl[iy, ix], vrd[iy, ix]], - ':', color='magenta', linewidth=2.5, - transform=data_proj) - # Show this is the SAME as lineplot - plt.quiver(xbl, ybl, u0r, v0r * yscale, - color='magenta', width=0.01, - headwidth=1.2, # headlength=1.0, headaxislength=0.7, - angles=angles, - scale_units='xy', scale=scale, - transform=data_proj) - - plt.suptitle('angles from "{}"'.format(cell_angle_boundpoints)) - -# # Also draw small lines pointing at the correct (TRUE, not ) angle. -# ny, nx = x0s.shape -# size_degrees = 1.0 -# angles = angles_cube.copy() -# angles.convert_units('radians') -# angles = angles.data -# lats = testcube.coord('latitude').copy() -# lats.convert_units('radians') -# lats = lats.points -# dxs = size_degrees * u0.copy() #* np.cos(angles) -# dys = size_degrees * u0.copy() # / np.sqrt(np.cos(lats)) -# x1s = x0s + dxs -# y1s = y0s + dys -## for iy in range(ny): -## for ix in range(nx): -## plt.plot([x0s[iy, ix], x1s[iy, ix]], -## [y0s[iy, ix], y1s[iy, ix]], -## 'o-', markersize=4., markeredgewidth=0., -## color='green', # scale_units='xy', scale=scale, -## transform=data_proj) -# plt.quiver(x0s, y0s, dxs, dys, -# color='green', linewidth=0.2, -# angles=angles, -# scale_units='xy', scale=scale * 0.6, -# transform=data_proj) - - - - ax.set_global() - plt.show() - - angles_cube.convert_units('degrees') - - self.assertArrayAllClose( - angles_cube.data, - [[33.421, 17.928, 0., -17.928, -33.421], - [41.981, 24.069, 0., -24.069, -41.981], - [56.624, 37.809, 0., -37.809, -56.624], - [79.940, 74.227, 0., -74.227, -79.940], - [107.313, 126.361, -180., -126.361, -107.313]], - atol=0.002) - - -if __name__ == "__main__": - tests.main() diff --git a/lib/iris/tests/unit/plot/test_2d_coords.py b/lib/iris/tests/unit/plot/test_2d_coords.py index bb27341baa..7397921861 100644 --- a/lib/iris/tests/unit/plot/test_2d_coords.py +++ b/lib/iris/tests/unit/plot/test_2d_coords.py @@ -22,20 +22,25 @@ # Import iris.tests first so that some things can be initialised before # importing anything else. import iris.tests as tests -import iris.coords as coords import numpy.ma as ma +import iris.coords as coords from iris.tests.stock import simple_2d_w_multidim_coords as cube_2dcoords from iris.tests.stock import simple_3d_w_multidim_coords as cube3d_2dcoords -from iris.tests.stock import lat_lon_cube +from iris.tests.stock import sample_2d_latlons +from iris.tests.stock import make_bounds_discontiguous_at_point -from orca_utils.plot_testing import testdata_2d_coords as testdata if tests.MPL_AVAILABLE: import iris.plot as iplt +def full2d_global(): + return sample_2d_latlons(transformed=True) + + +@tests.skip_data class Test_2d_coords_plot_defn_bound_mode(tests.IrisTest): def setUp(self): self.multidim_cube = cube_2dcoords() @@ -43,8 +48,8 @@ def setUp(self): # latlon_2d is a cube with 2d coords, 4 bounds per point, # discontiguities in the bounds but masked data at the discontiguities. - self.latlon_2d = testdata.full2d_global() - testdata.make_bounds_discontiguous_at_point(self.latlon_2d, 2, 2) + self.latlon_2d = full2d_global() + make_bounds_discontiguous_at_point(self.latlon_2d, 2, 2) # # Take a latlon cube with 1D coords, broadcast the coords into 2D # # ones, then add ONE of them back into the cube in place of original: @@ -106,9 +111,10 @@ def test_total_span_check(self): # # TODO Find out where I can put a catch for this (if necessary) # cube = self.mixed_dims # with self.assertRaises(ValueError): - # iplt._get_plot_defn_custom_coords_picked(cube, - # ('latitude', 'longitude'), - # self.mode) + # iplt._get_plot_defn_custom_coords_picked( + # cube, + # ('latitude', 'longitude'), + # self.mode) def test_map_common_not_enough_bounds(self): # Test that a lat-lon cube with 2d coords and 2 bounds per point @@ -129,16 +135,16 @@ def test_map_common_2d(self): result = iplt._map_common('pcolor', None, self.mode, cube, plot_defn) self.assertTrue(result) - # def test_discontiguous_masked(self): - # # Test that a contiguity check will raise a warning (not an error) for - # # discontiguous bounds but appropriately masked data. - # cube = self.latlon_2d - # coord = cube.coord('longitude') - # msg = 'The bounds of the longitude coordinate are not contiguous. ' \ - # 'However, data is masked where the discontiguity occurs so ' \ - # 'plotting anyway.' - # with self.assertWarnsRegexp(msg): - # iplt._check_contiguity_and_bounds(coord, cube.data) +# def test_discontiguous_masked(self): +# # Test that a contiguity check will raise a warning (not an error) for +# # discontiguous bounds but appropriately masked data. +# cube = self.latlon_2d +# coord = cube.coord('longitude') +# msg = 'The bounds of the longitude coordinate are not contiguous. ' \ +# 'However, data is masked where the discontiguity occurs so ' \ +# 'plotting anyway.' +# with self.assertWarnsRegexp(msg): +# iplt._check_contiguity_and_bounds(coord, cube.data) def test_discontiguous_unmasked(self): # Check that an error occurs when the contiguity check finds @@ -155,3 +161,7 @@ def test_draw_2d_from_bounds(self): cube = self.latlon_2d result = iplt._draw_2d_from_bounds('pcolormesh', cube) self.assertTrue(result) + + +if __name__ == '__main__': + tests.main() From ebfc4c10e6d2fdf8dc98d67c46dcdeae4fc3e729 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Fri, 27 Jul 2018 19:11:21 +0100 Subject: [PATCH 26/40] 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 27/40] 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 28/40] 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 29/40] 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 30/40] 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 31/40] 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 32/40] 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 33/40] 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 34/40] 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 35/40] 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 36/40] 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 37/40] 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 38/40] 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 39/40] 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 40/40] 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)