diff --git a/docs/iris/src/whatsnew/contributions_3.0.0/newfeature_2019-Dec-20_cache_area_weights.txt b/docs/iris/src/whatsnew/contributions_3.0.0/newfeature_2019-Dec-20_cache_area_weights.txt new file mode 100644 index 0000000000..b0cce56005 --- /dev/null +++ b/docs/iris/src/whatsnew/contributions_3.0.0/newfeature_2019-Dec-20_cache_area_weights.txt @@ -0,0 +1,2 @@ +* The area weights used when performing area weighted regridding are now + cached. diff --git a/lib/iris/analysis/_area_weighted.py b/lib/iris/analysis/_area_weighted.py index a7c0e41d4b..06f44dc951 100644 --- a/lib/iris/analysis/_area_weighted.py +++ b/lib/iris/analysis/_area_weighted.py @@ -41,36 +41,37 @@ def __init__(self, src_grid_cube, target_grid_cube, mdtol=1): .. Note:: - Both sourge and target cubes must have an XY grid defined by + Both source and target cubes must have an XY grid defined by separate X and Y dimensions with dimension coordinates. All of the XY dimension coordinates must also be bounded, and have the same cooordinate system. """ - # Snapshot the state of the cubes to ensure that the regridder is - # impervious to external changes to the original source cubes. + # Snapshot the state of the source cube to ensure that the regridder is + # impervious to external changes to the original cubes. self._src_grid = snapshot_grid(src_grid_cube) - self._target_grid = snapshot_grid(target_grid_cube) + # Missing data tolerance. if not (0 <= mdtol <= 1): msg = "Value for mdtol must be in range 0 - 1, got {}." raise ValueError(msg.format(mdtol)) self._mdtol = mdtol - # The need for an actual Cube is an implementation quirk caused by the - # current usage of the experimental regrid function. - self._target_grid_cube_cache = None - - @property - def _target_grid_cube(self): - if self._target_grid_cube_cache is None: - x, y = self._target_grid - data = np.empty((y.points.size, x.points.size)) - cube = iris.cube.Cube(data) - cube.add_dim_coord(y, 0) - cube.add_dim_coord(x, 1) - self._target_grid_cube_cache = cube - return self._target_grid_cube_cache + # Store regridding information + _regrid_info = eregrid._regrid_area_weighted_rectilinear_src_and_grid__prepare( + src_grid_cube, target_grid_cube + ) + ( + src_x, + src_y, + src_x_dim, + src_y_dim, + self.grid_x, + self.grid_y, + self.meshgrid_x, + self.meshgrid_y, + self.weights_info, + ) = _regrid_info def __call__(self, cube): """ @@ -92,11 +93,25 @@ def __call__(self, cube): area-weighted regridding. """ - if get_xy_dim_coords(cube) != self._src_grid: + src_x, src_y = get_xy_dim_coords(cube) + if (src_x, src_y) != self._src_grid: raise ValueError( "The given cube is not defined on the same " "source grid as this regridder." ) - return eregrid.regrid_area_weighted_rectilinear_src_and_grid( - cube, self._target_grid_cube, mdtol=self._mdtol + src_x_dim = cube.coord_dims(src_x)[0] + src_y_dim = cube.coord_dims(src_y)[0] + _regrid_info = ( + src_x, + src_y, + src_x_dim, + src_y_dim, + self.grid_x, + self.grid_y, + self.meshgrid_x, + self.meshgrid_y, + self.weights_info, + ) + return eregrid._regrid_area_weighted_rectilinear_src_and_grid__perform( + cube, _regrid_info, mdtol=self._mdtol ) diff --git a/lib/iris/analysis/_regrid.py b/lib/iris/analysis/_regrid.py index 19ed909fdc..0670c073ae 100644 --- a/lib/iris/analysis/_regrid.py +++ b/lib/iris/analysis/_regrid.py @@ -881,9 +881,9 @@ def _create_cube( def copy_coords(src_coords, add_method): for coord in src_coords: dims = src.coord_dims(coord) - if coord is src_x_coord: + if coord == src_x_coord: coord = grid_x_coord - elif coord is src_y_coord: + elif coord == src_y_coord: coord = grid_y_coord elif x_dim in dims or y_dim in dims: continue diff --git a/lib/iris/experimental/regrid.py b/lib/iris/experimental/regrid.py index 35d47a7a0f..e79e03cf72 100644 --- a/lib/iris/experimental/regrid.py +++ b/lib/iris/experimental/regrid.py @@ -402,20 +402,7 @@ def _weighted_mean_with_mdtol(data, weights, axis=None, mdtol=0): return res -def _regrid_area_weighted_array( - src_data, - x_dim, - y_dim, - src_x_bounds, - src_y_bounds, - grid_x_bounds, - grid_y_bounds, - grid_x_decreasing, - grid_y_decreasing, - area_func, - circular=False, - mdtol=0, -): +def _regrid_area_weighted_array(src_data, x_dim, y_dim, weights_info, mdtol=0): """ Regrid the given data from its source grid to a new grid using an area weighted mean to determine the resulting data values. @@ -434,30 +421,12 @@ def _regrid_area_weighted_array( The X dimension within `src_data`. * y_dim: The Y dimension within `src_data`. - * src_x_bounds: - A NumPy array of bounds along the X axis defining the source grid. - * src_y_bounds: - A NumPy array of bounds along the Y axis defining the source grid. - * grid_x_bounds: - A NumPy array of bounds along the X axis defining the new grid. - * grid_y_bounds: - A NumPy array of bounds along the Y axis defining the new grid. - * grid_x_decreasing: - Boolean indicating whether the X coordinate of the new grid is - in descending order. - * grid_y_decreasing: - Boolean indicating whether the Y coordinate of the new grid is - in descending order. - * area_func: - A function that returns an (p, q) array of weights given an (p, 2) - shaped array of Y bounds and an (q, 2) shaped array of X bounds. + * weights_info: + The area weights information to be used for area-weighted + regridding. Kwargs: - * circular: - A boolean indicating whether the `src_x_bounds` are periodic. Default - is False. - * mdtol: Tolerance of missing data. The value returned in each element of the returned array will be masked if the fraction of missing data exceeds @@ -473,116 +442,6 @@ def _regrid_area_weighted_array( grid. """ - - def _calculate_regrid_area_weighted_weights( - src_x_bounds, - src_y_bounds, - grid_x_bounds, - grid_y_bounds, - grid_x_decreasing, - grid_y_decreasing, - area_func, - circular=False, - ): - """ - Compute the area weights used for area-weighted regridding. - - """ - # Determine which grid bounds are within src extent. - y_within_bounds = _within_bounds( - src_y_bounds, grid_y_bounds, grid_y_decreasing - ) - x_within_bounds = _within_bounds( - src_x_bounds, grid_x_bounds, grid_x_decreasing - ) - - # Cache which src_bounds are within grid bounds - cached_x_bounds = [] - cached_x_indices = [] - max_x_indices = 0 - for (x_0, x_1) in grid_x_bounds: - if grid_x_decreasing: - x_0, x_1 = x_1, x_0 - x_bounds, x_indices = _cropped_bounds(src_x_bounds, x_0, x_1) - cached_x_bounds.append(x_bounds) - cached_x_indices.append(x_indices) - # Keep record of the largest slice - if isinstance(x_indices, slice): - x_indices_size = np.sum(x_indices.stop - x_indices.start) - else: # is tuple of indices - x_indices_size = len(x_indices) - if x_indices_size > max_x_indices: - max_x_indices = x_indices_size - - # Cache which y src_bounds areas and weights are within grid bounds - cached_y_indices = [] - cached_weights = [] - max_y_indices = 0 - for j, (y_0, y_1) in enumerate(grid_y_bounds): - # Reverse lower and upper if dest grid is decreasing. - if grid_y_decreasing: - y_0, y_1 = y_1, y_0 - y_bounds, y_indices = _cropped_bounds(src_y_bounds, y_0, y_1) - cached_y_indices.append(y_indices) - # Keep record of the largest slice - if isinstance(y_indices, slice): - y_indices_size = np.sum(y_indices.stop - y_indices.start) - else: # is tuple of indices - y_indices_size = len(y_indices) - if y_indices_size > max_y_indices: - max_y_indices = y_indices_size - - weights_i = [] - for i, (x_0, x_1) in enumerate(grid_x_bounds): - # Reverse lower and upper if dest grid is decreasing. - if grid_x_decreasing: - x_0, x_1 = x_1, x_0 - x_bounds = cached_x_bounds[i] - x_indices = cached_x_indices[i] - - # Determine whether element i, j overlaps with src and hence - # an area weight should be computed. - # If x_0 > x_1 then we want [0]->x_1 and x_0->[0] + mod in the case - # of wrapped longitudes. However if the src grid is not global - # (i.e. circular) this new cell would include a region outside of - # the extent of the src grid and thus the weight is therefore - # invalid. - outside_extent = x_0 > x_1 and not circular - if ( - outside_extent - or not y_within_bounds[j] - or not x_within_bounds[i] - ): - weights = False - else: - # Calculate weights based on areas of cropped bounds. - if isinstance(x_indices, tuple) and isinstance( - y_indices, tuple - ): - raise RuntimeError( - "Cannot handle split bounds " "in both x and y." - ) - weights = area_func(y_bounds, x_bounds) - weights_i.append(weights) - cached_weights.append(weights_i) - return ( - tuple(cached_x_indices), - tuple(cached_y_indices), - max_x_indices, - max_y_indices, - tuple(cached_weights), - ) - - weights_info = _calculate_regrid_area_weighted_weights( - src_x_bounds, - src_y_bounds, - grid_x_bounds, - grid_y_bounds, - grid_x_decreasing, - grid_y_decreasing, - area_func, - circular, - ) ( cached_x_indices, cached_y_indices, @@ -590,11 +449,6 @@ def _calculate_regrid_area_weighted_weights( max_y_indices, cached_weights, ) = weights_info - # Delete variables that are not needed and would not be available - # if _calculate_regrid_area_weighted_weights was refactored further - del src_x_bounds, src_y_bounds, grid_x_bounds, grid_y_bounds - del grid_x_decreasing, grid_y_decreasing - del area_func, circular # Ensure we have x_dim and y_dim. x_dim_orig = x_dim @@ -778,6 +632,25 @@ def regrid_area_weighted_rectilinear_src_and_grid( Returns: A new :class:`iris.cube.Cube` instance. + """ + regrid_info = _regrid_area_weighted_rectilinear_src_and_grid__prepare( + src_cube, grid_cube + ) + result = _regrid_area_weighted_rectilinear_src_and_grid__perform( + src_cube, regrid_info, mdtol + ) + return result + + +def _regrid_area_weighted_rectilinear_src_and_grid__prepare( + src_cube, grid_cube +): + """ + First (setup) part of 'regrid_area_weighted_rectilinear_src_and_grid'. + + Check inputs and calculate related info. The 'regrid info' returned + can be re-used over many 2d slices. + """ # Get the 1d monotonic (or scalar) src and grid coordinates. src_x, src_y = _get_xy_coords(src_cube) @@ -853,6 +726,9 @@ def regrid_area_weighted_rectilinear_src_and_grid( grid_x_bounds = _get_bounds_in_units(grid_x, x_units, dtype) grid_y_bounds = _get_bounds_in_units(grid_y, y_units, dtype) + # Create 2d meshgrids as required by _create_cube func. + meshgrid_x, meshgrid_y = _meshgrid(grid_x.points, grid_y.points) + # Determine whether target grid bounds are decreasing. This must # be determined prior to wrap_lons being called. grid_x_decreasing = grid_x_bounds[-1, 0] < grid_x_bounds[0, 0] @@ -881,11 +757,135 @@ def regrid_area_weighted_rectilinear_src_and_grid( else: area_func = _cartesian_area - # Calculate new data array for regridded cube. - new_data = _regrid_area_weighted_array( - src_cube.data, - src_x_dim, - src_y_dim, + def _calculate_regrid_area_weighted_weights( + src_x_bounds, + src_y_bounds, + grid_x_bounds, + grid_y_bounds, + grid_x_decreasing, + grid_y_decreasing, + area_func, + circular=False, + ): + """ + Compute the area weights used for area-weighted regridding. + + Args: + + * src_x_bounds: + A NumPy array of bounds along the X axis defining the source grid. + * src_y_bounds: + A NumPy array of bounds along the Y axis defining the source grid. + * grid_x_bounds: + A NumPy array of bounds along the X axis defining the new grid. + * grid_y_bounds: + A NumPy array of bounds along the Y axis defining the new grid. + * grid_x_decreasing: + Boolean indicating whether the X coordinate of the new grid is + in descending order. + * grid_y_decreasing: + Boolean indicating whether the Y coordinate of the new grid is + in descending order. + * area_func: + A function that returns an (p, q) array of weights given an (p, 2) + shaped array of Y bounds and an (q, 2) shaped array of X bounds. + + Kwargs: + + * circular: + A boolean indicating whether the `src_x_bounds` are periodic. + Default is False. + + Returns: + The area weights to be used for area-weighted regridding. + + """ + # Determine which grid bounds are within src extent. + y_within_bounds = _within_bounds( + src_y_bounds, grid_y_bounds, grid_y_decreasing + ) + x_within_bounds = _within_bounds( + src_x_bounds, grid_x_bounds, grid_x_decreasing + ) + + # Cache which src_bounds are within grid bounds + cached_x_bounds = [] + cached_x_indices = [] + max_x_indices = 0 + for (x_0, x_1) in grid_x_bounds: + if grid_x_decreasing: + x_0, x_1 = x_1, x_0 + x_bounds, x_indices = _cropped_bounds(src_x_bounds, x_0, x_1) + cached_x_bounds.append(x_bounds) + cached_x_indices.append(x_indices) + # Keep record of the largest slice + if isinstance(x_indices, slice): + x_indices_size = np.sum(x_indices.stop - x_indices.start) + else: # is tuple of indices + x_indices_size = len(x_indices) + if x_indices_size > max_x_indices: + max_x_indices = x_indices_size + + # Cache which y src_bounds areas and weights are within grid bounds + cached_y_indices = [] + cached_weights = [] + max_y_indices = 0 + for j, (y_0, y_1) in enumerate(grid_y_bounds): + # Reverse lower and upper if dest grid is decreasing. + if grid_y_decreasing: + y_0, y_1 = y_1, y_0 + y_bounds, y_indices = _cropped_bounds(src_y_bounds, y_0, y_1) + cached_y_indices.append(y_indices) + # Keep record of the largest slice + if isinstance(y_indices, slice): + y_indices_size = np.sum(y_indices.stop - y_indices.start) + else: # is tuple of indices + y_indices_size = len(y_indices) + if y_indices_size > max_y_indices: + max_y_indices = y_indices_size + + weights_i = [] + for i, (x_0, x_1) in enumerate(grid_x_bounds): + # Reverse lower and upper if dest grid is decreasing. + if grid_x_decreasing: + x_0, x_1 = x_1, x_0 + x_bounds = cached_x_bounds[i] + x_indices = cached_x_indices[i] + + # Determine whether element i, j overlaps with src and hence + # an area weight should be computed. + # If x_0 > x_1 then we want [0]->x_1 and x_0->[0] + mod in the case + # of wrapped longitudes. However if the src grid is not global + # (i.e. circular) this new cell would include a region outside of + # the extent of the src grid and thus the weight is therefore + # invalid. + outside_extent = x_0 > x_1 and not circular + if ( + outside_extent + or not y_within_bounds[j] + or not x_within_bounds[i] + ): + weights = False + else: + # Calculate weights based on areas of cropped bounds. + if isinstance(x_indices, tuple) and isinstance( + y_indices, tuple + ): + raise RuntimeError( + "Cannot handle split bounds " "in both x and y." + ) + weights = area_func(y_bounds, x_bounds) + weights_i.append(weights) + cached_weights.append(weights_i) + return ( + tuple(cached_x_indices), + tuple(cached_y_indices), + max_x_indices, + max_y_indices, + tuple(cached_weights), + ) + + weights_info = _calculate_regrid_area_weighted_weights( src_x_bounds, src_y_bounds, grid_x_bounds, @@ -894,12 +894,48 @@ def regrid_area_weighted_rectilinear_src_and_grid( grid_y_decreasing, area_func, circular, - mdtol, + ) + + return ( + src_x, + src_y, + src_x_dim, + src_y_dim, + grid_x, + grid_y, + meshgrid_x, + meshgrid_y, + weights_info, + ) + + +def _regrid_area_weighted_rectilinear_src_and_grid__perform( + src_cube, regrid_info, mdtol +): + """ + Second (regrid) part of 'regrid_area_weighted_rectilinear_src_and_grid'. + + Perform the prepared regrid calculation on a single 2d cube. + + """ + ( + src_x, + src_y, + src_x_dim, + src_y_dim, + grid_x, + grid_y, + meshgrid_x, + meshgrid_y, + weights_info, + ) = regrid_info + + # Calculate new data array for regridded cube. + new_data = _regrid_area_weighted_array( + src_cube.data, src_x_dim, src_y_dim, weights_info, mdtol, ) # Wrap up the data as a Cube. - # Create 2d meshgrids as required by _create_cube func. - meshgrid_x, meshgrid_y = _meshgrid(grid_x.points, grid_y.points) regrid_callback = RectilinearRegridder._regrid new_cube = RectilinearRegridder._create_cube( new_data, diff --git a/lib/iris/tests/unit/analysis/area_weighted/test_AreaWeightedRegridder.py b/lib/iris/tests/unit/analysis/area_weighted/test_AreaWeightedRegridder.py index a335044133..18312737fe 100644 --- a/lib/iris/tests/unit/analysis/area_weighted/test_AreaWeightedRegridder.py +++ b/lib/iris/tests/unit/analysis/area_weighted/test_AreaWeightedRegridder.py @@ -11,6 +11,7 @@ # Import iris.tests first so that some things can be initialised before # importing anything else. import iris.tests as tests +import iris.experimental.regrid as eregrid from unittest import mock @@ -30,11 +31,13 @@ def cube(self, x, y): lon = DimCoord(x, "longitude", units="degrees") cube.add_dim_coord(lat, 0) cube.add_dim_coord(lon, 1) + cube.coord("latitude").guess_bounds() + cube.coord("longitude").guess_bounds() return cube def grids(self): src = self.cube(np.linspace(20, 30, 3), np.linspace(10, 25, 4)) - target = self.cube(np.linspace(6, 18, 8), np.linspace(11, 22, 9)) + target = self.cube(np.linspace(22, 28, 8), np.linspace(11, 22, 9)) return src, target def extract_grid(self, cube): @@ -42,34 +45,51 @@ def extract_grid(self, cube): def check_mdtol(self, mdtol=None): src_grid, target_grid = self.grids() - if mdtol is None: - regridder = AreaWeightedRegridder(src_grid, target_grid) - mdtol = 1 - else: - regridder = AreaWeightedRegridder( - src_grid, target_grid, mdtol=mdtol - ) - - # Make a new cube to regrid with different data so we can - # distinguish between regridding the original src grid - # definition cube and the cube passed to the regridder. - src = src_grid.copy() - src.data += 10 - + # Get _regrid_info result + _regrid_info = eregrid._regrid_area_weighted_rectilinear_src_and_grid__prepare( + src_grid, target_grid + ) + self.assertEqual(len(_regrid_info), 9) with mock.patch( "iris.experimental.regrid." - "regrid_area_weighted_rectilinear_src_and_grid", - return_value=mock.sentinel.result, - ) as regrid: - result = regridder(src) - - self.assertEqual(regrid.call_count, 1) - _, args, kwargs = regrid.mock_calls[0] - - self.assertEqual(args[0], src) + "_regrid_area_weighted_rectilinear_src_and_grid__prepare", + return_value=_regrid_info, + ) as prepare: + with mock.patch( + "iris.experimental.regrid." + "_regrid_area_weighted_rectilinear_src_and_grid__perform", + return_value=mock.sentinel.result, + ) as perform: + # Setup the regridder + if mdtol is None: + regridder = AreaWeightedRegridder(src_grid, target_grid) + mdtol = 1 + else: + regridder = AreaWeightedRegridder( + src_grid, target_grid, mdtol=mdtol + ) + # Now regrid the source cube + src = src_grid + result = regridder(src) + + # Make a new cube to regrid with different data so we can + # distinguish between regridding the original src grid + # definition cube and the cube passed to the regridder. + src = src_grid.copy() + src.data += 10 + result = regridder(src) + + # Prepare: + self.assertEqual(prepare.call_count, 1) + _, args, kwargs = prepare.mock_calls[0] self.assertEqual( self.extract_grid(args[1]), self.extract_grid(target_grid) ) + + # Perform: + self.assertEqual(perform.call_count, 2) + _, args, kwargs = perform.mock_calls[1] + self.assertEqual(args[0], src) self.assertEqual(kwargs, {"mdtol": mdtol}) self.assertIs(result, mock.sentinel.result) @@ -102,6 +122,129 @@ def test_mismatched_src_coord_systems(self): with self.assertRaises(ValueError): AreaWeightedRegridder(src, target) + def test_src_and_target_are_the_same(self): + src = self.cube(np.linspace(20, 30, 3), np.linspace(10, 25, 4)) + target = self.cube(np.linspace(20, 30, 3), np.linspace(10, 25, 4)) + regridder = AreaWeightedRegridder(src, target) + result = regridder(src) + self.assertArrayAllClose(result.data, target.data) + + def test_multiple_src_on_same_grid(self): + coord_names = ["latitude", "longitude"] + src1 = self.cube(np.linspace(20, 32, 4), np.linspace(10, 22, 4)) + src2 = self.cube(np.linspace(20, 32, 4), np.linspace(10, 22, 4)) + src2.data *= 4 + self.assertArrayEqual(src1.data * 4, src2.data) + for name in coord_names: + # Remove coords system and units so it is no longer spherical. + src1.coord(name).coord_system = None + src1.coord(name).units = None + src2.coord(name).coord_system = None + src2.coord(name).units = None + + target = self.cube(np.linspace(20, 32, 2), np.linspace(10, 22, 2)) + # Ensure the bounds of the target cover the same range as the + # source. + target_lat_bounds = np.column_stack( + ( + src1.coord("latitude").bounds[[0, 1], [0, 1]], + src1.coord("latitude").bounds[[2, 3], [0, 1]], + ) + ) + target.coord("latitude").bounds = target_lat_bounds + target_lon_bounds = np.column_stack( + ( + src1.coord("longitude").bounds[[0, 1], [0, 1]], + src1.coord("longitude").bounds[[2, 3], [0, 1]], + ) + ) + target.coord("longitude").bounds = target_lon_bounds + for name in coord_names: + # Remove coords system and units so it is no longer spherical. + target.coord(name).coord_system = None + target.coord(name).units = None + + regridder = AreaWeightedRegridder(src1, target) + result1 = regridder(src1) + result2 = regridder(src2) + + reference1 = self.cube(np.linspace(20, 32, 2), np.linspace(10, 22, 2)) + reference1.data = np.array( + [ + [np.mean(src1.data[0:2, 0:2]), np.mean(src1.data[0:2, 2:4])], + [np.mean(src1.data[2:4, 0:2]), np.mean(src1.data[2:4, 2:4])], + ] + ) + reference1.coord("latitude").bounds = target_lat_bounds + reference1.coord("longitude").bounds = target_lon_bounds + + reference2 = self.cube(np.linspace(20, 32, 2), np.linspace(10, 22, 2)) + reference2.data = np.array( + [ + [np.mean(src2.data[0:2, 0:2]), np.mean(src2.data[0:2, 2:4])], + [np.mean(src2.data[2:4, 0:2]), np.mean(src2.data[2:4, 2:4])], + ] + ) + reference2.coord("latitude").bounds = target_lat_bounds + reference2.coord("longitude").bounds = target_lon_bounds + + for name in coord_names: + # Remove coords system and units so it is no longer spherical. + reference1.coord(name).coord_system = None + reference1.coord(name).units = None + reference2.coord(name).coord_system = None + reference2.coord(name).units = None + + # Compare the cubes rather than just the data. + self.assertEqual(result1, reference1) + self.assertEqual(result2, reference2) + + def test_src_data_different_dims(self): + src, target = self.grids() + regridder = AreaWeightedRegridder(src, target) + result = regridder(src) + expected_mean, expected_std = 4.772097735195653, 2.211698479817678 + self.assertArrayShapeStats(result, (9, 8), expected_mean, expected_std) + # New source cube with additional "levels" dimension + # Each level has identical x-y data so the mean and std stats remain + # identical when x, y and z dims are reordered + levels = DimCoord(np.arange(5), "model_level_number") + lat = src.coord("latitude") + lon = src.coord("longitude") + data = np.repeat(src.data[np.newaxis, ...], 5, axis=0) + src = Cube(data) + src.add_dim_coord(levels, 0) + src.add_dim_coord(lat, 1) + src.add_dim_coord(lon, 2) + result = regridder(src) + self.assertArrayShapeStats( + result, (5, 9, 8), expected_mean, expected_std + ) + # Check data with dims in different order + # Reshape src so that the coords are ordered [x, z, y], + # the mean and std statistics should be the same + data = np.moveaxis(src.data.copy(), 2, 0) + src = Cube(data) + src.add_dim_coord(lon, 0) + src.add_dim_coord(levels, 1) + src.add_dim_coord(lat, 2) + result = regridder(src) + self.assertArrayShapeStats( + result, (8, 5, 9), expected_mean, expected_std + ) + # Check data with dims in different order + # Reshape src so that the coords are ordered [y, x, z], + # the mean and std statistics should be the same + data = np.moveaxis(src.data.copy(), 2, 0) + src = Cube(data) + src.add_dim_coord(lat, 0) + src.add_dim_coord(lon, 1) + src.add_dim_coord(levels, 2) + result = regridder(src) + self.assertArrayShapeStats( + result, (9, 8, 5), expected_mean, expected_std + ) + if __name__ == "__main__": tests.main()