From 77f40fc7313068df08eb580562ee3c69da95058b Mon Sep 17 00:00:00 2001 From: Emma Hogan Date: Tue, 17 Dec 2019 12:35:26 +0000 Subject: [PATCH 1/8] Split regrid_area_weighted_rectilinear_src_and_grid into prepare and perform (#3601) * Split regrid_area_weighted_rectilinear_src_and_grid into prepare and perform --- lib/iris/experimental/regrid.py | 72 ++++++++++++++++++++++++++++++++- 1 file changed, 70 insertions(+), 2 deletions(-) diff --git a/lib/iris/experimental/regrid.py b/lib/iris/experimental/regrid.py index 35d47a7a0f..66717dd8f1 100644 --- a/lib/iris/experimental/regrid.py +++ b/lib/iris/experimental/regrid.py @@ -778,6 +778,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 +872,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,6 +903,54 @@ def regrid_area_weighted_rectilinear_src_and_grid( else: area_func = _cartesian_area + return ( + src_x, + src_y, + src_x_dim, + src_y_dim, + src_x_bounds, + src_y_bounds, + grid_x, + grid_y, + grid_x_bounds, + grid_y_bounds, + grid_x_decreasing, + grid_y_decreasing, + meshgrid_x, + meshgrid_y, + area_func, + circular, + ) + + +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, + src_x_bounds, + src_y_bounds, + grid_x, + grid_y, + grid_x_bounds, + grid_y_bounds, + grid_x_decreasing, + grid_y_decreasing, + meshgrid_x, + meshgrid_y, + area_func, + circular, + ) = regrid_info + # Calculate new data array for regridded cube. new_data = _regrid_area_weighted_array( src_cube.data, @@ -898,8 +968,6 @@ def regrid_area_weighted_rectilinear_src_and_grid( ) # 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, From 7f02e17a5797d43829392a06a92ff867cef224a3 Mon Sep 17 00:00:00 2001 From: Emma Hogan Date: Tue, 17 Dec 2019 15:24:57 +0000 Subject: [PATCH 2/8] Test for AreaWeightedRegridder (#3604) * Add tests for AreaWeightedRegridder --- .../test_AreaWeightedRegridder.py | 62 +++++++++++++++++++ 1 file changed, 62 insertions(+) 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..5ab24ee89d 100644 --- a/lib/iris/tests/unit/analysis/area_weighted/test_AreaWeightedRegridder.py +++ b/lib/iris/tests/unit/analysis/area_weighted/test_AreaWeightedRegridder.py @@ -102,6 +102,68 @@ 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)) + for name in ["latitude", "longitude"]: + src.coord(name).guess_bounds() + target.coord(name).guess_bounds() + 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 + # Ensure contiguous bounds exists. + src1.coord(name).guess_bounds() + src2.coord(name).guess_bounds() + + target = self.cube(np.linspace(20, 32, 2), np.linspace(10, 22, 2)) + 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 + # Ensure the bounds of the target cover the same range as the + # source. + target.coord(name).bounds = np.column_stack( + ( + src1.coord(name).bounds[[0, 1], [0, 1]], + src1.coord(name).bounds[[2, 3], [0, 1]], + ) + ) + + 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])], + ] + ) + 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])], + ] + ) + + self.assertArrayEqual(result1.data, reference1.data) + self.assertArrayEqual(result2.data, reference2.data) + if __name__ == "__main__": tests.main() From 6356ece4c15421423bc42f59125af9df8edf4bce Mon Sep 17 00:00:00 2001 From: Emma Hogan Date: Fri, 20 Dec 2019 09:27:32 +0000 Subject: [PATCH 3/8] PI-2472: Refactor AreaWeightedRegridder (#3606) * Add prepare and perform to the AreaWeightedRegridder --- lib/iris/analysis/_area_weighted.py | 31 +++++++-- .../test_AreaWeightedRegridder.py | 66 +++++++++++++++---- 2 files changed, 79 insertions(+), 18 deletions(-) diff --git a/lib/iris/analysis/_area_weighted.py b/lib/iris/analysis/_area_weighted.py index a7c0e41d4b..a6ab8586cb 100644 --- a/lib/iris/analysis/_area_weighted.py +++ b/lib/iris/analysis/_area_weighted.py @@ -41,16 +41,22 @@ 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. + # impervious to external changes to the original cubes. self._src_grid = snapshot_grid(src_grid_cube) self._target_grid = snapshot_grid(target_grid_cube) + + # Store the x_dim and y_dim of the source cube. + x, y = get_xy_dim_coords(src_grid_cube) + self._src_x_dim = src_grid_cube.coord_dims(x) + self._src_y_dim = src_grid_cube.coord_dims(y) + # Missing data tolerance. if not (0 <= mdtol <= 1): msg = "Value for mdtol must be in range 0 - 1, got {}." @@ -61,6 +67,10 @@ def __init__(self, src_grid_cube, target_grid_cube, mdtol=1): # current usage of the experimental regrid function. self._target_grid_cube_cache = None + self._regrid_info = eregrid._regrid_area_weighted_rectilinear_src_and_grid__prepare( + src_grid_cube, self._target_grid_cube + ) + @property def _target_grid_cube(self): if self._target_grid_cube_cache is None: @@ -92,11 +102,22 @@ def __call__(self, cube): area-weighted regridding. """ - if get_xy_dim_coords(cube) != self._src_grid: + if get_xy_dim_coords(cube) != self._src_grid or not ( + _xy_data_dims_are_equal(cube, self._src_x_dim, self._src_y_dim) + ): 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 + return eregrid._regrid_area_weighted_rectilinear_src_and_grid__perform( + cube, self._regrid_info, mdtol=self._mdtol ) + + +def _xy_data_dims_are_equal(cube, x_dim, y_dim): + """ + Return whether the data dimensions of the x and y coordinates on the + the cube are equal to the values ``x_dim`` and ``y_dim``, respectively. + """ + x1, y1 = get_xy_dim_coords(cube) + return cube.coord_dims(x1) == x_dim and cube.coord_dims(y1) == y_dim 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 5ab24ee89d..b4dc426c6d 100644 --- a/lib/iris/tests/unit/analysis/area_weighted/test_AreaWeightedRegridder.py +++ b/lib/iris/tests/unit/analysis/area_weighted/test_AreaWeightedRegridder.py @@ -42,13 +42,17 @@ 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 - ) + with mock.patch( + "iris.experimental.regrid." + "_regrid_area_weighted_rectilinear_src_and_grid__prepare" + ) as prepare: + 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 @@ -58,18 +62,22 @@ def check_mdtol(self, mdtol=None): with mock.patch( "iris.experimental.regrid." - "regrid_area_weighted_rectilinear_src_and_grid", + "_regrid_area_weighted_rectilinear_src_and_grid__perform", return_value=mock.sentinel.result, - ) as regrid: + ) as perform: result = regridder(src) - self.assertEqual(regrid.call_count, 1) - _, args, kwargs = regrid.mock_calls[0] - - self.assertEqual(args[0], 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, 1) + _, args, kwargs = perform.mock_calls[0] + self.assertEqual(args[0], src) self.assertEqual(kwargs, {"mdtol": mdtol}) self.assertIs(result, mock.sentinel.result) @@ -164,6 +172,38 @@ def test_multiple_src_on_same_grid(self): self.assertArrayEqual(result1.data, reference1.data) self.assertArrayEqual(result2.data, reference2.data) + def test_mismatched_data_dims(self): + coord_names = ["latitude", "longitude"] + x = np.linspace(20, 32, 4) + y = np.linspace(10, 22, 4) + src1 = self.cube(x, y) + + data = np.arange(len(y) * len(x)).reshape(len(x), len(y)) + src2 = Cube(data) + lat = DimCoord(y, "latitude", units="degrees") + lon = DimCoord(x, "longitude", units="degrees") + # Add dim coords in opposite order to self.cube. + src2.add_dim_coord(lat, 1) + src2.add_dim_coord(lon, 0) + for name in coord_names: + # Ensure contiguous bounds exists. + src1.coord(name).guess_bounds() + src2.coord(name).guess_bounds() + + target = self.cube(np.linspace(20, 32, 2), np.linspace(10, 22, 2)) + for name in coord_names: + # Ensure the bounds of the target cover the same range as the + # source. + target.coord(name).bounds = np.column_stack( + ( + src1.coord(name).bounds[[0, 1], [0, 1]], + src1.coord(name).bounds[[2, 3], [0, 1]], + ) + ) + + regridder = AreaWeightedRegridder(src1, target) + self.assertRaises(ValueError, regridder, src2) + if __name__ == "__main__": tests.main() From 75d01598caadce7969f4f669e83ba911ada6f0e5 Mon Sep 17 00:00:00 2001 From: Emma Hogan Date: Fri, 20 Dec 2019 14:12:44 +0000 Subject: [PATCH 4/8] Move calculation of area weights from perform to prepare (#3617) --- lib/iris/experimental/regrid.py | 230 ++++++++++++++++---------------- 1 file changed, 115 insertions(+), 115 deletions(-) diff --git a/lib/iris/experimental/regrid.py b/lib/iris/experimental/regrid.py index 66717dd8f1..c5715fcb25 100644 --- a/lib/iris/experimental/regrid.py +++ b/lib/iris/experimental/regrid.py @@ -413,6 +413,7 @@ def _regrid_area_weighted_array( grid_x_decreasing, grid_y_decreasing, area_func, + weights_info, circular=False, mdtol=0, ): @@ -451,6 +452,7 @@ def _regrid_area_weighted_array( * 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 to be used for area-weighted regridding. Kwargs: @@ -473,116 +475,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 +482,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 @@ -903,6 +790,116 @@ def _regrid_area_weighted_rectilinear_src_and_grid__prepare( else: area_func = _cartesian_area + 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, + ) + return ( src_x, src_y, @@ -919,6 +916,7 @@ def _regrid_area_weighted_rectilinear_src_and_grid__prepare( meshgrid_x, meshgrid_y, area_func, + weights_info, circular, ) @@ -948,6 +946,7 @@ def _regrid_area_weighted_rectilinear_src_and_grid__perform( meshgrid_x, meshgrid_y, area_func, + weights_info, circular, ) = regrid_info @@ -963,6 +962,7 @@ def _regrid_area_weighted_rectilinear_src_and_grid__perform( grid_x_decreasing, grid_y_decreasing, area_func, + weights_info, circular, mdtol, ) From 37d507c309ca49a4659d0847cb344987d0c39236 Mon Sep 17 00:00:00 2001 From: Emma Hogan Date: Fri, 20 Dec 2019 16:09:15 +0000 Subject: [PATCH 5/8] Add docstring and what's new (#3619) * Add docstring and what's new * Actually add docstring --- ...feature_2019-Dec-20_cache_area_weights.txt | 2 ++ lib/iris/experimental/regrid.py | 29 +++++++++++++++++++ 2 files changed, 31 insertions(+) create mode 100644 docs/iris/src/whatsnew/contributions_3.0.0/newfeature_2019-Dec-20_cache_area_weights.txt 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/experimental/regrid.py b/lib/iris/experimental/regrid.py index c5715fcb25..b2637a8fc6 100644 --- a/lib/iris/experimental/regrid.py +++ b/lib/iris/experimental/regrid.py @@ -803,6 +803,35 @@ def _calculate_regrid_area_weighted_weights( """ 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( From 841bb838891b14f847e393934d37bff228fe4390 Mon Sep 17 00:00:00 2001 From: Emma Hogan Date: Fri, 20 Dec 2019 18:12:18 +0000 Subject: [PATCH 6/8] Remove unused variables in regridding code (#3620) --- lib/iris/experimental/regrid.py | 71 +++------------------------------ 1 file changed, 5 insertions(+), 66 deletions(-) diff --git a/lib/iris/experimental/regrid.py b/lib/iris/experimental/regrid.py index b2637a8fc6..e79e03cf72 100644 --- a/lib/iris/experimental/regrid.py +++ b/lib/iris/experimental/regrid.py @@ -402,21 +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, - weights_info, - 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. @@ -435,31 +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 to be used for area-weighted regridding. + * 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 @@ -934,19 +901,11 @@ def _calculate_regrid_area_weighted_weights( src_y, src_x_dim, src_y_dim, - src_x_bounds, - src_y_bounds, grid_x, grid_y, - grid_x_bounds, - grid_y_bounds, - grid_x_decreasing, - grid_y_decreasing, meshgrid_x, meshgrid_y, - area_func, weights_info, - circular, ) @@ -964,36 +923,16 @@ def _regrid_area_weighted_rectilinear_src_and_grid__perform( src_y, src_x_dim, src_y_dim, - src_x_bounds, - src_y_bounds, grid_x, grid_y, - grid_x_bounds, - grid_y_bounds, - grid_x_decreasing, - grid_y_decreasing, meshgrid_x, meshgrid_y, - area_func, weights_info, - circular, ) = regrid_info # Calculate new data array for regridded cube. new_data = _regrid_area_weighted_array( - src_cube.data, - src_x_dim, - src_y_dim, - src_x_bounds, - src_y_bounds, - grid_x_bounds, - grid_y_bounds, - grid_x_decreasing, - grid_y_decreasing, - area_func, - weights_info, - circular, - mdtol, + src_cube.data, src_x_dim, src_y_dim, weights_info, mdtol, ) # Wrap up the data as a Cube. From c09c4698bdc5f7416bcddbd8e32da7107b7ba584 Mon Sep 17 00:00:00 2001 From: Emma Hogan Date: Tue, 24 Dec 2019 09:52:49 +0000 Subject: [PATCH 7/8] PI-2472: Fix issue related to comparison / equality in regridding code (#3624) Looks good. --- lib/iris/analysis/_regrid.py | 4 +- .../test_AreaWeightedRegridder.py | 41 ++++++++++++++----- 2 files changed, 33 insertions(+), 12 deletions(-) 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/tests/unit/analysis/area_weighted/test_AreaWeightedRegridder.py b/lib/iris/tests/unit/analysis/area_weighted/test_AreaWeightedRegridder.py index b4dc426c6d..1e494de259 100644 --- a/lib/iris/tests/unit/analysis/area_weighted/test_AreaWeightedRegridder.py +++ b/lib/iris/tests/unit/analysis/area_weighted/test_AreaWeightedRegridder.py @@ -137,18 +137,26 @@ def test_multiple_src_on_same_grid(self): src2.coord(name).guess_bounds() 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 - # Ensure the bounds of the target cover the same range as the - # source. - target.coord(name).bounds = np.column_stack( - ( - src1.coord(name).bounds[[0, 1], [0, 1]], - src1.coord(name).bounds[[2, 3], [0, 1]], - ) - ) regridder = AreaWeightedRegridder(src1, target) result1 = regridder(src1) @@ -161,6 +169,9 @@ def test_multiple_src_on_same_grid(self): [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( [ @@ -168,9 +179,19 @@ def test_multiple_src_on_same_grid(self): [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 - self.assertArrayEqual(result1.data, reference1.data) - self.assertArrayEqual(result2.data, reference2.data) + 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_mismatched_data_dims(self): coord_names = ["latitude", "longitude"] From fb5170fa28d403f89e85fb6c0b9c9d17d32e455f Mon Sep 17 00:00:00 2001 From: abooton Date: Tue, 7 Jan 2020 10:48:54 +0000 Subject: [PATCH 8/8] PI-2472: Update area weighted regridder to accept different dim data (#3625) * Update area weighted regridder to accept data with dims in different orders * Add extra unittest to test_AreaWeightedRegridder --- lib/iris/analysis/_area_weighted.py | 68 ++++----- .../test_AreaWeightedRegridder.py | 142 ++++++++++-------- 2 files changed, 112 insertions(+), 98 deletions(-) diff --git a/lib/iris/analysis/_area_weighted.py b/lib/iris/analysis/_area_weighted.py index a6ab8586cb..06f44dc951 100644 --- a/lib/iris/analysis/_area_weighted.py +++ b/lib/iris/analysis/_area_weighted.py @@ -47,15 +47,9 @@ def __init__(self, src_grid_cube, target_grid_cube, mdtol=1): the same cooordinate system. """ - # Snapshot the state of the cubes to ensure that the regridder is + # 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) - - # Store the x_dim and y_dim of the source cube. - x, y = get_xy_dim_coords(src_grid_cube) - self._src_x_dim = src_grid_cube.coord_dims(x) - self._src_y_dim = src_grid_cube.coord_dims(y) # Missing data tolerance. if not (0 <= mdtol <= 1): @@ -63,24 +57,21 @@ def __init__(self, src_grid_cube, target_grid_cube, mdtol=1): 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 - - self._regrid_info = eregrid._regrid_area_weighted_rectilinear_src_and_grid__prepare( - src_grid_cube, self._target_grid_cube + # Store regridding information + _regrid_info = eregrid._regrid_area_weighted_rectilinear_src_and_grid__prepare( + src_grid_cube, target_grid_cube ) - - @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 + ( + 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): """ @@ -102,22 +93,25 @@ def __call__(self, cube): area-weighted regridding. """ - if get_xy_dim_coords(cube) != self._src_grid or not ( - _xy_data_dims_are_equal(cube, self._src_x_dim, self._src_y_dim) - ): + 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." ) + 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, self._regrid_info, mdtol=self._mdtol + cube, _regrid_info, mdtol=self._mdtol ) - - -def _xy_data_dims_are_equal(cube, x_dim, y_dim): - """ - Return whether the data dimensions of the x and y coordinates on the - the cube are equal to the values ``x_dim`` and ``y_dim``, respectively. - """ - x1, y1 = get_xy_dim_coords(cube) - return cube.coord_dims(x1) == x_dim and cube.coord_dims(y1) == y_dim 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 1e494de259..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,30 +45,39 @@ def extract_grid(self, cube): def check_mdtol(self, mdtol=None): src_grid, target_grid = self.grids() + # 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__prepare" + "_regrid_area_weighted_rectilinear_src_and_grid__prepare", + return_value=_regrid_info, ) as prepare: - 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 - - with mock.patch( - "iris.experimental.regrid." - "_regrid_area_weighted_rectilinear_src_and_grid__perform", - return_value=mock.sentinel.result, - ) as perform: - result = regridder(src) + 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) @@ -75,8 +87,8 @@ def check_mdtol(self, mdtol=None): ) # Perform: - self.assertEqual(perform.call_count, 1) - _, args, kwargs = perform.mock_calls[0] + 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) @@ -113,9 +125,6 @@ def test_mismatched_src_coord_systems(self): 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)) - for name in ["latitude", "longitude"]: - src.coord(name).guess_bounds() - target.coord(name).guess_bounds() regridder = AreaWeightedRegridder(src, target) result = regridder(src) self.assertArrayAllClose(result.data, target.data) @@ -132,9 +141,6 @@ def test_multiple_src_on_same_grid(self): src1.coord(name).units = None src2.coord(name).coord_system = None src2.coord(name).units = None - # Ensure contiguous bounds exists. - src1.coord(name).guess_bounds() - src2.coord(name).guess_bounds() 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 @@ -193,37 +199,51 @@ def test_multiple_src_on_same_grid(self): self.assertEqual(result1, reference1) self.assertEqual(result2, reference2) - def test_mismatched_data_dims(self): - coord_names = ["latitude", "longitude"] - x = np.linspace(20, 32, 4) - y = np.linspace(10, 22, 4) - src1 = self.cube(x, y) - - data = np.arange(len(y) * len(x)).reshape(len(x), len(y)) - src2 = Cube(data) - lat = DimCoord(y, "latitude", units="degrees") - lon = DimCoord(x, "longitude", units="degrees") - # Add dim coords in opposite order to self.cube. - src2.add_dim_coord(lat, 1) - src2.add_dim_coord(lon, 0) - for name in coord_names: - # Ensure contiguous bounds exists. - src1.coord(name).guess_bounds() - src2.coord(name).guess_bounds() - - target = self.cube(np.linspace(20, 32, 2), np.linspace(10, 22, 2)) - for name in coord_names: - # Ensure the bounds of the target cover the same range as the - # source. - target.coord(name).bounds = np.column_stack( - ( - src1.coord(name).bounds[[0, 1], [0, 1]], - src1.coord(name).bounds[[2, 3], [0, 1]], - ) - ) - - regridder = AreaWeightedRegridder(src1, target) - self.assertRaises(ValueError, regridder, src2) + 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__":