From d1f8fe3c6deb3a811bebbd6f6dedcfe023ef5ff6 Mon Sep 17 00:00:00 2001 From: Alastair Gemmell Date: Tue, 27 Aug 2019 17:19:21 +0100 Subject: [PATCH 1/7] Move PointInCell class and associated classes and methods from experiemental to analysis. WIP - Still need to fix tests. --- lib/iris/analysis/__init__.py | 66 ++++- lib/iris/analysis/_regrid.py | 426 +++++++++++++++++++++++++++ lib/iris/experimental/regrid.py | 495 ++------------------------------ 3 files changed, 511 insertions(+), 476 deletions(-) diff --git a/lib/iris/analysis/__init__.py b/lib/iris/analysis/__init__.py index cc0ddc5788..3eede7f9e8 100644 --- a/lib/iris/analysis/__init__.py +++ b/lib/iris/analysis/__init__.py @@ -62,7 +62,7 @@ from iris.analysis._area_weighted import AreaWeightedRegridder from iris.analysis._interpolation import (EXTRAPOLATION_MODES, RectilinearInterpolator) -from iris.analysis._regrid import RectilinearRegridder +from iris.analysis._regrid import RectilinearRegridder, CurvilinearRegridder import iris.coords from iris.exceptions import LazyAggregatorError import iris._lazy_data @@ -2554,3 +2554,67 @@ def regridder(self, src_cube, target_grid): from iris.analysis.trajectory import \ UnstructuredNearestNeigbourRegridder return UnstructuredNearestNeigbourRegridder(src_cube, target_grid) + + +class PointInCell(object): + """ + This class describes the point-in-cell regridding scheme for use + typically with :meth:`iris.cube.Cube.regrid()`. + + The PointInCell regridder can regrid data from a source grid of any + dimensionality and in any coordinate system. + The location of each source point is specified by X and Y coordinates + mapped over the same cube dimensions, aka "grid dimensions" : the grid may + have any dimensionality. The X and Y coordinates must also have the same, + defined coord_system. + The weights, if specified, must have the same shape as the X and Y + coordinates. + The output grid can be any 'normal' XY grid, specified by *separate* X + and Y coordinates : That is, X and Y have two different cube dimensions. + The output X and Y coordinates must also have a common, specified + coord_system. + + """ + def __init__(self, weights=None): + """ + Point-in-cell regridding scheme suitable for regridding over one + or more orthogonal coordinates. + + Optional Args: + + * weights: + A :class:`numpy.ndarray` instance that defines the weights + for the grid cells of the source grid. Must have the same shape + as the data of the source grid. + If unspecified, equal weighting is assumed. + + """ + self.weights = weights + + def regridder(self, src_grid, target_grid): + """ + Creates a point-in-cell regridder to perform regridding from the + source grid to the target grid. + + Typically you should use :meth:`iris.cube.Cube.regrid` for + regridding a cube. There are, however, some situations when + constructing your own regridder is preferable. These are detailed in + the :ref:`user guide `. + + Args: + + * src_grid: + The :class:`~iris.cube.Cube` defining the source grid. + * target_grid: + The :class:`~iris.cube.Cube` defining the target grid. + + Returns: + A callable with the interface: + + `callable(cube)` + + where `cube` is a cube with the same grid as `src_grid` + that is to be regridded to the `target_grid`. + + """ + return CurvilinearRegridder(src_grid, target_grid, self.weights) diff --git a/lib/iris/analysis/_regrid.py b/lib/iris/analysis/_regrid.py index 83561a2b94..719fb64e26 100644 --- a/lib/iris/analysis/_regrid.py +++ b/lib/iris/analysis/_regrid.py @@ -34,6 +34,314 @@ from iris.util import _meshgrid +def _transform_xy_arrays(crs_from, x, y, crs_to): + """ + Transform 2d points between cartopy coordinate reference systems. + + NOTE: copied private function from iris.analysis.cartography. + + Args: + + * crs_from, crs_to (:class:`cartopy.crs.Projection`): + The coordinate reference systems. + * x, y (arrays): + point locations defined in 'crs_from'. + + Returns: + x, y : Arrays of locations defined in 'crs_to'. + + """ + pts = crs_to.transform_points(crs_from, x, y) + return pts[..., 0], pts[..., 1] + + +def _regrid_weighted_curvilinear_to_rectilinear__prepare( + src_cube, weights, grid_cube): + """ + First (setup) part of 'regrid_weighted_curvilinear_to_rectilinear'. + + Check inputs and calculate the sparse regrid matrix and related info. + The 'regrid info' returned can be re-used over many 2d slices. + + """ + if src_cube.aux_factories: + msg = 'All source cube derived coordinates will be ignored.' + warnings.warn(msg) + + # Get the source cube x and y 2D auxiliary coordinates. + sx, sy = src_cube.coord(axis='x'), src_cube.coord(axis='y') + # Get the target grid cube x and y dimension coordinates. + tx, ty = get_xy_dim_coords(grid_cube) + + if sx.units != sy.units: + msg = 'The source cube x ({!r}) and y ({!r}) coordinates must ' \ + 'have the same units.' + raise ValueError(msg.format(sx.name(), sy.name())) + + if src_cube.coord_dims(sx) != src_cube.coord_dims(sy): + msg = 'The source cube x ({!r}) and y ({!r}) coordinates must ' \ + 'map onto the same cube dimensions.' + raise ValueError(msg.format(sx.name(), sy.name())) + + if sx.coord_system != sy.coord_system: + msg = 'The source cube x ({!r}) and y ({!r}) coordinates must ' \ + 'have the same coordinate system.' + raise ValueError(msg.format(sx.name(), sy.name())) + + if sx.coord_system is None: + msg = ('The source X and Y coordinates must have a defined ' + 'coordinate system.') + raise ValueError(msg) + + if tx.units != ty.units: + msg = 'The target grid cube x ({!r}) and y ({!r}) coordinates must ' \ + 'have the same units.' + raise ValueError(msg.format(tx.name(), ty.name())) + + if tx.coord_system is None: + msg = ('The target X and Y coordinates must have a defined ' + 'coordinate system.') + raise ValueError(msg) + + if tx.coord_system != ty.coord_system: + msg = 'The target grid cube x ({!r}) and y ({!r}) coordinates must ' \ + 'have the same coordinate system.' + raise ValueError(msg.format(tx.name(), ty.name())) + + if weights is None: + weights = np.ones(sx.shape) + if weights.shape != sx.shape: + msg = ('Provided weights must have the same shape as the X and Y ' + 'coordinates.') + raise ValueError(msg) + + if not tx.has_bounds() or not tx.is_contiguous(): + msg = 'The target grid cube x ({!r})coordinate requires ' \ + 'contiguous bounds.' + raise ValueError(msg.format(tx.name())) + + if not ty.has_bounds() or not ty.is_contiguous(): + msg = 'The target grid cube y ({!r}) coordinate requires ' \ + 'contiguous bounds.' + raise ValueError(msg.format(ty.name())) + + def _src_align_and_flatten(coord): + # Return a flattened, unmasked copy of a coordinate's points array that + # will align with a flattened version of the source cube's data. + # + # PP-TODO: Should work with any cube dimensions for X and Y coords. + # Probably needs fixing anyway? + # + points = coord.points + if src_cube.coord_dims(coord) == (1, 0): + points = points.T + if points.shape != src_cube.shape: + msg = 'The shape of the points array of {!r} is not compatible ' \ + 'with the shape of {!r}.' + raise ValueError(msg.format(coord.name(), src_cube.name())) + return np.asarray(points.flatten()) + + # Align and flatten the coordinate points of the source space. + sx_points = _src_align_and_flatten(sx) + sy_points = _src_align_and_flatten(sy) + + # Transform source X and Y points into the target coord-system, if needed. + if sx.coord_system != tx.coord_system: + src_crs = sx.coord_system.as_cartopy_projection() + tgt_crs = tx.coord_system.as_cartopy_projection() + sx_points, sy_points = _transform_xy_arrays( + src_crs, sx_points, sy_points, tgt_crs) + # + # TODO: how does this work with scaled units ?? + # e.g. if crs is latlon, units could be degrees OR radians ? + # + + # Wrap modular values (e.g. longitudes) if required. + modulus = sx.units.modulus + if modulus is not None: + # Match the source cube x coordinate range to the target grid + # cube x coordinate range. + min_sx, min_tx = np.min(sx.points), np.min(tx.points) + if min_sx < 0 and min_tx >= 0: + indices = np.where(sx_points < 0) + # Ensure += doesn't raise a TypeError + if not np.can_cast(modulus, sx_points.dtype): + sx_points = sx_points.astype(type(modulus), casting='safe') + sx_points[indices] += modulus + elif min_sx >= 0 and min_tx < 0: + indices = np.where(sx_points > (modulus / 2)) + # Ensure -= doesn't raise a TypeError + if not np.can_cast(modulus, sx_points.dtype): + sx_points = sx_points.astype(type(modulus), casting='safe') + sx_points[indices] -= modulus + + # Create target grid cube x and y cell boundaries. + tx_depth, ty_depth = tx.points.size, ty.points.size + tx_dim, = grid_cube.coord_dims(tx) + ty_dim, = grid_cube.coord_dims(ty) + + tx_cells = np.concatenate((tx.bounds[:, 0], + tx.bounds[-1, 1].reshape(1))) + ty_cells = np.concatenate((ty.bounds[:, 0], + ty.bounds[-1, 1].reshape(1))) + + # Determine the target grid cube x and y cells that bound + # the source cube x and y points. + + def _regrid_indices(cells, depth, points): + # Calculate the minimum difference in cell extent. + extent = np.min(np.diff(cells)) + if extent == 0: + # Detected an dimension coordinate with an invalid + # zero length cell extent. + msg = 'The target grid cube {} ({!r}) coordinate contains ' \ + 'a zero length cell extent.' + axis, name = 'x', tx.name() + if points is sy_points: + axis, name = 'y', ty.name() + raise ValueError(msg.format(axis, name)) + elif extent > 0: + # The cells of the dimension coordinate are in ascending order. + indices = np.searchsorted(cells, points, side='right') - 1 + else: + # The cells of the dimension coordinate are in descending order. + # np.searchsorted() requires ascending order, so we require to + # account for this restriction. + cells = cells[::-1] + right = np.searchsorted(cells, points, side='right') + left = np.searchsorted(cells, points, side='left') + indices = depth - right + # Only those points that exactly match the left-hand cell bound + # will differ between 'left' and 'right'. Thus their appropriate + # target cell location requires to be recalculated to give the + # correct descending [upper, lower) interval cell, source to target + # regrid behaviour. + delta = np.where(left != right)[0] + if delta.size: + indices[delta] = depth - left[delta] + return indices + + x_indices = _regrid_indices(tx_cells, tx_depth, sx_points) + y_indices = _regrid_indices(ty_cells, ty_depth, sy_points) + + # Now construct a sparse M x N matix, where M is the flattened target + # space, and N is the flattened source space. The sparse matrix will then + # be populated with those source cube points that contribute to a specific + # target cube cell. + + # Determine the valid indices and their offsets in M x N space. + # Calculate the valid M offsets. + cols = np.where((y_indices >= 0) & (y_indices < ty_depth) & + (x_indices >= 0) & (x_indices < tx_depth))[0] + + # Reduce the indices to only those that are valid. + x_indices = x_indices[cols] + y_indices = y_indices[cols] + + # Calculate the valid N offsets. + if ty_dim < tx_dim: + rows = y_indices * tx.points.size + x_indices + else: + rows = x_indices * ty.points.size + y_indices + + # Calculate the associated valid weights. + weights_flat = weights.flatten() + data = weights_flat[cols] + + # Build our sparse M x N matrix of weights. + sparse_matrix = csc_matrix((data, (rows, cols)), + shape=(grid_cube.data.size, src_cube.data.size)) + + # Performing a sparse sum to collapse the matrix to (M, 1). + sum_weights = sparse_matrix.sum(axis=1).getA() + + # Determine the rows (flattened target indices) that have a + # contribution from one or more source points. + rows = np.nonzero(sum_weights) + + # NOTE: when source points are masked, this 'sum_weights' is possibly + # incorrect and needs re-calculating. Likewise 'rows' may cover target + # cells which happen to get no data. This is dealt with by adjusting as + # required in the '__perform' function, below. + + regrid_info = (sparse_matrix, sum_weights, rows, grid_cube) + return regrid_info + + +def _regrid_weighted_curvilinear_to_rectilinear__perform( + src_cube, regrid_info): + """ + Second (regrid) part of 'regrid_weighted_curvilinear_to_rectilinear'. + + Perform the prepared regrid calculation on a single 2d cube. + + """ + sparse_matrix, sum_weights, rows, grid_cube = regrid_info + + # Calculate the numerator of the weighted mean (M, 1). + is_masked = ma.isMaskedArray(src_cube.data) + if not is_masked: + data = src_cube.data + else: + # Use raw data array + data = src_cube.data.data + # Check if there are any masked source points to take account of. + is_masked = np.ma.is_masked(src_cube.data) + if is_masked: + # Zero any masked source points so they add nothing in output sums. + mask = src_cube.data.mask + data[mask] = 0.0 + # Calculate a new 'sum_weights' to allow for missing source points. + # N.B. it is more efficient to use the original once-calculated + # sparse matrix, but in this case we can't. + # Hopefully, this post-multiplying by the validities is less costly + # than repeating the whole sparse calculation. + valid_src_cells = ~mask.flat[:] + src_cell_validity_factors = sparse_diags( + np.array(valid_src_cells, dtype=int), + 0) + valid_weights = sparse_matrix * src_cell_validity_factors + sum_weights = valid_weights.sum(axis=1).getA() + # Work out where output cells are missing all contributions. + # This allows for where 'rows' contains output cells that have no + # data because of missing input points. + zero_sums = sum_weights == 0.0 + # Make sure we can still divide by sum_weights[rows]. + sum_weights[zero_sums] = 1.0 + + # Calculate sum in each target cell, over contributions from each source + # cell. + numerator = sparse_matrix * data.reshape(-1, 1) + + # Create a template for the weighted mean result. + weighted_mean = ma.masked_all(numerator.shape, dtype=numerator.dtype) + + # Calculate final results in all relevant places. + weighted_mean[rows] = numerator[rows] / sum_weights[rows] + if is_masked: + # Ensure masked points where relevant source cells were all missing. + if np.any(zero_sums): + # Make masked if it wasn't. + weighted_mean = np.ma.asarray(weighted_mean) + # Mask where contributing sums were zero. + weighted_mean[zero_sums] = np.ma.masked + + # Construct the final regridded weighted mean cube. + tx = grid_cube.coord(axis='x', dim_coords=True) + ty = grid_cube.coord(axis='y', dim_coords=True) + tx_dim, = grid_cube.coord_dims(tx) + ty_dim, = grid_cube.coord_dims(ty) + dim_coords_and_dims = list(zip((ty.copy(), tx.copy()), (ty_dim, tx_dim))) + cube = iris.cube.Cube(weighted_mean.reshape(grid_cube.shape), + dim_coords_and_dims=dim_coords_and_dims) + cube.metadata = copy.deepcopy(src_cube.metadata) + + for coord in src_cube.coords(dimensions=()): + cube.add_aux_coord(coord.copy()) + + return cube + + class RectilinearRegridder(object): """ This class provides support for performing nearest-neighbour or @@ -545,3 +853,121 @@ def __call__(self, src): sample_grid_x, sample_grid_y, regrid_callback) return result + + +class CurvilinearRegridder(object): + """ + This class provides support for performing point-in-cell regridding + between a curvilinear source grid and a rectilinear target grid. + + """ + def __init__(self, src_grid_cube, target_grid_cube, weights=None): + """ + Create a regridder for conversions between the source + and target grids. + + Args: + + * src_grid_cube: + The :class:`~iris.cube.Cube` providing the source grid. + * tgt_grid_cube: + The :class:`~iris.cube.Cube` providing the target grid. + + Optional Args: + + * weights: + A :class:`numpy.ndarray` instance that defines the weights + for the grid cells of the source grid. Must have the same shape + as the data of the source grid. + If unspecified, equal weighting is assumed. + + """ + # Validity checks. + if not isinstance(src_grid_cube, iris.cube.Cube): + raise TypeError("'src_grid_cube' must be a Cube") + if not isinstance(target_grid_cube, iris.cube.Cube): + raise TypeError("'target_grid_cube' must be a Cube") + # Snapshot the state of the cubes to ensure that the regridder + # is impervious to external changes to the original source cubes. + self._src_cube = src_grid_cube.copy() + self._target_cube = target_grid_cube.copy() + self.weights = weights + self._regrid_info = None + + @staticmethod + def _get_horizontal_coord(cube, axis): + """ + Gets the horizontal coordinate on the supplied cube along the + specified axis. + + Args: + + * cube: + An instance of :class:`iris.cube.Cube`. + * axis: + Locate coordinates on `cube` along this axis. + + Returns: + The horizontal coordinate on the specified axis of the supplied + cube. + + """ + coords = cube.coords(axis=axis, dim_coords=False) + if len(coords) != 1: + raise ValueError('Cube {!r} must contain a single 1D {} ' + 'coordinate.'.format(cube.name()), axis) + return coords[0] + + def __call__(self, src): + """ + Regrid the supplied :class:`~iris.cube.Cube` on to the target grid of + this :class:`_CurvilinearRegridder`. + + The given cube must be defined with the same grid as the source + grid used to create this :class:`_CurvilinearRegridder`. + + Args: + + * src: + A :class:`~iris.cube.Cube` to be regridded. + + Returns: + A cube defined with the horizontal dimensions of the target + and the other dimensions from this cube. The data values of + this cube will be converted to values on the new grid using + point-in-cell regridding. + + """ + # Validity checks. + if not isinstance(src, iris.cube.Cube): + raise TypeError("'src' must be a Cube") + + gx = self._get_horizontal_coord(self._src_cube, 'x') + gy = self._get_horizontal_coord(self._src_cube, 'y') + src_grid = (gx.copy(), gy.copy()) + sx = self._get_horizontal_coord(src, 'x') + sy = self._get_horizontal_coord(src, 'y') + if (sx, sy) != src_grid: + raise ValueError('The given cube is not defined on the same ' + 'source grid as this regridder.') + + # Call the regridder function. + # This includes repeating over any non-XY dimensions, because the + # underlying routine does not support this. + # FOR NOW: we will use cube.slices and merge to achieve this, + # though that is not a terribly efficient method ... + # TODO: create a template result cube and paste data slices into it, + # which would be more efficient. + result_slices = iris.cube.CubeList([]) + for slice_cube in src.slices(sx): + if self._regrid_info is None: + # Calculate the basic regrid info just once. + self._regrid_info = \ + _regrid_weighted_curvilinear_to_rectilinear__prepare( + slice_cube, self.weights, self._target_cube) + slice_result = \ + _regrid_weighted_curvilinear_to_rectilinear__perform( + slice_cube, self._regrid_info) + result_slices.append(slice_result) + result = result_slices.merge_cube() + return result diff --git a/lib/iris/experimental/regrid.py b/lib/iris/experimental/regrid.py index 91e48af61a..cd29476b12 100644 --- a/lib/iris/experimental/regrid.py +++ b/lib/iris/experimental/regrid.py @@ -32,16 +32,17 @@ import numpy as np import numpy.ma as ma import scipy.interpolate -from scipy.sparse import csc_matrix, diags as sparse_diags import six import iris.analysis.cartography from iris.analysis._interpolation import (get_xy_dim_coords, get_xy_coords, snapshot_grid) -from iris.analysis._regrid import RectilinearRegridder +from iris.analysis._regrid import (RectilinearRegridder, + _regrid_weighted_curvilinear_to_rectilinear__prepare, + _regrid_weighted_curvilinear_to_rectilinear__perform) import iris.coord_systems import iris.cube -from iris.util import _meshgrid, promote_aux_coord_to_dim_coord +from iris.util import _meshgrid _Version = namedtuple('Version', ('major', 'minor', 'micro')) @@ -780,27 +781,6 @@ def regrid_area_weighted_rectilinear_src_and_grid(src_cube, grid_cube, return new_cube -def _transform_xy_arrays(crs_from, x, y, crs_to): - """ - Transform 2d points between cartopy coordinate reference systems. - - NOTE: copied private function from iris.analysis.cartography. - - Args: - - * crs_from, crs_to (:class:`cartopy.crs.Projection`): - The coordinate reference systems. - * x, y (arrays): - point locations defined in 'crs_from'. - - Returns: - x, y : Arrays of locations defined in 'crs_to'. - - """ - pts = crs_to.transform_points(crs_from, x, y) - return pts[..., 0], pts[..., 1] - - def regrid_weighted_curvilinear_to_rectilinear(src_cube, weights, grid_cube): """ Return a new cube with the data values calculated using the weighted @@ -856,428 +836,17 @@ def regrid_weighted_curvilinear_to_rectilinear(src_cube, weights, grid_cube): return result -def _regrid_weighted_curvilinear_to_rectilinear__prepare( - src_cube, weights, grid_cube): - """ - First (setup) part of 'regrid_weighted_curvilinear_to_rectilinear'. - - Check inputs and calculate the sparse regrid matrix and related info. - The 'regrid info' returned can be re-used over many 2d slices. - - """ - if src_cube.aux_factories: - msg = 'All source cube derived coordinates will be ignored.' - warnings.warn(msg) - - # Get the source cube x and y 2D auxiliary coordinates. - sx, sy = src_cube.coord(axis='x'), src_cube.coord(axis='y') - # Get the target grid cube x and y dimension coordinates. - tx, ty = get_xy_dim_coords(grid_cube) - - if sx.units != sy.units: - msg = 'The source cube x ({!r}) and y ({!r}) coordinates must ' \ - 'have the same units.' - raise ValueError(msg.format(sx.name(), sy.name())) - - if src_cube.coord_dims(sx) != src_cube.coord_dims(sy): - msg = 'The source cube x ({!r}) and y ({!r}) coordinates must ' \ - 'map onto the same cube dimensions.' - raise ValueError(msg.format(sx.name(), sy.name())) - - if sx.coord_system != sy.coord_system: - msg = 'The source cube x ({!r}) and y ({!r}) coordinates must ' \ - 'have the same coordinate system.' - raise ValueError(msg.format(sx.name(), sy.name())) - - if sx.coord_system is None: - msg = ('The source X and Y coordinates must have a defined ' - 'coordinate system.') - raise ValueError(msg) - - if tx.units != ty.units: - msg = 'The target grid cube x ({!r}) and y ({!r}) coordinates must ' \ - 'have the same units.' - raise ValueError(msg.format(tx.name(), ty.name())) - - if tx.coord_system is None: - msg = ('The target X and Y coordinates must have a defined ' - 'coordinate system.') - raise ValueError(msg) - - if tx.coord_system != ty.coord_system: - msg = 'The target grid cube x ({!r}) and y ({!r}) coordinates must ' \ - 'have the same coordinate system.' - raise ValueError(msg.format(tx.name(), ty.name())) - - if weights is None: - weights = np.ones(sx.shape) - if weights.shape != sx.shape: - msg = ('Provided weights must have the same shape as the X and Y ' - 'coordinates.') - raise ValueError(msg) - - if not tx.has_bounds() or not tx.is_contiguous(): - msg = 'The target grid cube x ({!r})coordinate requires ' \ - 'contiguous bounds.' - raise ValueError(msg.format(tx.name())) - - if not ty.has_bounds() or not ty.is_contiguous(): - msg = 'The target grid cube y ({!r}) coordinate requires ' \ - 'contiguous bounds.' - raise ValueError(msg.format(ty.name())) - - def _src_align_and_flatten(coord): - # Return a flattened, unmasked copy of a coordinate's points array that - # will align with a flattened version of the source cube's data. - # - # PP-TODO: Should work with any cube dimensions for X and Y coords. - # Probably needs fixing anyway? - # - points = coord.points - if src_cube.coord_dims(coord) == (1, 0): - points = points.T - if points.shape != src_cube.shape: - msg = 'The shape of the points array of {!r} is not compatible ' \ - 'with the shape of {!r}.' - raise ValueError(msg.format(coord.name(), src_cube.name())) - return np.asarray(points.flatten()) - - # Align and flatten the coordinate points of the source space. - sx_points = _src_align_and_flatten(sx) - sy_points = _src_align_and_flatten(sy) - - # Transform source X and Y points into the target coord-system, if needed. - if sx.coord_system != tx.coord_system: - src_crs = sx.coord_system.as_cartopy_projection() - tgt_crs = tx.coord_system.as_cartopy_projection() - sx_points, sy_points = _transform_xy_arrays( - src_crs, sx_points, sy_points, tgt_crs) - # - # TODO: how does this work with scaled units ?? - # e.g. if crs is latlon, units could be degrees OR radians ? - # - - # Wrap modular values (e.g. longitudes) if required. - modulus = sx.units.modulus - if modulus is not None: - # Match the source cube x coordinate range to the target grid - # cube x coordinate range. - min_sx, min_tx = np.min(sx.points), np.min(tx.points) - if min_sx < 0 and min_tx >= 0: - indices = np.where(sx_points < 0) - # Ensure += doesn't raise a TypeError - if not np.can_cast(modulus, sx_points.dtype): - sx_points = sx_points.astype(type(modulus), casting='safe') - sx_points[indices] += modulus - elif min_sx >= 0 and min_tx < 0: - indices = np.where(sx_points > (modulus / 2)) - # Ensure -= doesn't raise a TypeError - if not np.can_cast(modulus, sx_points.dtype): - sx_points = sx_points.astype(type(modulus), casting='safe') - sx_points[indices] -= modulus - - # Create target grid cube x and y cell boundaries. - tx_depth, ty_depth = tx.points.size, ty.points.size - tx_dim, = grid_cube.coord_dims(tx) - ty_dim, = grid_cube.coord_dims(ty) - - tx_cells = np.concatenate((tx.bounds[:, 0], - tx.bounds[-1, 1].reshape(1))) - ty_cells = np.concatenate((ty.bounds[:, 0], - ty.bounds[-1, 1].reshape(1))) - - # Determine the target grid cube x and y cells that bound - # the source cube x and y points. - - def _regrid_indices(cells, depth, points): - # Calculate the minimum difference in cell extent. - extent = np.min(np.diff(cells)) - if extent == 0: - # Detected an dimension coordinate with an invalid - # zero length cell extent. - msg = 'The target grid cube {} ({!r}) coordinate contains ' \ - 'a zero length cell extent.' - axis, name = 'x', tx.name() - if points is sy_points: - axis, name = 'y', ty.name() - raise ValueError(msg.format(axis, name)) - elif extent > 0: - # The cells of the dimension coordinate are in ascending order. - indices = np.searchsorted(cells, points, side='right') - 1 - else: - # The cells of the dimension coordinate are in descending order. - # np.searchsorted() requires ascending order, so we require to - # account for this restriction. - cells = cells[::-1] - right = np.searchsorted(cells, points, side='right') - left = np.searchsorted(cells, points, side='left') - indices = depth - right - # Only those points that exactly match the left-hand cell bound - # will differ between 'left' and 'right'. Thus their appropriate - # target cell location requires to be recalculated to give the - # correct descending [upper, lower) interval cell, source to target - # regrid behaviour. - delta = np.where(left != right)[0] - if delta.size: - indices[delta] = depth - left[delta] - return indices - - x_indices = _regrid_indices(tx_cells, tx_depth, sx_points) - y_indices = _regrid_indices(ty_cells, ty_depth, sy_points) - - # Now construct a sparse M x N matix, where M is the flattened target - # space, and N is the flattened source space. The sparse matrix will then - # be populated with those source cube points that contribute to a specific - # target cube cell. - - # Determine the valid indices and their offsets in M x N space. - # Calculate the valid M offsets. - cols = np.where((y_indices >= 0) & (y_indices < ty_depth) & - (x_indices >= 0) & (x_indices < tx_depth))[0] - - # Reduce the indices to only those that are valid. - x_indices = x_indices[cols] - y_indices = y_indices[cols] - - # Calculate the valid N offsets. - if ty_dim < tx_dim: - rows = y_indices * tx.points.size + x_indices - else: - rows = x_indices * ty.points.size + y_indices - - # Calculate the associated valid weights. - weights_flat = weights.flatten() - data = weights_flat[cols] - - # Build our sparse M x N matrix of weights. - sparse_matrix = csc_matrix((data, (rows, cols)), - shape=(grid_cube.data.size, src_cube.data.size)) - - # Performing a sparse sum to collapse the matrix to (M, 1). - sum_weights = sparse_matrix.sum(axis=1).getA() - - # Determine the rows (flattened target indices) that have a - # contribution from one or more source points. - rows = np.nonzero(sum_weights) - - # NOTE: when source points are masked, this 'sum_weights' is possibly - # incorrect and needs re-calculating. Likewise 'rows' may cover target - # cells which happen to get no data. This is dealt with by adjusting as - # required in the '__perform' function, below. - - regrid_info = (sparse_matrix, sum_weights, rows, grid_cube) - return regrid_info - - -def _regrid_weighted_curvilinear_to_rectilinear__perform( - src_cube, regrid_info): - """ - Second (regrid) part of 'regrid_weighted_curvilinear_to_rectilinear'. - - Perform the prepared regrid calculation on a single 2d cube. - - """ - sparse_matrix, sum_weights, rows, grid_cube = regrid_info - - # Calculate the numerator of the weighted mean (M, 1). - is_masked = ma.isMaskedArray(src_cube.data) - if not is_masked: - data = src_cube.data - else: - # Use raw data array - data = src_cube.data.data - # Check if there are any masked source points to take account of. - is_masked = np.ma.is_masked(src_cube.data) - if is_masked: - # Zero any masked source points so they add nothing in output sums. - mask = src_cube.data.mask - data[mask] = 0.0 - # Calculate a new 'sum_weights' to allow for missing source points. - # N.B. it is more efficient to use the original once-calculated - # sparse matrix, but in this case we can't. - # Hopefully, this post-multiplying by the validities is less costly - # than repeating the whole sparse calculation. - valid_src_cells = ~mask.flat[:] - src_cell_validity_factors = sparse_diags( - np.array(valid_src_cells, dtype=int), - 0) - valid_weights = sparse_matrix * src_cell_validity_factors - sum_weights = valid_weights.sum(axis=1).getA() - # Work out where output cells are missing all contributions. - # This allows for where 'rows' contains output cells that have no - # data because of missing input points. - zero_sums = sum_weights == 0.0 - # Make sure we can still divide by sum_weights[rows]. - sum_weights[zero_sums] = 1.0 - - # Calculate sum in each target cell, over contributions from each source - # cell. - numerator = sparse_matrix * data.reshape(-1, 1) - - # Create a template for the weighted mean result. - weighted_mean = ma.masked_all(numerator.shape, dtype=numerator.dtype) - - # Calculate final results in all relevant places. - weighted_mean[rows] = numerator[rows] / sum_weights[rows] - if is_masked: - # Ensure masked points where relevant source cells were all missing. - if np.any(zero_sums): - # Make masked if it wasn't. - weighted_mean = np.ma.asarray(weighted_mean) - # Mask where contributing sums were zero. - weighted_mean[zero_sums] = np.ma.masked - - # Construct the final regridded weighted mean cube. - tx = grid_cube.coord(axis='x', dim_coords=True) - ty = grid_cube.coord(axis='y', dim_coords=True) - tx_dim, = grid_cube.coord_dims(tx) - ty_dim, = grid_cube.coord_dims(ty) - dim_coords_and_dims = list(zip((ty.copy(), tx.copy()), (ty_dim, tx_dim))) - cube = iris.cube.Cube(weighted_mean.reshape(grid_cube.shape), - dim_coords_and_dims=dim_coords_and_dims) - cube.metadata = copy.deepcopy(src_cube.metadata) - - for coord in src_cube.coords(dimensions=()): - cube.add_aux_coord(coord.copy()) - - return cube - - -class _CurvilinearRegridder(object): - """ - This class provides support for performing point-in-cell regridding - between a curvilinear source grid and a rectilinear target grid. - - """ - def __init__(self, src_grid_cube, target_grid_cube, weights=None): - """ - Create a regridder for conversions between the source - and target grids. - - Args: - - * src_grid_cube: - The :class:`~iris.cube.Cube` providing the source grid. - * tgt_grid_cube: - The :class:`~iris.cube.Cube` providing the target grid. - - Optional Args: - - * weights: - A :class:`numpy.ndarray` instance that defines the weights - for the grid cells of the source grid. Must have the same shape - as the data of the source grid. - If unspecified, equal weighting is assumed. - - """ - # Validity checks. - if not isinstance(src_grid_cube, iris.cube.Cube): - raise TypeError("'src_grid_cube' must be a Cube") - if not isinstance(target_grid_cube, iris.cube.Cube): - raise TypeError("'target_grid_cube' must be a Cube") - # Snapshot the state of the cubes to ensure that the regridder - # is impervious to external changes to the original source cubes. - self._src_cube = src_grid_cube.copy() - self._target_cube = target_grid_cube.copy() - self.weights = weights - self._regrid_info = None - - @staticmethod - def _get_horizontal_coord(cube, axis): - """ - Gets the horizontal coordinate on the supplied cube along the - specified axis. - - Args: - - * cube: - An instance of :class:`iris.cube.Cube`. - * axis: - Locate coordinates on `cube` along this axis. - - Returns: - The horizontal coordinate on the specified axis of the supplied - cube. - - """ - coords = cube.coords(axis=axis, dim_coords=False) - if len(coords) != 1: - raise ValueError('Cube {!r} must contain a single 1D {} ' - 'coordinate.'.format(cube.name()), axis) - return coords[0] - - def __call__(self, src): - """ - Regrid the supplied :class:`~iris.cube.Cube` on to the target grid of - this :class:`_CurvilinearRegridder`. - - The given cube must be defined with the same grid as the source - grid used to create this :class:`_CurvilinearRegridder`. - - Args: - - * src: - A :class:`~iris.cube.Cube` to be regridded. - - Returns: - A cube defined with the horizontal dimensions of the target - and the other dimensions from this cube. The data values of - this cube will be converted to values on the new grid using - point-in-cell regridding. - - """ - # Validity checks. - if not isinstance(src, iris.cube.Cube): - raise TypeError("'src' must be a Cube") - - gx = self._get_horizontal_coord(self._src_cube, 'x') - gy = self._get_horizontal_coord(self._src_cube, 'y') - src_grid = (gx.copy(), gy.copy()) - sx = self._get_horizontal_coord(src, 'x') - sy = self._get_horizontal_coord(src, 'y') - if (sx, sy) != src_grid: - raise ValueError('The given cube is not defined on the same ' - 'source grid as this regridder.') - - # Call the regridder function. - # This includes repeating over any non-XY dimensions, because the - # underlying routine does not support this. - # FOR NOW: we will use cube.slices and merge to achieve this, - # though that is not a terribly efficient method ... - # TODO: create a template result cube and paste data slices into it, - # which would be more efficient. - result_slices = iris.cube.CubeList([]) - for slice_cube in src.slices(sx): - if self._regrid_info is None: - # Calculate the basic regrid info just once. - self._regrid_info = \ - _regrid_weighted_curvilinear_to_rectilinear__prepare( - slice_cube, self.weights, self._target_cube) - slice_result = \ - _regrid_weighted_curvilinear_to_rectilinear__perform( - slice_cube, self._regrid_info) - result_slices.append(slice_result) - result = result_slices.merge_cube() - return result - - class PointInCell(object): """ This class describes the point-in-cell regridding scheme for use typically with :meth:`iris.cube.Cube.regrid()`. - The PointInCell regridder can regrid data from a source grid of any - dimensionality and in any coordinate system. - The location of each source point is specified by X and Y coordinates - mapped over the same cube dimensions, aka "grid dimensions" : the grid may - have any dimensionality. The X and Y coordinates must also have the same, - defined coord_system. - The weights, if specified, must have the same shape as the X and Y - coordinates. - The output grid can be any 'normal' XY grid, specified by *separate* X - and Y coordinates : That is, X and Y have two different cube dimensions. - The output X and Y coordinates must also have a common, specified - coord_system. + .. warning:: + + This class is now **disabled**. + + The functionality has been moved to + :class:`iris.analysis.PointInCell`. """ def __init__(self, weights=None): @@ -1285,44 +854,20 @@ def __init__(self, weights=None): Point-in-cell regridding scheme suitable for regridding over one or more orthogonal coordinates. - Optional Args: - - * weights: - A :class:`numpy.ndarray` instance that defines the weights - for the grid cells of the source grid. Must have the same shape - as the data of the source grid. - If unspecified, equal weighting is assumed. - - """ - self.weights = weights - - def regridder(self, src_grid, target_grid): - """ - Creates a point-in-cell regridder to perform regridding from the - source grid to the target grid. - - Typically you should use :meth:`iris.cube.Cube.regrid` for - regridding a cube. There are, however, some situations when - constructing your own regridder is preferable. These are detailed in - the :ref:`user guide `. - - Args: + .. warning:: - * src_grid: - The :class:`~iris.cube.Cube` defining the source grid. - * target_grid: - The :class:`~iris.cube.Cube` defining the target grid. + This class is now **disabled**. - Returns: - A callable with the interface: - - `callable(cube)` - - where `cube` is a cube with the same grid as `src_grid` - that is to be regridded to the `target_grid`. + The functionality has been moved to + :class:`iris.analysis.PointInCell`. """ - return _CurvilinearRegridder(src_grid, target_grid, self.weights) + raise Exception( + 'The class "iris.experimental.PointInCell" has been ' + 'moved, and is now in iris.analysis' + '\nPlease replace ' + '"iris.experimental.PointInCell" with ' + '"iris.analysis.PointInCell".') class _ProjectedUnstructuredRegridder(object): From 4dd104ee100240f98bfd8e7f7cdd14c72797ef38 Mon Sep 17 00:00:00 2001 From: Alastair Gemmell Date: Thu, 29 Aug 2019 13:59:30 +0100 Subject: [PATCH 2/7] Move PointInCell related tests from experiemental to analysis --- lib/iris/analysis/_regrid.py | 4 +++- lib/iris/experimental/regrid.py | 10 ++++++---- .../regrid/test_PointInCell.py | 6 +++--- .../regrid/test__CurvilinearRegridder.py | 8 ++++---- 4 files changed, 16 insertions(+), 12 deletions(-) rename lib/iris/tests/unit/{experimental => analysis}/regrid/test_PointInCell.py (89%) rename lib/iris/tests/unit/{experimental => analysis}/regrid/test__CurvilinearRegridder.py (98%) diff --git a/lib/iris/analysis/_regrid.py b/lib/iris/analysis/_regrid.py index 719fb64e26..897d9afe3d 100644 --- a/lib/iris/analysis/_regrid.py +++ b/lib/iris/analysis/_regrid.py @@ -1,4 +1,4 @@ -# (C) British Crown Copyright 2014 - 2018, Met Office +# (C) British Crown Copyright 2014 - 2019, Met Office # # This file is part of Iris. # @@ -33,6 +33,8 @@ import iris.cube from iris.util import _meshgrid +from scipy.sparse import csc_matrix, diags as sparse_diags + def _transform_xy_arrays(crs_from, x, y, crs_to): """ diff --git a/lib/iris/experimental/regrid.py b/lib/iris/experimental/regrid.py index cd29476b12..edcff3d3c4 100644 --- a/lib/iris/experimental/regrid.py +++ b/lib/iris/experimental/regrid.py @@ -1,4 +1,4 @@ -# (C) British Crown Copyright 2013 - 2018, Met Office +# (C) British Crown Copyright 2013 - 2019, Met Office # # This file is part of Iris. # @@ -37,9 +37,11 @@ import iris.analysis.cartography from iris.analysis._interpolation import (get_xy_dim_coords, get_xy_coords, snapshot_grid) -from iris.analysis._regrid import (RectilinearRegridder, - _regrid_weighted_curvilinear_to_rectilinear__prepare, - _regrid_weighted_curvilinear_to_rectilinear__perform) +from iris.analysis._regrid import ( + RectilinearRegridder, + _regrid_weighted_curvilinear_to_rectilinear__prepare, + _regrid_weighted_curvilinear_to_rectilinear__perform +) import iris.coord_systems import iris.cube from iris.util import _meshgrid diff --git a/lib/iris/tests/unit/experimental/regrid/test_PointInCell.py b/lib/iris/tests/unit/analysis/regrid/test_PointInCell.py similarity index 89% rename from lib/iris/tests/unit/experimental/regrid/test_PointInCell.py rename to lib/iris/tests/unit/analysis/regrid/test_PointInCell.py index 17be3845ec..4c09500fff 100644 --- a/lib/iris/tests/unit/experimental/regrid/test_PointInCell.py +++ b/lib/iris/tests/unit/analysis/regrid/test_PointInCell.py @@ -14,7 +14,7 @@ # # You should have received a copy of the GNU Lesser General Public License # along with Iris. If not, see . -"""Unit tests for :class:`iris.experimental.regrid.PointInCell`.""" +"""Unit tests for :class:`iris.analysis.PointInCell`.""" from __future__ import (absolute_import, division, print_function) from six.moves import (filter, input, map, range, zip) # noqa @@ -23,7 +23,7 @@ # importing anything else. import iris.tests as tests -from iris.experimental.regrid import PointInCell +from iris.analysis import PointInCell from iris.tests import mock @@ -31,7 +31,7 @@ class Test_regridder(tests.IrisTest): def test(self): point_in_cell = PointInCell(mock.sentinel.weights) - with mock.patch('iris.experimental.regrid._CurvilinearRegridder', + with mock.patch('iris.analysis._regrid.CurvilinearRegridder', return_value=mock.sentinel.regridder) as ecr: regridder = point_in_cell.regridder(mock.sentinel.src, mock.sentinel.target) diff --git a/lib/iris/tests/unit/experimental/regrid/test__CurvilinearRegridder.py b/lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py similarity index 98% rename from lib/iris/tests/unit/experimental/regrid/test__CurvilinearRegridder.py rename to lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py index 99dd0c483b..f89531b26f 100644 --- a/lib/iris/tests/unit/experimental/regrid/test__CurvilinearRegridder.py +++ b/lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py @@ -14,7 +14,7 @@ # # You should have received a copy of the GNU Lesser General Public License # along with Iris. If not, see . -"""Unit tests for :class:`iris.experimental.regrid._CurvilinearRegridder`.""" +"""Unit tests for :class:`iris.analysis._regrid._CurvilinearRegridder`.""" from __future__ import (absolute_import, division, print_function) from six.moves import (filter, input, map, range, zip) # noqa @@ -30,7 +30,7 @@ from iris.coords import AuxCoord, DimCoord from iris.coord_systems import GeogCS, RotatedGeogCS from iris.fileformats.pp import EARTH_RADIUS -from iris.experimental.regrid import _CurvilinearRegridder as Regridder +from iris.analysis._regrid import CurvilinearRegridder as Regridder from iris.tests import mock from iris.tests.stock import global_pp, lat_lon_cube @@ -57,10 +57,10 @@ def test_bad_grid_type(self): class Test___call__(tests.IrisTest): def setUp(self): self.func_setup = ( - 'iris.experimental.regrid.' + 'iris.analysis._regrid.' '_regrid_weighted_curvilinear_to_rectilinear__prepare') self.func_operate = ( - 'iris.experimental.regrid.' + 'iris.analysis.regrid.' '_regrid_weighted_curvilinear_to_rectilinear__perform') # Define a test source grid and target grid, basically the same. self.src_grid = global_pp() From b0c60d903eea03b0af667d5262a7e8f83e49b66e Mon Sep 17 00:00:00 2001 From: Alastair Gemmell Date: Wed, 4 Sep 2019 13:00:52 +0100 Subject: [PATCH 3/7] Fix broken PointInCell test, and move it from analysis/regrid to analysis --- lib/iris/tests/unit/analysis/{regrid => }/test_PointInCell.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) rename lib/iris/tests/unit/analysis/{regrid => }/test_PointInCell.py (96%) diff --git a/lib/iris/tests/unit/analysis/regrid/test_PointInCell.py b/lib/iris/tests/unit/analysis/test_PointInCell.py similarity index 96% rename from lib/iris/tests/unit/analysis/regrid/test_PointInCell.py rename to lib/iris/tests/unit/analysis/test_PointInCell.py index 4c09500fff..0a44267785 100644 --- a/lib/iris/tests/unit/analysis/regrid/test_PointInCell.py +++ b/lib/iris/tests/unit/analysis/test_PointInCell.py @@ -31,7 +31,7 @@ class Test_regridder(tests.IrisTest): def test(self): point_in_cell = PointInCell(mock.sentinel.weights) - with mock.patch('iris.analysis._regrid.CurvilinearRegridder', + with mock.patch('iris.analysis.CurvilinearRegridder', return_value=mock.sentinel.regridder) as ecr: regridder = point_in_cell.regridder(mock.sentinel.src, mock.sentinel.target) From d3941141a143727ee1b1fe8d33ee0a7e60aa5b87 Mon Sep 17 00:00:00 2001 From: Alastair Gemmell Date: Wed, 4 Sep 2019 13:21:53 +0100 Subject: [PATCH 4/7] Add whatsnew file --- ...04_move_experimental_point-in-cell_regridder_to_analysis.txt | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 docs/iris/src/whatsnew/contributions_2.3.0/newfeature_2019-Sep-04_move_experimental_point-in-cell_regridder_to_analysis.txt diff --git a/docs/iris/src/whatsnew/contributions_2.3.0/newfeature_2019-Sep-04_move_experimental_point-in-cell_regridder_to_analysis.txt b/docs/iris/src/whatsnew/contributions_2.3.0/newfeature_2019-Sep-04_move_experimental_point-in-cell_regridder_to_analysis.txt new file mode 100644 index 0000000000..65b3d8d568 --- /dev/null +++ b/docs/iris/src/whatsnew/contributions_2.3.0/newfeature_2019-Sep-04_move_experimental_point-in-cell_regridder_to_analysis.txt @@ -0,0 +1,2 @@ +* :class:`iris.experimental.regrid.PointInCell` moved to +:class:`iris.analysis.PointInCell` to make this regridding scheme public \ No newline at end of file From 3a104950f2ea8e79fab83a98cb9d4c0ff361aa77 Mon Sep 17 00:00:00 2001 From: Alastair Gemmell Date: Wed, 4 Sep 2019 17:13:10 +0100 Subject: [PATCH 5/7] Fix path in failing tests --- .../tests/unit/analysis/regrid/test__CurvilinearRegridder.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py b/lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py index f89531b26f..c02d05cb60 100644 --- a/lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py +++ b/lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py @@ -60,7 +60,7 @@ def setUp(self): 'iris.analysis._regrid.' '_regrid_weighted_curvilinear_to_rectilinear__prepare') self.func_operate = ( - 'iris.analysis.regrid.' + 'iris.analysis._regrid.' '_regrid_weighted_curvilinear_to_rectilinear__perform') # Define a test source grid and target grid, basically the same. self.src_grid = global_pp() From a7a513bd07906b3bc3292d1cf1428e6a6dcf4409 Mon Sep 17 00:00:00 2001 From: Alastair Gemmell Date: Wed, 4 Sep 2019 17:28:06 +0100 Subject: [PATCH 6/7] Update copyright file header --- .../tests/unit/analysis/regrid/test__CurvilinearRegridder.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py b/lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py index c02d05cb60..ad39a6b5c9 100644 --- a/lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py +++ b/lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py @@ -1,4 +1,4 @@ -# (C) British Crown Copyright 2015 - 2017, Met Office +# (C) British Crown Copyright 2015 - 2019, Met Office # # This file is part of Iris. # From 921cd7799290a0756d14af1cb720d8d0ab117b8a Mon Sep 17 00:00:00 2001 From: Alastair Gemmell Date: Mon, 9 Sep 2019 09:29:23 +0100 Subject: [PATCH 7/7] Minor PR tweaks in response to review comments --- lib/iris/analysis/_regrid.py | 236 +++++++++--------- .../regrid/test__CurvilinearRegridder.py | 2 +- .../tests/unit/analysis/test_PointInCell.py | 2 +- 3 files changed, 120 insertions(+), 120 deletions(-) diff --git a/lib/iris/analysis/_regrid.py b/lib/iris/analysis/_regrid.py index 897d9afe3d..f80d797642 100644 --- a/lib/iris/analysis/_regrid.py +++ b/lib/iris/analysis/_regrid.py @@ -344,6 +344,124 @@ def _regrid_weighted_curvilinear_to_rectilinear__perform( return cube +class CurvilinearRegridder(object): + """ + This class provides support for performing point-in-cell regridding + between a curvilinear source grid and a rectilinear target grid. + + """ + def __init__(self, src_grid_cube, target_grid_cube, weights=None): + """ + Create a regridder for conversions between the source + and target grids. + + Args: + + * src_grid_cube: + The :class:`~iris.cube.Cube` providing the source grid. + * tgt_grid_cube: + The :class:`~iris.cube.Cube` providing the target grid. + + Optional Args: + + * weights: + A :class:`numpy.ndarray` instance that defines the weights + for the grid cells of the source grid. Must have the same shape + as the data of the source grid. + If unspecified, equal weighting is assumed. + + """ + # Validity checks. + if not isinstance(src_grid_cube, iris.cube.Cube): + raise TypeError("'src_grid_cube' must be a Cube") + if not isinstance(target_grid_cube, iris.cube.Cube): + raise TypeError("'target_grid_cube' must be a Cube") + # Snapshot the state of the cubes to ensure that the regridder + # is impervious to external changes to the original source cubes. + self._src_cube = src_grid_cube.copy() + self._target_cube = target_grid_cube.copy() + self.weights = weights + self._regrid_info = None + + @staticmethod + def _get_horizontal_coord(cube, axis): + """ + Gets the horizontal coordinate on the supplied cube along the + specified axis. + + Args: + + * cube: + An instance of :class:`iris.cube.Cube`. + * axis: + Locate coordinates on `cube` along this axis. + + Returns: + The horizontal coordinate on the specified axis of the supplied + cube. + + """ + coords = cube.coords(axis=axis, dim_coords=False) + if len(coords) != 1: + raise ValueError('Cube {!r} must contain a single 1D {} ' + 'coordinate.'.format(cube.name()), axis) + return coords[0] + + def __call__(self, src): + """ + Regrid the supplied :class:`~iris.cube.Cube` on to the target grid of + this :class:`_CurvilinearRegridder`. + + The given cube must be defined with the same grid as the source + grid used to create this :class:`_CurvilinearRegridder`. + + Args: + + * src: + A :class:`~iris.cube.Cube` to be regridded. + + Returns: + A cube defined with the horizontal dimensions of the target + and the other dimensions from this cube. The data values of + this cube will be converted to values on the new grid using + point-in-cell regridding. + + """ + # Validity checks. + if not isinstance(src, iris.cube.Cube): + raise TypeError("'src' must be a Cube") + + gx = self._get_horizontal_coord(self._src_cube, 'x') + gy = self._get_horizontal_coord(self._src_cube, 'y') + src_grid = (gx.copy(), gy.copy()) + sx = self._get_horizontal_coord(src, 'x') + sy = self._get_horizontal_coord(src, 'y') + if (sx, sy) != src_grid: + raise ValueError('The given cube is not defined on the same ' + 'source grid as this regridder.') + + # Call the regridder function. + # This includes repeating over any non-XY dimensions, because the + # underlying routine does not support this. + # FOR NOW: we will use cube.slices and merge to achieve this, + # though that is not a terribly efficient method ... + # TODO: create a template result cube and paste data slices into it, + # which would be more efficient. + result_slices = iris.cube.CubeList([]) + for slice_cube in src.slices(sx): + if self._regrid_info is None: + # Calculate the basic regrid info just once. + self._regrid_info = \ + _regrid_weighted_curvilinear_to_rectilinear__prepare( + slice_cube, self.weights, self._target_cube) + slice_result = \ + _regrid_weighted_curvilinear_to_rectilinear__perform( + slice_cube, self._regrid_info) + result_slices.append(slice_result) + result = result_slices.merge_cube() + return result + + class RectilinearRegridder(object): """ This class provides support for performing nearest-neighbour or @@ -855,121 +973,3 @@ def __call__(self, src): sample_grid_x, sample_grid_y, regrid_callback) return result - - -class CurvilinearRegridder(object): - """ - This class provides support for performing point-in-cell regridding - between a curvilinear source grid and a rectilinear target grid. - - """ - def __init__(self, src_grid_cube, target_grid_cube, weights=None): - """ - Create a regridder for conversions between the source - and target grids. - - Args: - - * src_grid_cube: - The :class:`~iris.cube.Cube` providing the source grid. - * tgt_grid_cube: - The :class:`~iris.cube.Cube` providing the target grid. - - Optional Args: - - * weights: - A :class:`numpy.ndarray` instance that defines the weights - for the grid cells of the source grid. Must have the same shape - as the data of the source grid. - If unspecified, equal weighting is assumed. - - """ - # Validity checks. - if not isinstance(src_grid_cube, iris.cube.Cube): - raise TypeError("'src_grid_cube' must be a Cube") - if not isinstance(target_grid_cube, iris.cube.Cube): - raise TypeError("'target_grid_cube' must be a Cube") - # Snapshot the state of the cubes to ensure that the regridder - # is impervious to external changes to the original source cubes. - self._src_cube = src_grid_cube.copy() - self._target_cube = target_grid_cube.copy() - self.weights = weights - self._regrid_info = None - - @staticmethod - def _get_horizontal_coord(cube, axis): - """ - Gets the horizontal coordinate on the supplied cube along the - specified axis. - - Args: - - * cube: - An instance of :class:`iris.cube.Cube`. - * axis: - Locate coordinates on `cube` along this axis. - - Returns: - The horizontal coordinate on the specified axis of the supplied - cube. - - """ - coords = cube.coords(axis=axis, dim_coords=False) - if len(coords) != 1: - raise ValueError('Cube {!r} must contain a single 1D {} ' - 'coordinate.'.format(cube.name()), axis) - return coords[0] - - def __call__(self, src): - """ - Regrid the supplied :class:`~iris.cube.Cube` on to the target grid of - this :class:`_CurvilinearRegridder`. - - The given cube must be defined with the same grid as the source - grid used to create this :class:`_CurvilinearRegridder`. - - Args: - - * src: - A :class:`~iris.cube.Cube` to be regridded. - - Returns: - A cube defined with the horizontal dimensions of the target - and the other dimensions from this cube. The data values of - this cube will be converted to values on the new grid using - point-in-cell regridding. - - """ - # Validity checks. - if not isinstance(src, iris.cube.Cube): - raise TypeError("'src' must be a Cube") - - gx = self._get_horizontal_coord(self._src_cube, 'x') - gy = self._get_horizontal_coord(self._src_cube, 'y') - src_grid = (gx.copy(), gy.copy()) - sx = self._get_horizontal_coord(src, 'x') - sy = self._get_horizontal_coord(src, 'y') - if (sx, sy) != src_grid: - raise ValueError('The given cube is not defined on the same ' - 'source grid as this regridder.') - - # Call the regridder function. - # This includes repeating over any non-XY dimensions, because the - # underlying routine does not support this. - # FOR NOW: we will use cube.slices and merge to achieve this, - # though that is not a terribly efficient method ... - # TODO: create a template result cube and paste data slices into it, - # which would be more efficient. - result_slices = iris.cube.CubeList([]) - for slice_cube in src.slices(sx): - if self._regrid_info is None: - # Calculate the basic regrid info just once. - self._regrid_info = \ - _regrid_weighted_curvilinear_to_rectilinear__prepare( - slice_cube, self.weights, self._target_cube) - slice_result = \ - _regrid_weighted_curvilinear_to_rectilinear__perform( - slice_cube, self._regrid_info) - result_slices.append(slice_result) - result = result_slices.merge_cube() - return result diff --git a/lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py b/lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py index ad39a6b5c9..14997fc7a2 100644 --- a/lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py +++ b/lib/iris/tests/unit/analysis/regrid/test__CurvilinearRegridder.py @@ -14,7 +14,7 @@ # # You should have received a copy of the GNU Lesser General Public License # along with Iris. If not, see . -"""Unit tests for :class:`iris.analysis._regrid._CurvilinearRegridder`.""" +"""Unit tests for :class:`iris.analysis._regrid.CurvilinearRegridder`.""" from __future__ import (absolute_import, division, print_function) from six.moves import (filter, input, map, range, zip) # noqa diff --git a/lib/iris/tests/unit/analysis/test_PointInCell.py b/lib/iris/tests/unit/analysis/test_PointInCell.py index 0a44267785..2d4272712a 100644 --- a/lib/iris/tests/unit/analysis/test_PointInCell.py +++ b/lib/iris/tests/unit/analysis/test_PointInCell.py @@ -1,4 +1,4 @@ -# (C) British Crown Copyright 2015, Met Office +# (C) British Crown Copyright 2015 - 2019, Met Office # # This file is part of Iris. #