From 0847a4cbd9f5a8d0aabd0286dbf57b5cd2f911d3 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 20 Jun 2017 18:52:30 +0100 Subject: [PATCH 01/11] Integration tests for OceanSigmaZFactory -- lazy cases not currently working. --- .../tests/integration/aux_factory/__init__.py | 20 +++ .../aux_factory/test_OceanSigmaZFactory.py | 160 ++++++++++++++++++ lib/iris/tests/stock.py | 60 ++++++- 3 files changed, 238 insertions(+), 2 deletions(-) create mode 100644 lib/iris/tests/integration/aux_factory/__init__.py create mode 100644 lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py diff --git a/lib/iris/tests/integration/aux_factory/__init__.py b/lib/iris/tests/integration/aux_factory/__init__.py new file mode 100644 index 0000000000..f49ea330c1 --- /dev/null +++ b/lib/iris/tests/integration/aux_factory/__init__.py @@ -0,0 +1,20 @@ +# (C) British Crown Copyright 2017, 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 . +"""Integration tests for the :mod:`iris.aux_factory` package.""" + +from __future__ import (absolute_import, division, print_function) +from six.moves import (filter, input, map, range, zip) # noqa diff --git a/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py b/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py new file mode 100644 index 0000000000..b63508dfa1 --- /dev/null +++ b/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py @@ -0,0 +1,160 @@ +# (C) British Crown Copyright 2017, 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 . +""" +Integratation tests for the +`iris.aux_factory.OceanSigmaZFactory` class. + +""" + +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 itertools + +import numpy as np + +from iris._lazy_data import as_lazy_data +from iris.tests.stock import ocean_sigma_z as stock_sample_osz + + +class Test_sample(tests.IrisTest): + def setUp(self): + self.cube = stock_sample_osz() + # Snapshot result, printed with ... + # >>> np.set_printoptions(linewidth=180, + # formatter={'float':lambda x:'{:-09.3f}'.format(x)}) + # >>> print repr(coord.points) + self.basic_derived_result = np.array( + [[[[-0000.632, -0000.526, -0000.421, -0000.316], + [-0000.789, -0000.684, -0000.579, -0000.474], + [-0000.947, -0000.842, -0000.737, -0000.632]], + [[-0014.358, -0014.264, -0014.169, -0014.074], + [-0014.501, -0014.406, -0014.311, -0014.216], + [-0014.643, -0014.548, -0014.453, -0014.358]], + [[-0082.993, -0082.951, -0082.908, -0082.866], + [-0083.056, -0083.014, -0082.972, -0082.929], + [-0083.119, -0083.077, -0083.035, -0082.993]], + [[-0368.400, -0368.400, -0368.400, -0368.400], + [-0368.400, -0368.400, -0368.400, -0368.400], + [-0368.400, -0368.400, -0368.400, -0368.400]], + [[-1495.600, -1495.600, -1495.600, -1495.600], + [-1495.600, -1495.600, -1495.600, -1495.600], + [-1495.600, -1495.600, -1495.600, -1495.600]]], + + [[[-0000.842, -0000.737, -0000.632, -0000.526], + [-0001.000, -0000.895, -0000.789, -0000.684], + [-0001.158, -0001.053, -0000.947, -0000.842]], + [[-0014.548, -0014.453, -0014.358, -0014.264], + [-0014.690, -0014.595, -0014.501, -0014.406], + [-0014.832, -0014.737, -0014.643, -0014.548]], + [[-0083.077, -0083.035, -0082.993, -0082.951], + [-0083.140, -0083.098, -0083.056, -0083.014], + [-0083.203, -0083.161, -0083.119, -0083.077]], + [[-0368.400, -0368.400, -0368.400, -0368.400], + [-0368.400, -0368.400, -0368.400, -0368.400], + [-0368.400, -0368.400, -0368.400, -0368.400]], + [[-1495.600, -1495.600, -1495.600, -1495.600], + [-1495.600, -1495.600, -1495.600, -1495.600], + [-1495.600, -1495.600, -1495.600, -1495.600]]]]) + + self.derived_coord_name = \ + 'sea_surface_height_above_reference_ellipsoid' + + def _check_result(self, cube, expected_result=None, **kwargs): + if expected_result is None: + expected_result = self.basic_derived_result + coord = cube.coord(self.derived_coord_name) + result = coord.points + self.assertArrayAllClose(result, expected_result, atol=0.005, **kwargs) + + def test_basic(self): + self._check_result(self.cube) + + def _lazy_testcube(self): + cube = self.cube + for dep_name in ('depth', 'layer_depth', 'ocean_sigma_z_coordinate'): + coord = cube.coord(dep_name) + coord.points = as_lazy_data(coord.points, coord.shape) + return cube + + def test_lazy(self): + # Check same results when key coords are made lazy. + cube = self.cube + self.assertFalse( + cube.coord(self.derived_coord_name).has_lazy_points()) + cube = self._lazy_testcube() + self.assertTrue( + cube.coord(self.derived_coord_name).has_lazy_points()) + self._check_result(cube) + + def test_transpose(self): + # Check it works with all possible dimension orders. + for dims_list in itertools.permutations(range(self.cube.ndim)): + cube = self.cube.copy() + cube.transpose(dims_list) + expected = self.basic_derived_result.transpose(dims_list) + msg = 'Unexpected result when cube transposed by {}' + msg = msg.format(dims_list) + self._check_result(cube, expected, + err_msg=msg) + + def test_lazy_transpose(self): + # Check lazy calc works with all possible dimension orders. + for dims_list in itertools.permutations(range(self.cube.ndim)): + cube = self._lazy_testcube() + cube.transpose(dims_list) + expected = self.basic_derived_result.transpose(dims_list) + msg = 'Unexpected result when cube transposed by {}' + msg = msg.format(dims_list) + self._check_result(cube, expected, + err_msg=msg) + + def test_no_sigma(self): + # Check it still works when 'sigma' is removed ... + # Set all sigma points to zero + snapshot the resulting derived points. + trial_cube = self.cube.copy() + trial_cube.coord('ocean_sigma_z_coordinate').points[:] = 0.0 + expected = trial_cube.coord(self.derived_coord_name).points + # Remove sigma altogether + check the result is the same. + cube = self.cube + cube.remove_coord('ocean_sigma_z_coordinate') + self._check_result(cube, expected) + + def test_no_eta(self): + # Check it still works when 'eta' is removed ... + + # Set all sigma points to zero + snapshot the resulting derived points. + trial_cube = self.cube.copy() + trial_cube.coord('sea_surface_height').points[:] = 0.0 + expected = trial_cube.coord(self.derived_coord_name).points + # Check this has no variation between the two timepoints. + self.assertArrayAllClose(expected[0], expected[1]) + # Take first time, as no sigma --> result *has* no time dimension. + expected = expected[0] + + # Remove sigma altogether + check the result is the same. + cube = self.cube + cube.remove_coord('sea_surface_height') + self._check_result(cube, expected) + + +if __name__ == "__main__": + tests.main() diff --git a/lib/iris/tests/stock.py b/lib/iris/tests/stock.py index df63472690..25df8eda33 100644 --- a/lib/iris/tests/stock.py +++ b/lib/iris/tests/stock.py @@ -1,4 +1,4 @@ -# (C) British Crown Copyright 2010 - 2016, Met Office +# (C) British Crown Copyright 2010 - 2017, Met Office # # This file is part of Iris. # @@ -28,10 +28,12 @@ import numpy as np import numpy.ma as ma +from cf_units import Unit from iris.cube import Cube import iris.aux_factory import iris.coords import iris.coords as icoords +from iris.coords import DimCoord, AuxCoord import iris.tests as tests from iris.coord_systems import GeogCS, RotatedGeogCS @@ -482,7 +484,6 @@ def realistic_4d(): level_height_bnds, model_level_pts, sigma_pts, sigma_bnds, time_pts, \ _source_pts, forecast_period_pts, data, orography = arrays - ll_cs = RotatedGeogCS(37.5, 177.5, ellipsoid=GeogCS(6371229.0)) lat = icoords.DimCoord(lat_pts, standard_name='grid_latitude', units='degrees', @@ -558,3 +559,58 @@ def global_grib2(): path = tests.get_data_path(('GRIB', 'global_t', 'global.grib2')) cube = iris.load_cube(path) return cube + + +def ocean_sigma_z(): + """ + Return a sample cube with an + :class:`iris.aux_factory.OceanSigmaZFactory` vertical coordinate. + + This is a fairly small cube with real coordinate arrays. The coordinate + values are derived from the sample data linked at + https://github.com/SciTools/iris/pull/509#issuecomment-23565381. + + """ + co_time = DimCoord([0.0, 1.0], standard_name='time', units='') + co_lats = DimCoord([-58.1, -52.7, -46.9], + standard_name='latitude', units=Unit('degrees')) + co_lons = DimCoord([65.1, 72.9, 83.7, 96.5], + standard_name='longitude', units=Unit('degrees')) + co_ssh = AuxCoord([[[-0.63157895, -0.52631579, -0.42105263, -0.31578947], + [-0.78947368, -0.68421053, -0.57894737, -0.47368421], + [-0.94736842, -0.84210526, -0.73684211, -0.63157895]], + [[-0.84210526, -0.73684211, -0.63157895, -0.52631579], + [-1.00000000, -0.89473684, -0.78947368, -0.68421053], + [-1.15789474, -1.05263158, -0.94736842, -0.84210526]]], + standard_name=u'sea_surface_height', units=Unit('m')) + + co_sigma = AuxCoord([0., -0.1, -0.6, -1., -1.], + standard_name=u'ocean_sigma_z_coordinate', + units=Unit('1'), + attributes={'positive': 'up'}) + + co_zlay = AuxCoord([-137.2, -137.3, -137.4, -368.4, -1495.6], + long_name='layer_depth', units=Unit('m')) + co_depth = AuxCoord([[1625.7, 3921.2, 4106.4, 5243.5], + [3615.4, 4942.6, 3883.6, 4823.1], + [3263.2, 2816.3, 2741.8, 3883.6]], + standard_name=u'depth', units=Unit('m')) + co_depthc = DimCoord(137.9, long_name='depth_c', units=Unit('m')) + co_nsigma = DimCoord(3, long_name='nsigma') + + cube = Cube(np.zeros((2, 5, 3, 4))) + cube.add_dim_coord(co_time, 0) + cube.add_dim_coord(co_lats, 2) + cube.add_dim_coord(co_lons, 3) + cube.add_aux_coord(co_zlay, 1) + cube.add_aux_coord(co_sigma, 1) + cube.add_aux_coord(co_ssh, (0, 2, 3)) + cube.add_aux_coord(co_depth, (2, 3)) + cube.add_aux_coord(co_depthc) + cube.add_aux_coord(co_nsigma) + + fact = iris.aux_factory.OceanSigmaZFactory( + depth=co_depth, eta=co_ssh, depth_c=co_depthc, zlev=co_zlay, + sigma=co_sigma, nsigma=co_nsigma) + cube.add_aux_factory(fact) + return cube From 431e41ea4cf1a5247dcfb2244d0e0d752f0ccd3e Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 21 Jun 2017 13:59:45 +0100 Subject: [PATCH 02/11] Refactor osz (not yet lazy). --- lib/iris/aux_factory.py | 113 ++++++++++++++++++++++++++-------------- 1 file changed, 73 insertions(+), 40 deletions(-) diff --git a/lib/iris/aux_factory.py b/lib/iris/aux_factory.py index efd6607fc6..8729744424 100644 --- a/lib/iris/aux_factory.py +++ b/lib/iris/aux_factory.py @@ -307,31 +307,6 @@ def _remap_with_bounds(self, dependency_dims, derived_dims): nd_values_by_key[key] = nd_values return nd_values_by_key - def _shape(self, nd_values_by_key): - nd_values = sorted(nd_values_by_key.values(), - key=lambda value: value.ndim) - shape = list(nd_values.pop().shape) - for array in nd_values: - for i, size in enumerate(array.shape): - if size > 1: - # NB. If there's an inconsistency it can only come - # from a mismatch in the number of bounds (the Cube - # ensures the other dimensions must match). - # But we can't afford to raise an error now - it'd - # break Cube.derived_coords. Instead, we let the - # error happen when the derived coordinate's bounds - # are accessed. - shape[i] = size - return shape - - def _dtype(self, arrays_by_key, **other_args): - dummy_args = {} - for key, array in six.iteritems(arrays_by_key): - dummy_args[key] = np.zeros(1, dtype=array.dtype) - dummy_args.update(other_args) - dummy_data = self._derive(**dummy_args) - return dummy_data.dtype - class HybridHeightFactory(AuxCoordFactory): """ @@ -783,7 +758,7 @@ def dependencies(self): return dict(sigma=self.sigma, eta=self.eta, depth=self.depth, depth_c=self.depth_c, nsigma=self.nsigma, zlev=self.zlev) - def _derive(self, sigma, eta, depth, depth_c, + def _old_derive(self, sigma, eta, depth, depth_c, zlev, shape, nsigma_slice): # Perform the ocean sigma over z coordinate nsigma slice. if eta.ndim: @@ -800,6 +775,73 @@ def _derive(self, sigma, eta, depth, depth_c, return result + def _derive(self, sigma, eta, depth, depth_c, + zlev, nsigma, coord_dims_func): + i_levels_dim, = coord_dims_func(self.dependencies['zlev']) + allshapes = np.array( + [el.shape + for el in (sigma, eta, depth, depth_c, zlev) + if el.ndim]) + shape = list(np.max(allshapes, axis=0)) + ndims = len(shape) + # Make a slice tuple to index the first nsigma z-levels. + nsigma_slice = [slice(None)] * ndims + nsigma_slice[i_levels_dim] = slice(0, int(nsigma)) + nsigma_slice = tuple(nsigma_slice) + # Perform the ocean sigma over z coordinate nsigma slice. + if eta.ndim: + eta = eta[nsigma_slice] + if sigma.ndim: + sigma = sigma[nsigma_slice] + if depth.ndim: + depth = depth[nsigma_slice] + # Note that, this performs a point-wise minimum. + temp = eta + sigma * (np.minimum(depth_c, depth) + eta) + # Calculate the final derived result. + result = np.ones(shape, dtype=temp.dtype) * zlev + result[nsigma_slice] = temp + + return result + + def _new_derive(self, sigma, eta, depth, depth_c, + zlev, nsigma, coord_dims_func): + i_levels_dim, = coord_dims_func(self.dependencies['zlev']) + allshapes = np.array( + [el.shape + for el in (sigma, eta, depth, depth_c, zlev) + if el.ndim]) + result_shape = list(np.max(allshapes, axis=0)) + ndims = len(result_shape) + # Make a slice tuple to index the first nsigma z-levels. + z_slices_nsigma = [slice(None)] * ndims + z_slices_nsigma[i_levels_dim] = slice(0, int(nsigma)) + z_slices_nsigma = tuple(z_slices_nsigma) + # Make a slice tuple to index the remaining z-levels. + z_slices_rest = [slice(None)] * ndims + z_slices_rest[i_levels_dim] = slice(int(nsigma), None) + z_slices_rest = tuple(z_slices_rest) + # Perform the ocean sigma over z coordinate nsigma slice. + if eta.ndim: + eta = eta[z_slices_nsigma] + if sigma.ndim: + sigma = sigma[z_slices_nsigma] + if depth.ndim: + depth = depth[z_slices_nsigma] + # Note that, this performs a point-wise minimum. + nsigma_levs = eta + sigma * (np.minimum(depth_c, depth) + eta) + # Expand shape to nsigma z-levels, as it can have lower dimensionality. + ones_full_result = np.ones(result_shape, dtype=np.int16) + ones_nsigma_result = ones_full_result[z_slices_nsigma] + result_nsigma_levs = nsigma_levs * ones_nsigma_result + # Likewise, expand zlev to a full result shape. + zlev = zlev * ones_full_result + # From this, take the 'remaining' levels for the result. + result_rest_levs = zlev[z_slices_rest] + # Combine nsigma and 'rest' levels for final result. + result = np.concatenate([result_nsigma_levs, result_rest_levs], + axis=i_levels_dim) + return result + def make_coord(self, coord_dims_func): """ Returns a new :class:`iris.coords.AuxCoord` as defined by this factory. @@ -817,30 +859,21 @@ def make_coord(self, coord_dims_func): # Build the points array. nd_points_by_key = self._remap(dependency_dims, derived_dims) - points_shape = self._shape(nd_points_by_key) - - # Calculate the nsigma slice. - nsigma_slice = [slice(None)] * len(derived_dims) - dim, = dependency_dims['zlev'] - index = derived_dims.index(dim) - nsigma_slice[index] = slice(0, int(nd_points_by_key['nsigma'])) - nsigma_slice = tuple(nsigma_slice) + nsigma, = nd_points_by_key['nsigma'] points = self._derive(nd_points_by_key['sigma'], nd_points_by_key['eta'], nd_points_by_key['depth'], nd_points_by_key['depth_c'], nd_points_by_key['zlev'], - points_shape, - nsigma_slice) + nsigma, + coord_dims_func) bounds = None if self.zlev.nbounds or (self.sigma and self.sigma.nbounds): # Build the bounds array. nd_values_by_key = self._remap_with_bounds(dependency_dims, derived_dims) - bounds_shape = self._shape(nd_values_by_key) - nsigma_slice_bounds = nsigma_slice + (slice(None),) # Define the function here to obtain a closure. def calc_bounds(): @@ -867,8 +900,8 @@ def calc_bounds(): nd_values_by_key['depth'], nd_values_by_key['depth_c'], nd_values_by_key['zlev'], - bounds_shape, - nsigma_slice_bounds) + nsigma, + coord_dims_func) bounds = calc_bounds() coord = iris.coords.AuxCoord(points, From 90b257331da83cc84870e2eafb55f90c037fcdc4 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 21 Jun 2017 15:34:15 +0100 Subject: [PATCH 03/11] Add test with extra cube dims. --- .../aux_factory/test_OceanSigmaZFactory.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py b/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py index b63508dfa1..9f3c08c784 100644 --- a/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py +++ b/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py @@ -33,6 +33,7 @@ from iris._lazy_data import as_lazy_data from iris.tests.stock import ocean_sigma_z as stock_sample_osz +import iris.util class Test_sample(tests.IrisTest): @@ -127,6 +128,18 @@ def test_lazy_transpose(self): self._check_result(cube, expected, err_msg=msg) + def test_extra_dims(self): + # Insert extra cube dimensions + check it still works. + cube = self.cube + cube = iris.util.new_axis(cube) + cube = iris.util.new_axis(cube) + cube = iris.util.new_axis(cube) + # N.B. shape is now (1, 1, 1, t, z, y, x) + cube.transpose((0, 3, 1, 4, 5, 2, 6)) + # N.B. shape is now (1, t, 1, z, y, 1, x) + # Should get same original result, as derived dims are the same. + self._check_result(cube) + def test_no_sigma(self): # Check it still works when 'sigma' is removed ... # Set all sigma points to zero + snapshot the resulting derived points. From d434705f591468e9bc3141c2310e332083b26f09 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 21 Jun 2017 15:36:15 +0100 Subject: [PATCH 04/11] Fixed calculation; all working except lazy integration tests. --- lib/iris/aux_factory.py | 52 +++++++++++++++++++++++++++++++++++------ 1 file changed, 45 insertions(+), 7 deletions(-) diff --git a/lib/iris/aux_factory.py b/lib/iris/aux_factory.py index 8729744424..9ec1831e06 100644 --- a/lib/iris/aux_factory.py +++ b/lib/iris/aux_factory.py @@ -307,6 +307,23 @@ def _remap_with_bounds(self, dependency_dims, derived_dims): nd_values_by_key[key] = nd_values return nd_values_by_key +# def _shape(self, nd_values_by_key): +# nd_values = sorted(nd_values_by_key.values(), +# key=lambda value: value.ndim) +# shape = list(nd_values.pop().shape) +# for array in nd_values: +# for i, size in enumerate(array.shape): +# if size > 1: +# # NB. If there's an inconsistency it can only come +# # from a mismatch in the number of bounds (the Cube +# # ensures the other dimensions must match). +# # But we can't afford to raise an error now - it'd +# # break Cube.derived_coords. Instead, we let the +# # error happen when the derived coordinate's bounds +# # are accessed. +# shape[i] = size +# return shape + class HybridHeightFactory(AuxCoordFactory): """ @@ -758,7 +775,7 @@ def dependencies(self): return dict(sigma=self.sigma, eta=self.eta, depth=self.depth, depth_c=self.depth_c, nsigma=self.nsigma, zlev=self.zlev) - def _old_derive(self, sigma, eta, depth, depth_c, + def _original_derive(self, sigma, eta, depth, depth_c, zlev, shape, nsigma_slice): # Perform the ocean sigma over z coordinate nsigma slice. if eta.ndim: @@ -775,15 +792,27 @@ def _old_derive(self, sigma, eta, depth, depth_c, return result - def _derive(self, sigma, eta, depth, depth_c, + def _oldstyle_derive(self, sigma, eta, depth, depth_c, zlev, nsigma, coord_dims_func): - i_levels_dim, = coord_dims_func(self.dependencies['zlev']) + derived_cubedims = self.derived_dims(coord_dims_func) + i_levels_cubedim, = coord_dims_func(self.dependencies['zlev']) + i_levels_dim = i_levels_cubedim - sum( + i_dim not in derived_cubedims + for i_dim in range(i_levels_cubedim)) allshapes = np.array( [el.shape for el in (sigma, eta, depth, depth_c, zlev) if el.ndim]) shape = list(np.max(allshapes, axis=0)) ndims = len(shape) + + # N.B. it is not a problem forming 'shape': following checks out ok. + # derived_dims = self.derived_dims(coord_dims_func) + # dependency_dims = self._dependency_dims(coord_dims_func) + # nd_points_by_key = self._remap(dependency_dims, derived_dims) + # orig_shape = self._shape(nd_points_by_key) + # assert orig_shape == shape + # Make a slice tuple to index the first nsigma z-levels. nsigma_slice = [slice(None)] * ndims nsigma_slice[i_levels_dim] = slice(0, int(nsigma)) @@ -803,9 +832,18 @@ def _derive(self, sigma, eta, depth, depth_c, return result - def _new_derive(self, sigma, eta, depth, depth_c, + def _derive(self, sigma, eta, depth, depth_c, zlev, nsigma, coord_dims_func): - i_levels_dim, = coord_dims_func(self.dependencies['zlev']) + # Calculate the index of the 'z' dimension in the inputs. + # Get the cube dimension... + i_levels_cubedim, = coord_dims_func(self.dependencies['zlev']) + # ...subtract any previous dimensions not used in the calculation. + derived_cubedims = self.derived_dims(coord_dims_func) + i_levels_dim = i_levels_cubedim - sum( + i_dim not in derived_cubedims + for i_dim in range(i_levels_cubedim)) + + # Calculate the result shape as a combination of all the inputs. allshapes = np.array( [el.shape for el in (sigma, eta, depth, depth_c, zlev) @@ -829,7 +867,7 @@ def _new_derive(self, sigma, eta, depth, depth_c, depth = depth[z_slices_nsigma] # Note that, this performs a point-wise minimum. nsigma_levs = eta + sigma * (np.minimum(depth_c, depth) + eta) - # Expand shape to nsigma z-levels, as it can have lower dimensionality. + # Expand to full shape, as it may sometimes have lower dimensionality. ones_full_result = np.ones(result_shape, dtype=np.int16) ones_nsigma_result = ones_full_result[z_slices_nsigma] result_nsigma_levs = nsigma_levs * ones_nsigma_result @@ -837,7 +875,7 @@ def _new_derive(self, sigma, eta, depth, depth_c, zlev = zlev * ones_full_result # From this, take the 'remaining' levels for the result. result_rest_levs = zlev[z_slices_rest] - # Combine nsigma and 'rest' levels for final result. + # Combine nsigma and 'rest' levels for the final result. result = np.concatenate([result_nsigma_levs, result_rest_levs], axis=i_levels_dim) return result From e1b021c52741455f828308364330df09c72b5354 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 21 Jun 2017 15:49:14 +0100 Subject: [PATCH 05/11] Enable all-lazy operation; all tests working. --- lib/iris/aux_factory.py | 13 +++++---- .../aux_factory/test_OceanSigmaZFactory.py | 28 +++++++++++++------ 2 files changed, 27 insertions(+), 14 deletions(-) diff --git a/lib/iris/aux_factory.py b/lib/iris/aux_factory.py index 9ec1831e06..e3742d4ef7 100644 --- a/lib/iris/aux_factory.py +++ b/lib/iris/aux_factory.py @@ -26,6 +26,7 @@ from abc import ABCMeta, abstractmethod, abstractproperty import warnings +import dask.array as da import numpy as np from iris._cube_coord_common import CFVariableMixin @@ -188,7 +189,7 @@ def _nd_bounds(coord, dims, ndim): # Transpose to be consistent with the Cube. sorted_pairs = sorted(enumerate(dims), key=lambda pair: pair[1]) transpose_order = [pair[0] for pair in sorted_pairs] + [len(dims)] - bounds = coord.core_bounds() + bounds = coord.lazy_bounds() if dims: bounds = bounds.transpose(transpose_order) @@ -216,7 +217,7 @@ def _nd_points(coord, dims, ndim): # Transpose to be consistent with the Cube. sorted_pairs = sorted(enumerate(dims), key=lambda pair: pair[1]) transpose_order = [pair[0] for pair in sorted_pairs] - points = coord.core_points() + points = coord.lazy_points() if dims and transpose_order != list(range(len(dims))): points = points.transpose(transpose_order) @@ -866,7 +867,7 @@ def _derive(self, sigma, eta, depth, depth_c, if depth.ndim: depth = depth[z_slices_nsigma] # Note that, this performs a point-wise minimum. - nsigma_levs = eta + sigma * (np.minimum(depth_c, depth) + eta) + nsigma_levs = eta + sigma * (da.minimum(depth_c, depth) + eta) # Expand to full shape, as it may sometimes have lower dimensionality. ones_full_result = np.ones(result_shape, dtype=np.int16) ones_nsigma_result = ones_full_result[z_slices_nsigma] @@ -876,7 +877,7 @@ def _derive(self, sigma, eta, depth, depth_c, # From this, take the 'remaining' levels for the result. result_rest_levs = zlev[z_slices_rest] # Combine nsigma and 'rest' levels for the final result. - result = np.concatenate([result_nsigma_levs, result_rest_levs], + result = da.concatenate([result_nsigma_levs, result_rest_levs], axis=i_levels_dim) return result @@ -1424,8 +1425,8 @@ def dependencies(self): b=self.b, depth_c=self.depth_c) def _derive(self, s, eta, depth, a, b, depth_c): - c = ((1 - b) * np.sinh(a * s) / np.sinh(a) + b * - (np.tanh(a * (s + 0.5)) / (2 * np.tanh(0.5 * a)) - 0.5)) + c = ((1 - b) * da.sinh(a * s) / da.sinh(a) + b * + (da.tanh(a * (s + 0.5)) / (2 * da.tanh(0.5 * a)) - 0.5)) return eta * (1 + s) + depth_c * s + (depth - depth_c) * c def make_coord(self, coord_dims_func): diff --git a/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py b/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py index 9f3c08c784..9942a20e25 100644 --- a/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py +++ b/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py @@ -96,14 +96,24 @@ def _lazy_testcube(self): coord.points = as_lazy_data(coord.points, coord.shape) return cube - def test_lazy(self): + def test_nonlazy_cube_has_lazy_derived(self): # Check same results when key coords are made lazy. cube = self.cube - self.assertFalse( - cube.coord(self.derived_coord_name).has_lazy_points()) + self.assertEqual( + cube.coord('depth').has_lazy_points(), + False) + self.assertEqual( + cube.coord(self.derived_coord_name).has_lazy_points(), + True) + + def test_lazy_cube_same_result(self): cube = self._lazy_testcube() - self.assertTrue( - cube.coord(self.derived_coord_name).has_lazy_points()) + self.assertEqual( + cube.coord('depth').has_lazy_points(), + True) + self.assertEqual( + cube.coord(self.derived_coord_name).has_lazy_points(), + True) self._check_result(cube) def test_transpose(self): @@ -120,7 +130,7 @@ def test_transpose(self): def test_lazy_transpose(self): # Check lazy calc works with all possible dimension orders. for dims_list in itertools.permutations(range(self.cube.ndim)): - cube = self._lazy_testcube() + cube = self._lazy_testcube().copy() cube.transpose(dims_list) expected = self.basic_derived_result.transpose(dims_list) msg = 'Unexpected result when cube transposed by {}' @@ -129,7 +139,7 @@ def test_lazy_transpose(self): err_msg=msg) def test_extra_dims(self): - # Insert extra cube dimensions + check it still works. + # Insert some extra cube dimensions + check it still works. cube = self.cube cube = iris.util.new_axis(cube) cube = iris.util.new_axis(cube) @@ -142,10 +152,12 @@ def test_extra_dims(self): def test_no_sigma(self): # Check it still works when 'sigma' is removed ... + # Set all sigma points to zero + snapshot the resulting derived points. trial_cube = self.cube.copy() trial_cube.coord('ocean_sigma_z_coordinate').points[:] = 0.0 expected = trial_cube.coord(self.derived_coord_name).points + # Remove sigma altogether + check the result is the same. cube = self.cube cube.remove_coord('ocean_sigma_z_coordinate') @@ -163,7 +175,7 @@ def test_no_eta(self): # Take first time, as no sigma --> result *has* no time dimension. expected = expected[0] - # Remove sigma altogether + check the result is the same. + # Remove eta altogether + check the result is the same. cube = self.cube cube.remove_coord('sea_surface_height') self._check_result(cube, expected) From 507a4df23b657d0f1cfec94656a5229eebdc6833 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 22 Jun 2017 10:46:59 +0100 Subject: [PATCH 06/11] Fix nasty misuses of list.append. --- lib/iris/aux_factory.py | 104 ++++++---------------------------------- 1 file changed, 14 insertions(+), 90 deletions(-) diff --git a/lib/iris/aux_factory.py b/lib/iris/aux_factory.py index e3742d4ef7..d55427348a 100644 --- a/lib/iris/aux_factory.py +++ b/lib/iris/aux_factory.py @@ -308,23 +308,6 @@ def _remap_with_bounds(self, dependency_dims, derived_dims): nd_values_by_key[key] = nd_values return nd_values_by_key -# def _shape(self, nd_values_by_key): -# nd_values = sorted(nd_values_by_key.values(), -# key=lambda value: value.ndim) -# shape = list(nd_values.pop().shape) -# for array in nd_values: -# for i, size in enumerate(array.shape): -# if size > 1: -# # NB. If there's an inconsistency it can only come -# # from a mismatch in the number of bounds (the Cube -# # ensures the other dimensions must match). -# # But we can't afford to raise an error now - it'd -# # break Cube.derived_coords. Instead, we let the -# # error happen when the derived coordinate's bounds -# # are accessed. -# shape[i] = size -# return shape - class HybridHeightFactory(AuxCoordFactory): """ @@ -437,9 +420,8 @@ def calc_bounds(): 'These are being disregarded.', UserWarning, stacklevel=2) orography_pts = nd_points_by_key['orography'] - orography_pts_shape = list(orography_pts.shape) - orography = orography_pts.reshape( - orography_pts_shape.append(1)) + bds_shape = list(orography_pts.shape) + [1] + orography = orography_pts.reshape(bds_shape) return self._derive(delta, sigma, orography) bounds = calc_bounds() @@ -627,10 +609,9 @@ def calc_bounds(): 'These are being disregarded.') surface_air_pressure_pts = nd_points_by_key[ 'surface_air_pressure'] - surface_air_pressure_pts_shape = list( - surface_air_pressure_pts.shape) + bds_shape = list(surface_air_pressure_pts.shape) + [1] surface_air_pressure = surface_air_pressure_pts.reshape( - surface_air_pressure_pts_shape.append(1)) + bds_shape) return self._derive(delta, sigma, surface_air_pressure) bounds = calc_bounds() @@ -776,63 +757,6 @@ def dependencies(self): return dict(sigma=self.sigma, eta=self.eta, depth=self.depth, depth_c=self.depth_c, nsigma=self.nsigma, zlev=self.zlev) - def _original_derive(self, sigma, eta, depth, depth_c, - zlev, shape, nsigma_slice): - # Perform the ocean sigma over z coordinate nsigma slice. - if eta.ndim: - eta = eta[nsigma_slice] - if sigma.ndim: - sigma = sigma[nsigma_slice] - if depth.ndim: - depth = depth[nsigma_slice] - # Note that, this performs a point-wise minimum. - temp = eta + sigma * (np.minimum(depth_c, depth) + eta) - # Calculate the final derived result. - result = np.ones(shape, dtype=temp.dtype) * zlev - result[nsigma_slice] = temp - - return result - - def _oldstyle_derive(self, sigma, eta, depth, depth_c, - zlev, nsigma, coord_dims_func): - derived_cubedims = self.derived_dims(coord_dims_func) - i_levels_cubedim, = coord_dims_func(self.dependencies['zlev']) - i_levels_dim = i_levels_cubedim - sum( - i_dim not in derived_cubedims - for i_dim in range(i_levels_cubedim)) - allshapes = np.array( - [el.shape - for el in (sigma, eta, depth, depth_c, zlev) - if el.ndim]) - shape = list(np.max(allshapes, axis=0)) - ndims = len(shape) - - # N.B. it is not a problem forming 'shape': following checks out ok. - # derived_dims = self.derived_dims(coord_dims_func) - # dependency_dims = self._dependency_dims(coord_dims_func) - # nd_points_by_key = self._remap(dependency_dims, derived_dims) - # orig_shape = self._shape(nd_points_by_key) - # assert orig_shape == shape - - # Make a slice tuple to index the first nsigma z-levels. - nsigma_slice = [slice(None)] * ndims - nsigma_slice[i_levels_dim] = slice(0, int(nsigma)) - nsigma_slice = tuple(nsigma_slice) - # Perform the ocean sigma over z coordinate nsigma slice. - if eta.ndim: - eta = eta[nsigma_slice] - if sigma.ndim: - sigma = sigma[nsigma_slice] - if depth.ndim: - depth = depth[nsigma_slice] - # Note that, this performs a point-wise minimum. - temp = eta + sigma * (np.minimum(depth_c, depth) + eta) - # Calculate the final derived result. - result = np.ones(shape, dtype=temp.dtype) * zlev - result[nsigma_slice] = temp - - return result - def _derive(self, sigma, eta, depth, depth_c, zlev, nsigma, coord_dims_func): # Calculate the index of the 'z' dimension in the inputs. @@ -931,8 +855,8 @@ def calc_bounds(): 'These are being disregarded.'.format(key, name) warnings.warn(msg, UserWarning, stacklevel=2) # Swap bounds with points. - shape = list(nd_points_by_key[key].shape) - bounds = nd_points_by_key[key].reshape(shape.append(1)) + bds_shape = list(nd_points_by_key[key].shape) + [1] + bounds = nd_points_by_key[key].reshape(bds_shape) nd_values_by_key[key] = bounds return self._derive(nd_values_by_key['sigma'], nd_values_by_key['eta'], @@ -1100,8 +1024,8 @@ def calc_bounds(): 'These are being disregarded.'.format(key, name) warnings.warn(msg, UserWarning, stacklevel=2) # Swap bounds with points. - shape = list(nd_points_by_key[key].shape) - bounds = nd_points_by_key[key].reshape(shape.append(1)) + bds_shape = list(nd_points_by_key[key].shape) + [1] + bounds = nd_points_by_key[key].reshape(bds_shape) nd_values_by_key[key] = bounds return self._derive(nd_values_by_key['sigma'], nd_values_by_key['eta'], @@ -1288,8 +1212,8 @@ def calc_bounds(): 'These are being disregarded.'.format(key, name) warnings.warn(msg, UserWarning, stacklevel=2) # Swap bounds with points. - shape = list(nd_points_by_key[key].shape) - bounds = nd_points_by_key[key].reshape(shape.append(1)) + bds_shape = list(nd_points_by_key[key].shape) + [1] + bounds = nd_points_by_key[key].reshape(bds_shape) nd_values_by_key[key] = bounds return self._derive(nd_values_by_key['s'], nd_values_by_key['c'], @@ -1476,8 +1400,8 @@ def calc_bounds(): 'These are being disregarded.'.format(key, name) warnings.warn(msg, UserWarning, stacklevel=2) # Swap bounds with points. - shape = list(nd_points_by_key[key].shape) - bounds = nd_points_by_key[key].reshape(shape.append(1)) + bds_shape = list(nd_points_by_key[key].shape) + [1] + bounds = nd_points_by_key[key].reshape(bds_shape) nd_values_by_key[key] = bounds return self._derive(nd_values_by_key['s'], nd_values_by_key['eta'], @@ -1667,8 +1591,8 @@ def calc_bounds(): 'These are being disregarded.'.format(key, name) warnings.warn(msg, UserWarning, stacklevel=2) # Swap bounds with points. - shape = list(nd_points_by_key[key].shape) - bounds = nd_points_by_key[key].reshape(shape.append(1)) + bds_shape = list(nd_points_by_key[key].shape) + [1] + bounds = nd_points_by_key[key].reshape(bds_shape) nd_values_by_key[key] = bounds return self._derive(nd_values_by_key['s'], nd_values_by_key['c'], From f962c1bc78bffef121a7fe1200c591c647b0d93c Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 22 Jun 2017 11:48:21 +0100 Subject: [PATCH 07/11] Adjust testing for fixes to aux_factory code. --- lib/iris/tests/test_hybrid.py | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/lib/iris/tests/test_hybrid.py b/lib/iris/tests/test_hybrid.py index a3240d8578..6c83df9502 100644 --- a/lib/iris/tests/test_hybrid.py +++ b/lib/iris/tests/test_hybrid.py @@ -151,14 +151,15 @@ def test_bounded_orography(self): # Make sure altitude still works OK if orography was messed # with *after* altitude was created. - altitude = self.cube.coord('altitude') orog.bounds = np.zeros(orog.shape + (4,)) - self.assertIsInstance(altitude.bounds, np.ndarray) - # Make sure altitude.bounds now raises an error. - exp_emsg = 'operands could not be broadcast together' - with self.assertRaisesRegexp(ValueError, exp_emsg): - self.cube.coord('altitude') + # Check that altitude derivation now issues a warning. + msg = 'Orography.* bounds.* being disregarded' + with warnings.catch_warnings(): + # Cause all warnings to raise Exceptions + warnings.simplefilter("error") + with self.assertRaisesRegexp(UserWarning, msg): + self.cube.coord('altitude') @tests.skip_data @@ -230,15 +231,15 @@ def test_bounded_surface_pressure(self): # Make sure pressure still works OK if surface pressure was messed # with *after* pressure was created. - pressure = self.cube.coord('air_pressure') surface_pressure.bounds = np.zeros(surface_pressure.shape + (4,)) - self.assertIsInstance(pressure.bounds, np.ndarray) - - # Make sure pressure.bounds now raises an error. - exp_emsg = 'operands could not be broadcast together' - with self.assertRaisesRegexp(ValueError, exp_emsg): - self.cube.coord('air_pressure') + # Check that air_pressure derivation now issues a warning. + msg = 'Surface pressure.* bounds.* being disregarded' + with warnings.catch_warnings(): + # Cause all warnings to raise Exceptions + warnings.simplefilter("error") + with self.assertRaisesRegexp(UserWarning, msg): + self.cube.coord('air_pressure') if __name__ == "__main__": tests.main() From 884faa0343a50ffb33f643139eec0b571eb43e9a Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 22 Jun 2017 12:14:23 +0100 Subject: [PATCH 08/11] CML changes: missing altitude points NaN --> mask. --- .../analysis/regrid/linear_circular_grid.cml | 8 +- .../regrid/linear_masked_altitude.cml | 62 +++++++-------- .../analysis/regrid/linear_non_circular.cml | 8 +- .../regrid/linear_partial_overlap.cml | 48 ++++++------ .../analysis/regrid/nearest_circular_grid.cml | 8 +- .../regrid/nearest_masked_altitude.cml | 24 +++--- .../analysis/regrid/nearest_non_circular.cml | 8 +- .../regrid/nearest_partial_overlap.cml | 48 ++++++------ .../results/analysis/regrid/no_overlap.cml | 78 +++++++++---------- 9 files changed, 146 insertions(+), 146 deletions(-) diff --git a/lib/iris/tests/results/analysis/regrid/linear_circular_grid.cml b/lib/iris/tests/results/analysis/regrid/linear_circular_grid.cml index 09e68c9f4c..08c5a81e1f 100644 --- a/lib/iris/tests/results/analysis/regrid/linear_circular_grid.cml +++ b/lib/iris/tests/results/analysis/regrid/linear_circular_grid.cml @@ -7,10 +7,10 @@ - + diff --git a/lib/iris/tests/results/analysis/regrid/linear_masked_altitude.cml b/lib/iris/tests/results/analysis/regrid/linear_masked_altitude.cml index f4e9914aa6..e9b0aab378 100644 --- a/lib/iris/tests/results/analysis/regrid/linear_masked_altitude.cml +++ b/lib/iris/tests/results/analysis/regrid/linear_masked_altitude.cml @@ -8,60 +8,60 @@ + [[424.423, 398.043, --, --, --], + [368.688, 343.878, --, --, --], + [375.091, 347.861, --, --, --], + [446.161, 414.22, --, --, --]]]" shape="(2, 4, 5)" standard_name="altitude" units="Unit('m')" value_type="float32"> diff --git a/lib/iris/tests/results/analysis/regrid/linear_non_circular.cml b/lib/iris/tests/results/analysis/regrid/linear_non_circular.cml index 32a80bcce5..ded7ffd106 100644 --- a/lib/iris/tests/results/analysis/regrid/linear_non_circular.cml +++ b/lib/iris/tests/results/analysis/regrid/linear_non_circular.cml @@ -7,10 +7,10 @@ - + diff --git a/lib/iris/tests/results/analysis/regrid/linear_partial_overlap.cml b/lib/iris/tests/results/analysis/regrid/linear_partial_overlap.cml index e55e259bc8..147f5c69a1 100644 --- a/lib/iris/tests/results/analysis/regrid/linear_partial_overlap.cml +++ b/lib/iris/tests/results/analysis/regrid/linear_partial_overlap.cml @@ -6,54 +6,54 @@ - + [[--, --, 367.726, 355.63], + [--, --, 340.443, 407.574], + [--, --, 336.602, 419.093], + [--, --, 376.39, 341.021]]]" shape="(2, 4, 4)" standard_name="altitude" units="Unit('m')" value_type="float32"> diff --git a/lib/iris/tests/results/analysis/regrid/nearest_circular_grid.cml b/lib/iris/tests/results/analysis/regrid/nearest_circular_grid.cml index 7f040a6ff7..d73639f66d 100644 --- a/lib/iris/tests/results/analysis/regrid/nearest_circular_grid.cml +++ b/lib/iris/tests/results/analysis/regrid/nearest_circular_grid.cml @@ -7,10 +7,10 @@ - + diff --git a/lib/iris/tests/results/analysis/regrid/nearest_masked_altitude.cml b/lib/iris/tests/results/analysis/regrid/nearest_masked_altitude.cml index 83ceaea0a9..cfe47b3f8b 100644 --- a/lib/iris/tests/results/analysis/regrid/nearest_masked_altitude.cml +++ b/lib/iris/tests/results/analysis/regrid/nearest_masked_altitude.cml @@ -14,14 +14,14 @@ [[325.116, 337.95], [289.584, 302.472], - [nan, nan], - [nan, nan], + [--, --], + [--, --], [424.168, 436.85]], [[325.116, 337.95], [289.584, 302.472], - [nan, nan], - [nan, nan], + [--, --], + [--, --], [424.168, 436.85]], [[440.66, 453.316], @@ -39,14 +39,14 @@ [[337.95, 357.202], [302.472, 321.806], - [nan, nan], - [nan, nan], + [--, --], + [--, --], [436.85, 455.873]], [[337.95, 357.202], [302.472, 321.806], - [nan, nan], - [nan, nan], + [--, --], + [--, --], [436.85, 455.873]], [[453.316, 472.302], @@ -54,13 +54,13 @@ [315.637, 334.94], [315.637, 334.94], [272.726, 292.128]]]]" id="9041e969" points="[[[418.698, 379.592, 202.986, 202.986, 333.61], - [329.929, 294.417, nan, nan, 428.923], - [329.929, 294.417, nan, nan, 428.923], + [329.929, 294.417, --, --, 428.923], + [329.929, 294.417, --, --, 428.923], [445.406, 398.973, 307.595, 307.595, 264.642]], [[434.57, 395.539, 219.272, 219.272, 349.646], - [345.971, 310.528, nan, nan, 444.776], - [345.971, 310.528, nan, nan, 444.776], + [345.971, 310.528, --, --, 444.776], + [345.971, 310.528, --, --, 444.776], [461.227, 414.883, 323.68, 323.68, 280.81]]]" shape="(2, 4, 5)" standard_name="altitude" units="Unit('m')" value_type="float32"> diff --git a/lib/iris/tests/results/analysis/regrid/nearest_non_circular.cml b/lib/iris/tests/results/analysis/regrid/nearest_non_circular.cml index 1d8e87c7e4..5b842c8811 100644 --- a/lib/iris/tests/results/analysis/regrid/nearest_non_circular.cml +++ b/lib/iris/tests/results/analysis/regrid/nearest_non_circular.cml @@ -7,10 +7,10 @@ - + diff --git a/lib/iris/tests/results/analysis/regrid/nearest_partial_overlap.cml b/lib/iris/tests/results/analysis/regrid/nearest_partial_overlap.cml index 351e86a2ac..3763379f50 100644 --- a/lib/iris/tests/results/analysis/regrid/nearest_partial_overlap.cml +++ b/lib/iris/tests/results/analysis/regrid/nearest_partial_overlap.cml @@ -6,54 +6,54 @@ - + [[--, --, 395.539, 349.646], + [--, --, 310.528, 444.776], + [--, --, 310.528, 444.776], + [--, --, 414.883, 280.81]]]" shape="(2, 4, 4)" standard_name="altitude" units="Unit('m')" value_type="float32"> diff --git a/lib/iris/tests/results/analysis/regrid/no_overlap.cml b/lib/iris/tests/results/analysis/regrid/no_overlap.cml index 5168c88ca5..a4503e2532 100644 --- a/lib/iris/tests/results/analysis/regrid/no_overlap.cml +++ b/lib/iris/tests/results/analysis/regrid/no_overlap.cml @@ -6,54 +6,54 @@ - + [[--, --, --, --], + [--, --, --, --], + [--, --, --, --], + [--, --, --, --]]]" shape="(2, 4, 4)" standard_name="altitude" units="Unit('m')" value_type="float32"> From 48cf11f8ee4031bf772e28f904c4038714c329b1 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 22 Jun 2017 12:44:36 +0100 Subject: [PATCH 09/11] Clarify need for some integration testcases. --- .../integration/aux_factory/test_OceanSigmaZFactory.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py b/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py index 9942a20e25..a9245e0333 100644 --- a/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py +++ b/lib/iris/tests/integration/aux_factory/test_OceanSigmaZFactory.py @@ -151,7 +151,9 @@ def test_extra_dims(self): self._check_result(cube) def test_no_sigma(self): - # Check it still works when 'sigma' is removed ... + # Check it still works when 'sigma' is removed. + # NOTE: the unit test for this does not cover all cases because it + # doesn't provide a time dimension. # Set all sigma points to zero + snapshot the resulting derived points. trial_cube = self.cube.copy() @@ -164,7 +166,9 @@ def test_no_sigma(self): self._check_result(cube, expected) def test_no_eta(self): - # Check it still works when 'eta' is removed ... + # Check it still works when 'eta' is removed. + # NOTE: the unit test for this does not cover all cases because it + # doesn't provide a time dimension. # Set all sigma points to zero + snapshot the resulting derived points. trial_cube = self.cube.copy() From 8dab554bfb24aef375d87efbb0bf8aeb2111bc2f Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Fri, 23 Jun 2017 17:19:28 +0100 Subject: [PATCH 10/11] Review changes. --- lib/iris/aux_factory.py | 424 +++++++++++++++++++++------------------- 1 file changed, 228 insertions(+), 196 deletions(-) diff --git a/lib/iris/aux_factory.py b/lib/iris/aux_factory.py index d55427348a..6d5427ba65 100644 --- a/lib/iris/aux_factory.py +++ b/lib/iris/aux_factory.py @@ -111,7 +111,7 @@ def arg_text(item): def derived_dims(self, coord_dims_func): """ - Returns the virtual dim-mapping for the derived coordinate. + Returns the cube dimensions for the derived coordinate. Args: @@ -120,6 +120,10 @@ def derived_dims(self, coord_dims_func): to a given coordinate. See :meth:`iris.cube.Cube.coord_dims()`. + Returns: + + A sorted list of cube dimension numbers. + """ # Which dimensions are relevant? # e.g. If sigma -> [1] and orog -> [2, 3] then result = [1, 2, 3] @@ -174,8 +178,15 @@ def _dependency_dims(self, coord_dims_func): @staticmethod def _nd_bounds(coord, dims, ndim): """ - Returns the coord's bounds in Cube-orientation and - broadcastable to N dimensions. + Return a lazy bounds array for a dependency coordinate, 'coord'. + + The result is aligned to the first 'ndim' cube dimensions, and + expanded to the full ('ndim'+1)-dimensional shape. + + The value of 'ndim' must be >= the highest cube dimension of the + dependency coordinate. + + The extra final result dimension ('ndim'-th) is the bounds dimension. Example: coord.shape == (70,) @@ -190,7 +201,7 @@ def _nd_bounds(coord, dims, ndim): sorted_pairs = sorted(enumerate(dims), key=lambda pair: pair[1]) transpose_order = [pair[0] for pair in sorted_pairs] + [len(dims)] bounds = coord.lazy_bounds() - if dims: + if dims and transpose_order != list(range(len(dims))): bounds = bounds.transpose(transpose_order) # Figure out the n-dimensional shape. @@ -203,8 +214,13 @@ def _nd_bounds(coord, dims, ndim): @staticmethod def _nd_points(coord, dims, ndim): """ - Returns the coord's points in Cube-orientation and - broadcastable to N dimensions. + Return a lazy points array for a dependency coordinate, 'coord'. + + The result is aligned to the first 'ndim' cube dimensions, and + expanded to the full 'ndim'-dimensional shape. + + The value of 'ndim' must be >= the highest cube dimension of the + dependency coordinate. Example: coord.shape == (4, 3) @@ -236,6 +252,16 @@ def _nd_points(coord, dims, ndim): return points def _remap(self, dependency_dims, derived_dims): + """ + Return a mapping from dependency names to coordinate points arrays. + + For dependencies that are present, the values are all expanded and + aligned to the same dimensions, which is the full set of all the + dependency dimensions. + These non-missing values are all lazy arrays. + Missing dependencies, however, are assigned a scalar value of 0.0. + + """ if derived_dims: ndim = max(derived_dims) + 1 else: @@ -264,6 +290,20 @@ def _remap(self, dependency_dims, derived_dims): return nd_points_by_key def _remap_with_bounds(self, dependency_dims, derived_dims): + """ + Return a mapping from dependency names to coordinate bounds arrays. + + For dependencies that are present, the values are all expanded and + aligned to the same dimensions, which is the full set of all the + dependency dimensions, plus an extra bounds dimension. + These non-missing values are all lazy arrays. + Missing dependencies, however, are assigned a scalar value of 0.0. + + Where a dependency coordinate has no bounds, then the associated value + is taken from its points array, but reshaped to have an extra bounds + dimension of length 1. + + """ if derived_dims: ndim = max(derived_dims) + 1 else: @@ -404,26 +444,23 @@ def make_coord(self, coord_dims_func): # Build the bounds array. nd_values_by_key = self._remap_with_bounds(dependency_dims, derived_dims) - - # Define the function here to obtain a closure. - def calc_bounds(): - delta = nd_values_by_key['delta'] - sigma = nd_values_by_key['sigma'] - orography = nd_values_by_key['orography'] - ok_bound_shapes = [(), (1,), (2,)] - if delta.shape[-1:] not in ok_bound_shapes: - raise ValueError('Invalid delta coordinate bounds.') - if sigma.shape[-1:] not in ok_bound_shapes: - raise ValueError('Invalid sigma coordinate bounds.') - if orography.shape[-1:] not in [(), (1,)]: - warnings.warn('Orography coordinate has bounds. ' - 'These are being disregarded.', - UserWarning, stacklevel=2) - orography_pts = nd_points_by_key['orography'] - bds_shape = list(orography_pts.shape) + [1] - orography = orography_pts.reshape(bds_shape) - return self._derive(delta, sigma, orography) - bounds = calc_bounds() + delta = nd_values_by_key['delta'] + sigma = nd_values_by_key['sigma'] + orography = nd_values_by_key['orography'] + ok_bound_shapes = [(), (1,), (2,)] + if delta.shape[-1:] not in ok_bound_shapes: + raise ValueError('Invalid delta coordinate bounds.') + if sigma.shape[-1:] not in ok_bound_shapes: + raise ValueError('Invalid sigma coordinate bounds.') + if orography.shape[-1:] not in [(), (1,)]: + warnings.warn('Orography coordinate has bounds. ' + 'These are being disregarded.', + UserWarning, stacklevel=2) + orography_pts = nd_points_by_key['orography'] + bds_shape = list(orography_pts.shape) + [1] + orography = orography_pts.reshape(bds_shape) + + bounds = self._derive(delta, sigma, orography) hybrid_height = iris.coords.AuxCoord(points, standard_name=self.standard_name, @@ -593,27 +630,24 @@ def make_coord(self, coord_dims_func): # Build the bounds array. nd_values_by_key = self._remap_with_bounds(dependency_dims, derived_dims) - - # Define the function here to obtain a closure. - def calc_bounds(): - delta = nd_values_by_key['delta'] - sigma = nd_values_by_key['sigma'] - surface_air_pressure = nd_values_by_key['surface_air_pressure'] - ok_bound_shapes = [(), (1,), (2,)] - if delta.shape[-1:] not in ok_bound_shapes: - raise ValueError('Invalid delta coordinate bounds.') - if sigma.shape[-1:] not in ok_bound_shapes: - raise ValueError('Invalid sigma coordinate bounds.') - if surface_air_pressure.shape[-1:] not in [(), (1,)]: - warnings.warn('Surface pressure coordinate has bounds. ' - 'These are being disregarded.') - surface_air_pressure_pts = nd_points_by_key[ - 'surface_air_pressure'] - bds_shape = list(surface_air_pressure_pts.shape) + [1] - surface_air_pressure = surface_air_pressure_pts.reshape( - bds_shape) - return self._derive(delta, sigma, surface_air_pressure) - bounds = calc_bounds() + delta = nd_values_by_key['delta'] + sigma = nd_values_by_key['sigma'] + surface_air_pressure = nd_values_by_key['surface_air_pressure'] + ok_bound_shapes = [(), (1,), (2,)] + if delta.shape[-1:] not in ok_bound_shapes: + raise ValueError('Invalid delta coordinate bounds.') + if sigma.shape[-1:] not in ok_bound_shapes: + raise ValueError('Invalid sigma coordinate bounds.') + if surface_air_pressure.shape[-1:] not in [(), (1,)]: + warnings.warn('Surface pressure coordinate has bounds. ' + 'These are being disregarded.') + surface_air_pressure_pts = nd_points_by_key[ + 'surface_air_pressure'] + bds_shape = list(surface_air_pressure_pts.shape) + [1] + surface_air_pressure = surface_air_pressure_pts.reshape( + bds_shape) + + bounds = self._derive(delta, sigma, surface_air_pressure) hybrid_pressure = iris.coords.AuxCoord( points, standard_name=self.standard_name, long_name=self.long_name, @@ -759,30 +793,32 @@ def dependencies(self): def _derive(self, sigma, eta, depth, depth_c, zlev, nsigma, coord_dims_func): - # Calculate the index of the 'z' dimension in the inputs. - # Get the cube dimension... - i_levels_cubedim, = coord_dims_func(self.dependencies['zlev']) - # ...subtract any previous dimensions not used in the calculation. + # Calculate the index of the 'z' dimension in the input arrays. + # First find the cube 'z' dimension ... + [cube_z_dim] = coord_dims_func(self.dependencies['zlev']) + # ... then calculate the corresponding dependency dimension. derived_cubedims = self.derived_dims(coord_dims_func) - i_levels_dim = i_levels_cubedim - sum( - i_dim not in derived_cubedims - for i_dim in range(i_levels_cubedim)) + z_dim = derived_cubedims.index(cube_z_dim) # Calculate the result shape as a combination of all the inputs. + # Note: all the inputs have the same number of dimensions >= 1, except + # for any missing dependencies, which have scalar values. allshapes = np.array( [el.shape for el in (sigma, eta, depth, depth_c, zlev) - if el.ndim]) + if el.ndim > 0]) result_shape = list(np.max(allshapes, axis=0)) ndims = len(result_shape) + # Make a slice tuple to index the first nsigma z-levels. z_slices_nsigma = [slice(None)] * ndims - z_slices_nsigma[i_levels_dim] = slice(0, int(nsigma)) + z_slices_nsigma[z_dim] = slice(0, int(nsigma)) z_slices_nsigma = tuple(z_slices_nsigma) # Make a slice tuple to index the remaining z-levels. z_slices_rest = [slice(None)] * ndims - z_slices_rest[i_levels_dim] = slice(int(nsigma), None) + z_slices_rest[z_dim] = slice(int(nsigma), None) z_slices_rest = tuple(z_slices_rest) + # Perform the ocean sigma over z coordinate nsigma slice. if eta.ndim: eta = eta[z_slices_nsigma] @@ -792,17 +828,29 @@ def _derive(self, sigma, eta, depth, depth_c, depth = depth[z_slices_nsigma] # Note that, this performs a point-wise minimum. nsigma_levs = eta + sigma * (da.minimum(depth_c, depth) + eta) - # Expand to full shape, as it may sometimes have lower dimensionality. - ones_full_result = np.ones(result_shape, dtype=np.int16) - ones_nsigma_result = ones_full_result[z_slices_nsigma] - result_nsigma_levs = nsigma_levs * ones_nsigma_result - # Likewise, expand zlev to a full result shape. - zlev = zlev * ones_full_result - # From this, take the 'remaining' levels for the result. - result_rest_levs = zlev[z_slices_rest] + + # Make a result-shaped lazy "ones" array for expanding partial results. + # Note: for the 'chunks' arg, we try to use [1, 1, ... ny, nx]. + # This calculation could be assuming too much in some cases, as we + # don't actually check the dimensions of our dependencies anywhere. + result_chunks = result_shape + if len(result_shape) > 1: + result_chunks = [1] * len(result_shape) + result_chunks[-2:] = result_shape[-2:] + ones_full_result = da.ones(result_shape, chunks=result_chunks, + dtype=zlev.dtype) + + # Expand nsigma_levs to its full required shape : needed as the + # calculated result may have a fixed size of 1 in some dimensions. + result_nsigma_levs = nsigma_levs * ones_full_result[z_slices_nsigma] + + # Likewise, expand zlev to its full required shape. + result_rest_levs = (zlev[z_slices_rest] * + ones_full_result[z_slices_rest]) + # Combine nsigma and 'rest' levels for the final result. result = da.concatenate([result_nsigma_levs, result_rest_levs], - axis=i_levels_dim) + axis=z_dim) return result def make_coord(self, coord_dims_func): @@ -823,7 +871,7 @@ def make_coord(self, coord_dims_func): # Build the points array. nd_points_by_key = self._remap(dependency_dims, derived_dims) - nsigma, = nd_points_by_key['nsigma'] + [nsigma] = nd_points_by_key['nsigma'] points = self._derive(nd_points_by_key['sigma'], nd_points_by_key['eta'], nd_points_by_key['depth'], @@ -837,35 +885,32 @@ def make_coord(self, coord_dims_func): # Build the bounds array. nd_values_by_key = self._remap_with_bounds(dependency_dims, derived_dims) - - # Define the function here to obtain a closure. - def calc_bounds(): - valid_shapes = [(), (1,), (2,)] - for key in ('sigma', 'zlev'): - if nd_values_by_key[key].shape[-1:] not in valid_shapes: - name = self.dependencies[key].name() - msg = 'Invalid bounds for {} ' \ - 'coordinate {!r}.'.format(key, name) - raise ValueError(msg) - valid_shapes.pop() - for key in ('eta', 'depth', 'depth_c', 'nsigma'): - if nd_values_by_key[key].shape[-1:] not in valid_shapes: - name = self.dependencies[key].name() - msg = 'The {} coordinate {!r} has bounds. ' \ - 'These are being disregarded.'.format(key, name) - warnings.warn(msg, UserWarning, stacklevel=2) - # Swap bounds with points. - bds_shape = list(nd_points_by_key[key].shape) + [1] - bounds = nd_points_by_key[key].reshape(bds_shape) - nd_values_by_key[key] = bounds - return self._derive(nd_values_by_key['sigma'], - nd_values_by_key['eta'], - nd_values_by_key['depth'], - nd_values_by_key['depth_c'], - nd_values_by_key['zlev'], - nsigma, - coord_dims_func) - bounds = calc_bounds() + valid_shapes = [(), (1,), (2,)] + for key in ('sigma', 'zlev'): + if nd_values_by_key[key].shape[-1:] not in valid_shapes: + name = self.dependencies[key].name() + msg = 'Invalid bounds for {} ' \ + 'coordinate {!r}.'.format(key, name) + raise ValueError(msg) + valid_shapes.pop() + for key in ('eta', 'depth', 'depth_c', 'nsigma'): + if nd_values_by_key[key].shape[-1:] not in valid_shapes: + name = self.dependencies[key].name() + msg = 'The {} coordinate {!r} has bounds. ' \ + 'These are being disregarded.'.format(key, name) + warnings.warn(msg, UserWarning, stacklevel=2) + # Swap bounds with points. + bds_shape = list(nd_points_by_key[key].shape) + [1] + bounds = nd_points_by_key[key].reshape(bds_shape) + nd_values_by_key[key] = bounds + + bounds = self._derive(nd_values_by_key['sigma'], + nd_values_by_key['eta'], + nd_values_by_key['depth'], + nd_values_by_key['depth_c'], + nd_values_by_key['zlev'], + nsigma, + coord_dims_func) coord = iris.coords.AuxCoord(points, standard_name=self.standard_name, @@ -1006,32 +1051,28 @@ def make_coord(self, coord_dims_func): # Build the bounds array. nd_values_by_key = self._remap_with_bounds(dependency_dims, derived_dims) - - # Define the function here to obtain a closure. - def calc_bounds(): - valid_shapes = [(), (1,), (2,)] - key = 'sigma' + valid_shapes = [(), (1,), (2,)] + key = 'sigma' + if nd_values_by_key[key].shape[-1:] not in valid_shapes: + name = self.dependencies[key].name() + msg = 'Invalid bounds for {} ' \ + 'coordinate {!r}.'.format(key, name) + raise ValueError(msg) + valid_shapes.pop() + for key in ('eta', 'depth'): if nd_values_by_key[key].shape[-1:] not in valid_shapes: name = self.dependencies[key].name() - msg = 'Invalid bounds for {} ' \ - 'coordinate {!r}.'.format(key, name) - raise ValueError(msg) - valid_shapes.pop() - for key in ('eta', 'depth'): - if nd_values_by_key[key].shape[-1:] not in valid_shapes: - name = self.dependencies[key].name() - msg = 'The {} coordinate {!r} has bounds. ' \ - 'These are being disregarded.'.format(key, name) - warnings.warn(msg, UserWarning, stacklevel=2) - # Swap bounds with points. - bds_shape = list(nd_points_by_key[key].shape) + [1] - bounds = nd_points_by_key[key].reshape(bds_shape) - nd_values_by_key[key] = bounds - return self._derive(nd_values_by_key['sigma'], - nd_values_by_key['eta'], - nd_values_by_key['depth']) - - bounds = calc_bounds() + msg = 'The {} coordinate {!r} has bounds. ' \ + 'These are being disregarded.'.format(key, name) + warnings.warn(msg, UserWarning, stacklevel=2) + # Swap bounds with points. + bds_shape = list(nd_points_by_key[key].shape) + [1] + bounds = nd_points_by_key[key].reshape(bds_shape) + nd_values_by_key[key] = bounds + + bounds = self._derive(nd_values_by_key['sigma'], + nd_values_by_key['eta'], + nd_values_by_key['depth']) coord = iris.coords.AuxCoord(points, standard_name=self.standard_name, @@ -1194,33 +1235,30 @@ def make_coord(self, coord_dims_func): # Build the bounds array. nd_values_by_key = self._remap_with_bounds(dependency_dims, derived_dims) - - # Define the function here to obtain a closure. - def calc_bounds(): - valid_shapes = [(), (1,), (2,)] - key = 's' + valid_shapes = [(), (1,), (2,)] + key = 's' + if nd_values_by_key[key].shape[-1:] not in valid_shapes: + name = self.dependencies[key].name() + msg = 'Invalid bounds for {} ' \ + 'coordinate {!r}.'.format(key, name) + raise ValueError(msg) + valid_shapes.pop() + for key in ('eta', 'depth', 'depth_c'): if nd_values_by_key[key].shape[-1:] not in valid_shapes: name = self.dependencies[key].name() - msg = 'Invalid bounds for {} ' \ - 'coordinate {!r}.'.format(key, name) - raise ValueError(msg) - valid_shapes.pop() - for key in ('eta', 'depth', 'depth_c'): - if nd_values_by_key[key].shape[-1:] not in valid_shapes: - name = self.dependencies[key].name() - msg = 'The {} coordinate {!r} has bounds. ' \ - 'These are being disregarded.'.format(key, name) - warnings.warn(msg, UserWarning, stacklevel=2) - # Swap bounds with points. - bds_shape = list(nd_points_by_key[key].shape) + [1] - bounds = nd_points_by_key[key].reshape(bds_shape) - nd_values_by_key[key] = bounds - return self._derive(nd_values_by_key['s'], - nd_values_by_key['c'], - nd_values_by_key['eta'], - nd_values_by_key['depth'], - nd_values_by_key['depth_c']) - bounds = calc_bounds() + msg = 'The {} coordinate {!r} has bounds. ' \ + 'These are being disregarded.'.format(key, name) + warnings.warn(msg, UserWarning, stacklevel=2) + # Swap bounds with points. + bds_shape = list(nd_points_by_key[key].shape) + [1] + bounds = nd_points_by_key[key].reshape(bds_shape) + nd_values_by_key[key] = bounds + + bounds = self._derive(nd_values_by_key['s'], + nd_values_by_key['c'], + nd_values_by_key['eta'], + nd_values_by_key['depth'], + nd_values_by_key['depth_c']) coord = iris.coords.AuxCoord(points, standard_name=self.standard_name, @@ -1382,34 +1420,31 @@ def make_coord(self, coord_dims_func): # Build the bounds array. nd_values_by_key = self._remap_with_bounds(dependency_dims, derived_dims) - - # Define the function here to obtain a closure. - def calc_bounds(): - valid_shapes = [(), (1,), (2,)] - key = 's' + valid_shapes = [(), (1,), (2,)] + key = 's' + if nd_values_by_key[key].shape[-1:] not in valid_shapes: + name = self.dependencies[key].name() + msg = 'Invalid bounds for {} ' \ + 'coordinate {!r}.'.format(key, name) + raise ValueError(msg) + valid_shapes.pop() + for key in ('eta', 'depth', 'a', 'b', 'depth_c'): if nd_values_by_key[key].shape[-1:] not in valid_shapes: name = self.dependencies[key].name() - msg = 'Invalid bounds for {} ' \ - 'coordinate {!r}.'.format(key, name) - raise ValueError(msg) - valid_shapes.pop() - for key in ('eta', 'depth', 'a', 'b', 'depth_c'): - if nd_values_by_key[key].shape[-1:] not in valid_shapes: - name = self.dependencies[key].name() - msg = 'The {} coordinate {!r} has bounds. ' \ - 'These are being disregarded.'.format(key, name) - warnings.warn(msg, UserWarning, stacklevel=2) - # Swap bounds with points. - bds_shape = list(nd_points_by_key[key].shape) + [1] - bounds = nd_points_by_key[key].reshape(bds_shape) - nd_values_by_key[key] = bounds - return self._derive(nd_values_by_key['s'], - nd_values_by_key['eta'], - nd_values_by_key['depth'], - nd_values_by_key['a'], - nd_values_by_key['b'], - nd_values_by_key['depth_c']) - bounds = calc_bounds() + msg = 'The {} coordinate {!r} has bounds. ' \ + 'These are being disregarded.'.format(key, name) + warnings.warn(msg, UserWarning, stacklevel=2) + # Swap bounds with points. + bds_shape = list(nd_points_by_key[key].shape) + [1] + bounds = nd_points_by_key[key].reshape(bds_shape) + nd_values_by_key[key] = bounds + + bounds = self._derive(nd_values_by_key['s'], + nd_values_by_key['eta'], + nd_values_by_key['depth'], + nd_values_by_key['a'], + nd_values_by_key['b'], + nd_values_by_key['depth_c']) coord = iris.coords.AuxCoord(points, standard_name=self.standard_name, @@ -1573,33 +1608,30 @@ def make_coord(self, coord_dims_func): # Build the bounds array. nd_values_by_key = self._remap_with_bounds(dependency_dims, derived_dims) - - # Define the function here to obtain a closure. - def calc_bounds(): - valid_shapes = [(), (1,), (2,)] - key = 's' + valid_shapes = [(), (1,), (2,)] + key = 's' + if nd_values_by_key[key].shape[-1:] not in valid_shapes: + name = self.dependencies[key].name() + msg = 'Invalid bounds for {} ' \ + 'coordinate {!r}.'.format(key, name) + raise ValueError(msg) + valid_shapes.pop() + for key in ('eta', 'depth', 'depth_c'): if nd_values_by_key[key].shape[-1:] not in valid_shapes: name = self.dependencies[key].name() - msg = 'Invalid bounds for {} ' \ - 'coordinate {!r}.'.format(key, name) - raise ValueError(msg) - valid_shapes.pop() - for key in ('eta', 'depth', 'depth_c'): - if nd_values_by_key[key].shape[-1:] not in valid_shapes: - name = self.dependencies[key].name() - msg = 'The {} coordinate {!r} has bounds. ' \ - 'These are being disregarded.'.format(key, name) - warnings.warn(msg, UserWarning, stacklevel=2) - # Swap bounds with points. - bds_shape = list(nd_points_by_key[key].shape) + [1] - bounds = nd_points_by_key[key].reshape(bds_shape) - nd_values_by_key[key] = bounds - return self._derive(nd_values_by_key['s'], - nd_values_by_key['c'], - nd_values_by_key['eta'], - nd_values_by_key['depth'], - nd_values_by_key['depth_c']) - bounds = calc_bounds() + msg = 'The {} coordinate {!r} has bounds. ' \ + 'These are being disregarded.'.format(key, name) + warnings.warn(msg, UserWarning, stacklevel=2) + # Swap bounds with points. + bds_shape = list(nd_points_by_key[key].shape) + [1] + bounds = nd_points_by_key[key].reshape(bds_shape) + nd_values_by_key[key] = bounds + + bounds = self._derive(nd_values_by_key['s'], + nd_values_by_key['c'], + nd_values_by_key['eta'], + nd_values_by_key['depth'], + nd_values_by_key['depth_c']) coord = iris.coords.AuxCoord(points, standard_name=self.standard_name, From 3a4a536322ccb6539dda8abc0ecc8e3d078515ef Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 27 Jun 2017 12:52:40 +0100 Subject: [PATCH 11/11] Reroute assertRaisesRegexp to prevent Python3 deprecation undermining warnings tests. --- lib/iris/tests/__init__.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/lib/iris/tests/__init__.py b/lib/iris/tests/__init__.py index d59b2a1678..b420bba961 100644 --- a/lib/iris/tests/__init__.py +++ b/lib/iris/tests/__init__.py @@ -512,6 +512,20 @@ def assertXMLElement(self, obj, reference_filename): def assertArrayEqual(self, a, b, err_msg=''): np.testing.assert_array_equal(a, b, err_msg=err_msg) + def assertRaisesRegexp(self, *args, **kwargs): + """ + Emulate the old :meth:`unittest.TestCase.assertRaisesRegexp`. + + Because the original function is now deprecated in Python 3. + Now calls :meth:`six.assertRaisesRegex()` (no final "p") instead. + It is the same, except for providing an additional 'msg' argument. + + """ + # Note: invoke via parent class to avoid recursion as, in Python 2, + # "six.assertRaisesRegex" calls getattr(self, 'assertRaisesRegexp'). + return six.assertRaisesRegex(super(IrisTest_nometa, self), + *args, **kwargs) + def _assertMaskedArray(self, assertion, a, b, strict, **kwargs): # Define helper function to extract unmasked values as a 1d # array.