diff --git a/lib/iris/analysis/_interpolate_backdoor.py b/lib/iris/analysis/_interpolate_backdoor.py index 9b7c0faf90..c7392fc0f6 100644 --- a/lib/iris/analysis/_interpolate_backdoor.py +++ b/lib/iris/analysis/_interpolate_backdoor.py @@ -63,35 +63,3 @@ def _warn_deprecated(msg=None): if msg is None: msg = _INTERPOLATE_DEPRECATION_WARNING iris_warn_deprecated(msg) - - -def nearest_neighbour_indices(cube, sample_points): - msg = (_INTERPOLATE_DEPRECATION_WARNING + '\n' + - 'Please replace usage of ' - 'iris.analysis.interpolate.nearest_neighbour_indices() ' - 'with iris.coords.Coord.nearest_neighbour_index()).') - _warn_deprecated(msg) - return _interp.nearest_neighbour_indices(cube, sample_points) - -nearest_neighbour_indices.__doc__ = _interp.nearest_neighbour_indices.__doc__ - - -def extract_nearest_neighbour(cube, sample_points): - msg = (_INTERPOLATE_DEPRECATION_WARNING + '\n' + - 'Please replace usage of ' - 'iris.analysis.interpolate.extract_nearest_neighbour() with ' - 'iris.cube.Cube.interpolate(..., scheme=iris.analysis.Nearest()).') - _warn_deprecated(msg) - return _interp.extract_nearest_neighbour(cube, sample_points) - -extract_nearest_neighbour.__doc__ = _interp.extract_nearest_neighbour.__doc__ - - -def regrid(source_cube, grid_cube, mode='bilinear', **kwargs): - msg = (_INTERPOLATE_DEPRECATION_WARNING + '\n' + - 'Please replace usage of iris.analysis.interpolate.regrid() ' - 'with iris.cube.Cube.regrid().') - _warn_deprecated(msg) - return _interp.regrid(source_cube, grid_cube, mode=mode, **kwargs) - -regrid.__doc__ = _interp.regrid.__doc__ diff --git a/lib/iris/analysis/_interpolate_private.py b/lib/iris/analysis/_interpolate_private.py index e2d68966ad..d928800767 100644 --- a/lib/iris/analysis/_interpolate_private.py +++ b/lib/iris/analysis/_interpolate_private.py @@ -44,348 +44,208 @@ import iris.exceptions -def nearest_neighbour_indices(cube, sample_points): +def linear(cube, sample_points, extrapolation_mode='linear'): """ - Returns the indices to select the data value(s) closest to the given coordinate point values. + Return a cube of the linearly interpolated points given the desired + sample points. - The sample_points mapping does not have to include coordinate values corresponding to all data - dimensions. Any dimensions unspecified will default to a full slice. + Given a list of tuple pairs mapping coordinates (or coordinate names) + to their desired values, return a cube with linearly interpolated values. + If more than one coordinate is specified, the linear interpolation will be + carried out in sequence, thus providing n-linear interpolation + (bi-linear, tri-linear, etc.). + + If the input cube's data is masked, the result cube will have a data + mask interpolated to the new sample points + + .. testsetup:: + + import numpy as np For example: - >>> cube = iris.load_cube(iris.sample_data_path('ostia_monthly.nc')) - >>> iris.analysis.interpolate.nearest_neighbour_indices(cube, [('latitude', 0), ('longitude', 10)]) - (slice(None, None, None), 9, 12) - >>> iris.analysis.interpolate.nearest_neighbour_indices(cube, [('latitude', 0)]) - (slice(None, None, None), 9, slice(None, None, None)) + >>> cube = iris.load_cube(iris.sample_data_path('air_temp.pp')) + >>> sample_points = [('latitude', np.linspace(-90, 90, 10)), + ... ('longitude', np.linspace(-180, 180, 20))] + >>> iris.analysis.interpolate.linear(cube, sample_points) + + + .. note:: + + By definition, linear interpolation requires all coordinates to + be 1-dimensional. + + .. note:: + + If a specified coordinate is single valued its value will be + extrapolated to the desired sample points by assuming a gradient of + zero. Args: - * cube: - An :class:`iris.cube.Cube`. + * cube + The cube to be interpolated. + * sample_points - A list of tuple pairs mapping coordinate instances or unique coordinate names in the cube to point values. + List of one or more tuple pairs mapping coordinate to desired + points to interpolate. Points may be a scalar or a numpy array + of values. Multi-dimensional coordinates are not supported. - Returns: - The tuple of indices which will select the point in the cube closest to the supplied coordinate values. + Kwargs: + + * extrapolation_mode - string - one of 'linear', 'nan' or 'error' + + * If 'linear' the point will be calculated by extending the + gradient of closest two points. + * If 'nan' the extrapolation point will be put as a NaN. + * If 'error' a value error will be raised notifying of the + attempted extrapolation. .. note:: - Nearest neighbour interpolation of multidimensional coordinates is not - yet supported. + If the source cube's data, or any of its resampled coordinates, + have an integer data type they will be promoted to a floating + point data type in the result. .. deprecated:: 1.10 The module :mod:`iris.analysis.interpolate` is deprecated. Please replace usage of - :func:`iris.analysis.interpolate.nearest_neighbour_indices` - with :meth:`iris.coords.Coord.nearest_neighbour_index`. + :func:`iris.analysis.interpolate.linear` + with :meth:`iris.cube.Cube.interpolate` using the scheme + :class:`iris.analysis.Linear`. """ if isinstance(sample_points, dict): - msg = ('Providing a dictionary to specify points is deprecated. ' - 'Please provide a list of (coordinate, values) pairs.') - warn_deprecated(msg) sample_points = list(sample_points.items()) - if sample_points: - try: - coord, values = sample_points[0] - except ValueError: - raise ValueError('Sample points must be a list of (coordinate, value) pairs. Got %r.' % sample_points) - - points = [] - for coord, values in sample_points: - if isinstance(coord, six.string_types): - coord = cube.coord(coord) - else: - coord = cube.coord(coord) - points.append((coord, values)) - sample_points = points - - # Build up a list of indices to span the cube. - indices = [slice(None, None)] * cube.ndim - - # Build up a dictionary which maps the cube's data dimensions to a list (which will later - # be populated by coordinates in the sample points list) - dim_to_coord_map = {} - for i in range(cube.ndim): - dim_to_coord_map[i] = [] - - # Iterate over all of the specifications provided by sample_points - for coord, point in sample_points: - data_dim = cube.coord_dims(coord) - - # If no data dimension then we don't need to make any modifications to indices. - if not data_dim: - continue - elif len(data_dim) > 1: - raise iris.exceptions.CoordinateMultiDimError("Nearest neighbour interpolation of multidimensional " - "coordinates is not supported.") - data_dim = data_dim[0] - - dim_to_coord_map[data_dim].append(coord) + # catch the case where a user passes a single (coord/name, value) pair rather than a list of pairs + if sample_points and not (isinstance(sample_points[0], collections.Container) and not isinstance(sample_points[0], six.string_types)): + raise TypeError('Expecting the sample points to be a list of tuple pairs representing (coord, points), got a list of %s.' % type(sample_points[0])) - #calculate the nearest neighbour - min_index = coord.nearest_neighbour_index(point) + scheme = Linear(extrapolation_mode) + return cube.interpolate(sample_points, scheme) - if getattr(coord, 'circular', False): - warnings.warn("Nearest neighbour on a circular coordinate may not be picking the nearest point.", DeprecationWarning) - # If the dimension has already been interpolated then assert that the index from this coordinate - # agrees with the index already calculated, otherwise we have a contradicting specification - if indices[data_dim] != slice(None, None) and min_index != indices[data_dim]: - raise ValueError('The coordinates provided (%s) over specify dimension %s.' % - (', '.join([coord.name() for coord in dim_to_coord_map[data_dim]]), data_dim)) - - indices[data_dim] = min_index +def _interp1d_rolls_y(): + """ + Determines if :class:`scipy.interpolate.interp1d` rolls its array `y` by + comparing the shape of y passed into interp1d to the shape of its internal + representation of y. - return tuple(indices) + SciPy v0.13.x+ no longer rolls the axis of its internal representation + of y so we test for this occurring to prevent us subsequently + extrapolating along the wrong axis. + For further information on this change see, for example: + * https://github.com/scipy/scipy/commit/0d906d0fc54388464603c63119b9e35c9a9c4601 + (the commit that introduced the change in behaviour). + * https://github.com/scipy/scipy/issues/2621 + (a discussion on the change - note the issue is not resolved + at time of writing). -def extract_nearest_neighbour(cube, sample_points): """ - Returns a new cube using data value(s) closest to the given coordinate point values. + y = np.arange(12).reshape(3, 4) + f = interp1d(np.arange(3), y, axis=0) + # If the initial shape of y and the shape internal to interp1d are *not* + # the same then scipy.interp1d rolls y. + return y.shape != f.y.shape - The sample_points mapping does not have to include coordinate values corresponding to all data - dimensions. Any dimensions unspecified will default to a full slice. - For example: - - >>> cube = iris.load_cube(iris.sample_data_path('ostia_monthly.nc')) - >>> iris.analysis.interpolate.extract_nearest_neighbour(cube, [('latitude', 0), ('longitude', 10)]) - - >>> iris.analysis.interpolate.extract_nearest_neighbour(cube, [('latitude', 0)]) - +class Linear1dExtrapolator(object): + """ + Extension class to :class:`scipy.interpolate.interp1d` to provide linear extrapolation. - Args: + See also: :mod:`scipy.interpolate`. - * cube: - An :class:`iris.cube.Cube`. - * sample_points - A list of tuple pairs mapping coordinate instances or unique coordinate names in the cube to point values. + .. deprecated :: 1.10 - Returns: - A cube that represents uninterpolated data as near to the given points as possible. + """ + roll_y = _interp1d_rolls_y() - .. deprecated:: 1.10 + def __init__(self, interpolator): + """ + Given an already created :class:`scipy.interpolate.interp1d` instance, return a callable object + which supports linear extrapolation. - The module :mod:`iris.analysis.interpolate` is deprecated. - Please replace usage of - :func:`iris.analysis.interpolate.extract_nearest_neighbour` - with :meth:`iris.cube.Cube.interpolate` using the scheme - :class:`iris.analysis.Nearest`. + .. deprecated :: 1.10 - """ - return cube[nearest_neighbour_indices(cube, sample_points)] + """ + self._interpolator = interpolator + self.x = interpolator.x + # Store the y values given to the interpolator. + self.y = interpolator.y + """ + The y values given to the interpolator object. + .. note:: These are stored with the interpolator.axis last. -def regrid(source_cube, grid_cube, mode='bilinear', **kwargs): - """ - Returns a new cube with values derived from the source_cube on the horizontal grid specified - by the grid_cube. + """ + # Roll interpolator.axis to the end if scipy no longer does it for us. + if not self.roll_y: + self.y = np.rollaxis(self.y, self._interpolator.axis, self.y.ndim) - Fundamental input requirements: - 1) Both cubes must have a CoordSystem. - 2) The source 'x' and 'y' coordinates must not share data dimensions with any other coordinates. + def all_points_in_range(self, requested_x): + """Given the x points, do all of the points sit inside the interpolation range.""" + test = (requested_x >= self.x[0]) & (requested_x <= self.x[-1]) + if isinstance(test, np.ndarray): + test = test.all() + return test - In addition, the algorithm currently used requires: - 3) Both CS instances must be compatible: - i.e. of the same type, with the same attribute values, and with compatible coordinates. - 4) No new data dimensions can be created. - 5) Source cube coordinates to map to a single dimension. + def __call__(self, requested_x): + if not self.all_points_in_range(requested_x): + # cast requested_x to a numpy array if it is not already. + if not isinstance(requested_x, np.ndarray): + requested_x = np.array(requested_x) - Args: + # we need to catch the special case of providing a single value... + remember_that_i_was_0d = requested_x.ndim == 0 - * source_cube: - An instance of :class:`iris.cube.Cube` which supplies the source data and metadata. - * grid_cube: - An instance of :class:`iris.cube.Cube` which supplies the horizontal grid definition. + requested_x = requested_x.flatten() - Kwargs: + gt = np.where(requested_x > self.x[-1])[0] + lt = np.where(requested_x < self.x[0])[0] + ok = np.where( (requested_x >= self.x[0]) & (requested_x <= self.x[-1]) )[0] - * mode (string): - Regridding interpolation algorithm to be applied, which may be one of the following: + data_shape = list(self.y.shape) + data_shape[-1] = len(requested_x) + result = np.empty(data_shape, dtype=self._interpolator(self.x[0]).dtype) - * 'bilinear' for bi-linear interpolation (default), see :func:`iris.analysis.interpolate.linear`. - * 'nearest' for nearest neighbour interpolation. + # Make a variable to represent the slice into the resultant data. (This will be updated in each of gt, lt & ok) + interpolator_result_index = [slice(None, None)] * self.y.ndim - Returns: - A new :class:`iris.cube.Cube` instance. + if len(ok) != 0: + interpolator_result_index[-1] = ok - .. note:: + r = self._interpolator(requested_x[ok]) + # Reshape the properly formed array to put the interpolator.axis last i.e. dims 0, 1, 2 -> 0, 2, 1 if axis = 1 + axes = list(range(r.ndim)) + del axes[self._interpolator.axis] + axes.append(self._interpolator.axis) - The masked status of values are currently ignored. See :func:\ -`~iris.experimental.regrid.regrid_bilinear_rectilinear_src_and_grid` - for regrid support with mask awareness. + result[interpolator_result_index] = r.transpose(axes) - .. deprecated:: 1.10 + if len(lt) != 0: + interpolator_result_index[-1] = lt - Please use :meth:`iris.cube.Cube.regrid` instead, with an appropriate - regridding scheme: + grad = (self.y[..., 1:2] - self.y[..., 0:1]) / (self.x[1] - self.x[0]) + result[interpolator_result_index] = self.y[..., 0:1] + (requested_x[lt] - self.x[0]) * grad - * For mode='bilinear', simply use the :class:`~iris.analysis.Linear` - scheme. + if len(gt) != 0: + interpolator_result_index[-1] = gt - * For mode='nearest', use the :class:`~iris.analysis.Nearest` scheme, - with extrapolation_mode='extrapolate', but be aware of the - following possible differences: + grad = (self.y[..., -1:] - self.y[..., -2:-1]) / (self.x[-1] - self.x[-2]) + result[interpolator_result_index] = self.y[..., -1:] + (requested_x[gt] - self.x[-1]) * grad - * Any missing result points, i.e. those which match source points - which are masked or NaN, are returned as as NaN values by this - routine. The 'Nearest' scheme, however, represents missing - results as masked points in a masked array. - *Which* points are missing is unchanged. + axes = list(range(len(interpolator_result_index))) + axes.insert(self._interpolator.axis, axes.pop(axes[-1])) + result = result.transpose(axes) - * Longitude wrapping for this routine is controlled by the - 'circular' property of the x coordinate. - The 'Nearest' scheme, however, *always* wraps any coords with - modular units, such as (correctly formed) longitudes. - Thus, behaviour can be different if "x_coord.circular" is - False : In that case, if the original non-longitude-wrapped - operation is required, it can be replicated by converting all - X and Y coordinates' units to '1' and removing their coordinate - systems. + if remember_that_i_was_0d: + new_shape = list(result.shape) + del new_shape[self._interpolator.axis] + result = result.reshape(new_shape) - """ - if mode == 'bilinear': - scheme = Linear(**kwargs) - return source_cube.regrid(grid_cube, scheme) - - # Condition 1 - source_cs = source_cube.coord_system(iris.coord_systems.CoordSystem) - grid_cs = grid_cube.coord_system(iris.coord_systems.CoordSystem) - if (source_cs is None) != (grid_cs is None): - raise ValueError("The source and grid cubes must both have a CoordSystem or both have None.") - - # Condition 2: We can only have one x coordinate and one y coordinate with the source CoordSystem, and those coordinates - # must be the only ones occupying their respective dimension - source_x = source_cube.coord(axis='x', coord_system=source_cs) - source_y = source_cube.coord(axis='y', coord_system=source_cs) - - source_x_dims = source_cube.coord_dims(source_x) - source_y_dims = source_cube.coord_dims(source_y) - - source_x_dim = None - if source_x_dims: - if len(source_x_dims) > 1: - raise ValueError('The source x coordinate may not describe more than one data dimension.') - source_x_dim = source_x_dims[0] - dim_sharers = ', '.join([coord.name() for coord in source_cube.coords(contains_dimension=source_x_dim) if coord is not source_x]) - if dim_sharers: - raise ValueError('No coordinates may share a dimension (dimension %s) with the x ' - 'coordinate, but (%s) do.' % (source_x_dim, dim_sharers)) - - source_y_dim = None - if source_y_dims: - if len(source_y_dims) > 1: - raise ValueError('The source y coordinate may not describe more than one data dimension.') - source_y_dim = source_y_dims[0] - dim_sharers = ', '.join([coord.name() for coord in source_cube.coords(contains_dimension=source_y_dim) if coord is not source_y]) - if dim_sharers: - raise ValueError('No coordinates may share a dimension (dimension %s) with the y ' - 'coordinate, but (%s) do.' % (source_y_dim, dim_sharers)) - - if source_x_dim is not None and source_y_dim == source_x_dim: - raise ValueError('The source x and y coords may not describe the same data dimension.') - - - # Condition 3 - # Check for compatible horizontal CSs. Currently that means they're exactly the same except for the coordinate - # values. - # The same kind of CS ... - compatible = (source_cs == grid_cs) - if compatible: - grid_x = grid_cube.coord(axis='x', coord_system=grid_cs) - grid_y = grid_cube.coord(axis='y', coord_system=grid_cs) - compatible = source_x.is_compatible(grid_x) and \ - source_y.is_compatible(grid_y) - if not compatible: - raise ValueError("The new grid must be defined on the same coordinate system, and have the same coordinate " - "metadata, as the source.") - - # Condition 4 - if grid_cube.coord_dims(grid_x) and not source_x_dims or \ - grid_cube.coord_dims(grid_y) and not source_y_dims: - raise ValueError("The new grid must not require additional data dimensions.") - - x_coord = grid_x.copy() - y_coord = grid_y.copy() - - - # - # Adjust the data array to match the new grid. - # - - # get the new shape of the data - new_shape = list(source_cube.shape) - if source_x_dims: - new_shape[source_x_dims[0]] = grid_x.shape[0] - if source_y_dims: - new_shape[source_y_dims[0]] = grid_y.shape[0] - - new_data = np.empty(new_shape, dtype=source_cube.data.dtype) - - # Prepare the index pattern which will be used to insert a single "column" of data. - # NB. A "column" is a slice constrained to a single XY point, which therefore extends over *all* the other axes. - # For an XYZ cube this means a column only extends over Z and corresponds to the normal definition of "column". - indices = [slice(None, None)] * new_data.ndim - - if mode == 'bilinear': - # Perform bilinear interpolation, passing through any keywords. - points_dict = [(source_x, list(x_coord.points)), (source_y, list(y_coord.points))] - new_data = linear(source_cube, points_dict, **kwargs).data - else: - # Perform nearest neighbour interpolation on each column in turn. - for iy, y in enumerate(y_coord.points): - for ix, x in enumerate(x_coord.points): - column_pos = [(source_x, x), (source_y, y)] - column_data = extract_nearest_neighbour(source_cube, column_pos).data - if source_y_dim is not None: - indices[source_y_dim] = iy - if source_x_dim is not None: - indices[source_x_dim] = ix - new_data[tuple(indices)] = column_data - - # Special case to make 0-dimensional results take the same form as NumPy - if new_data.shape == (): - new_data = new_data.flat[0] - - # Start with just the metadata and the re-sampled data... - new_cube = iris.cube.Cube(new_data) - new_cube.metadata = source_cube.metadata - new_cube.fill_value = source_cube.fill_value - - # ... and then copy across all the unaffected coordinates. - - # Record a mapping from old coordinate IDs to new coordinates, - # for subsequent use in creating updated aux_factories. - coord_mapping = {} - - def copy_coords(source_coords, add_method): - for coord in source_coords: - if coord is source_x or coord is source_y: - continue - dims = source_cube.coord_dims(coord) - new_coord = coord.copy() - add_method(new_coord, dims) - coord_mapping[id(coord)] = new_coord - - copy_coords(source_cube.dim_coords, new_cube.add_dim_coord) - copy_coords(source_cube.aux_coords, new_cube.add_aux_coord) - - for factory in source_cube.aux_factories: - new_cube.add_aux_factory(factory.updated(coord_mapping)) - - # Add the new coords - if source_x in source_cube.dim_coords: - new_cube.add_dim_coord(x_coord, source_x_dim) - else: - new_cube.add_aux_coord(x_coord, source_x_dims) - - if source_y in source_cube.dim_coords: - new_cube.add_dim_coord(y_coord, source_y_dim) - else: - new_cube.add_aux_coord(y_coord, source_y_dims) - - return new_cube + return result + else: + return self._interpolator(requested_x) diff --git a/lib/iris/tests/integration/test_regrid_equivalence.py b/lib/iris/tests/integration/test_regrid_equivalence.py index 9b65ea761c..a87c3684cd 100644 --- a/lib/iris/tests/integration/test_regrid_equivalence.py +++ b/lib/iris/tests/integration/test_regrid_equivalence.py @@ -30,7 +30,6 @@ import numpy as np import iris -from iris.analysis._interpolate_private import regrid from iris.analysis import Nearest from iris.cube import Cube from iris.coords import AuxCoord, DimCoord @@ -157,19 +156,8 @@ def test_wrapping_non_circular(self): [3., 4., 5.]] src_cube = grid_cube(src_x, src_y, data) dst_cube = grid_cube(dst_x, dst_y) - # Account for a behavioural difference in this case : - # The Nearest scheme does wrapping of modular coordinate values. - # Thus target of 352.0 --> -8.0, which is nearest to -10. - # This looks just like "circular" handling, but only because it happens - # to produce the same results *for nearest-neighbour in particular*. - if isinstance(self, TestInterpolateRegridNearest): - # interpolate.regrid --> Wrapping-free results (non-circular). - expected_result = [[3., 3., 4., 4., 5., 5., 5., 5.], - [3., 3., 4., 4., 5., 5., 5., 5.]] - else: - # cube regrid --> Wrapped results. - expected_result = [[4., 3., 4., 4., 5., 5., 3., 4.], - [4., 3., 4., 4., 5., 5., 3., 4.]] + expected_result = [[4., 3., 4., 4., 5., 5., 3., 4.], + [4., 3., 4., 4., 5., 5., 3., 4.]] _debug_data(src_cube, "noncircular SOURCE") result_cube = self.regrid(src_cube, dst_cube) _debug_data(result_cube, "noncircular RESULT") @@ -240,18 +228,6 @@ def test_source_nan(self): self.assertArrayEqual(result_cube.data, expected_result) -# perform identical tests on the old + new approaches -class TestInterpolateRegridNearest(MixinCheckingCode, tests.IrisTest): - def regrid(self, src_cube, dst_cube, - translate_nans_to_mask=False, **kwargs): - result = regrid(src_cube, dst_cube, mode='nearest') - data = result.data - if translate_nans_to_mask and np.any(np.isnan(data)): - data = np.ma.masked_array(data, mask=np.isnan(data)) - result.data = data - return result - - class TestCubeRegridNearest(MixinCheckingCode, tests.IrisTest): scheme = Nearest(extrapolation_mode='extrapolate')