From 82961f0f6aa85b924be822a3f70c7734ecd5a5de Mon Sep 17 00:00:00 2001 From: cpelley Date: Thu, 12 Jul 2018 15:24:33 +0100 Subject: [PATCH 1/3] Working zonal mean from circular/extrapolated source to target --- lib/iris/analysis/_regrid.py | 5 +- lib/iris/tests/integration/test_regridding.py | 97 +++++++++++++++++++ 2 files changed, 100 insertions(+), 2 deletions(-) diff --git a/lib/iris/analysis/_regrid.py b/lib/iris/analysis/_regrid.py index eb7a16d075..71304a5dde 100644 --- a/lib/iris/analysis/_regrid.py +++ b/lib/iris/analysis/_regrid.py @@ -231,8 +231,9 @@ def _regrid(src_data, x_dim, y_dim, data = np.empty(shape, dtype=dtype) # The interpolation class requires monotonically increasing - # coordinates, so flip the coordinate(s) and data if the aren't. - reverse_x = src_x_coord.points[0] > src_x_coord.points[1] + # coordinates, so flip the coordinate(s) and data if they aren't. + reverse_x = (src_x_coord.points[0] > src_x_coord.points[1] if + src_x_coord.points.size > 1 else False) reverse_y = src_y_coord.points[0] > src_y_coord.points[1] flip_index = [slice(None)] * src_data.ndim if reverse_x: diff --git a/lib/iris/tests/integration/test_regridding.py b/lib/iris/tests/integration/test_regridding.py index ef521aa395..dd4aecdabd 100644 --- a/lib/iris/tests/integration/test_regridding.py +++ b/lib/iris/tests/integration/test_regridding.py @@ -109,5 +109,102 @@ def test_nearest(self): self.assertArrayShapeStats(res, (1, 6, 3, 4), 315.890808, 11.000724) +def plot_cube(source): + import iris.plot as iplt + iplt.plt.switch_backend('TkAgg') + iplt.pcolormesh(source) + iplt.plt.gca().coastlines() + iplt.plt.colorbar(orientation='horizontal') + iplt.show() + + +@tests.skip_data +class TestZonalMean(tests.IrisTest): + def setUp(self): + np.random.seed(0) + self.src = iris.cube.Cube(np.random.random_integers(0, 10, (140, 1))) + s_crs = iris.coord_systems.GeogCS(6371229.0) + sy_coord = iris.coords.DimCoord( + np.linspace(-90, 90, 140), standard_name='latitude', + units='degrees', coord_system=s_crs) + sx_coord = iris.coords.DimCoord( + -180, bounds=[-180, 180], standard_name='longitude', + units='degrees', circular=True, coord_system=s_crs) + self.src.add_dim_coord(sy_coord, 0) + self.src.add_dim_coord(sx_coord, 1) + + def test_linear_same_crs_global(self): + # Regrid the zonal mean onto an identical coordinate system target, but + # on a different set of longitudes - which should result in no change. + points = [-150, -90, -30, 30, 90, 150] + bounds = [[-180, -120], [-120, -60], [-60, 0], [0, 60], [60, 120], + [120, 180]] + sx_coord = self.src.coord(axis='x') + sy_coord = self.src.coord(axis='y') + x_coord = sx_coord.copy(points, bounds=bounds) + grid = iris.cube.Cube( + np.zeros([sy_coord.points.size, x_coord.points.size])) + grid.add_dim_coord(sy_coord, 0) + grid.add_dim_coord(x_coord, 1) + + res = self.src.regrid(grid, iris.analysis.Linear()) + + # Ensure data remains unchanged. + # (the same along each column) + self.assertTrue( + np.array([(res.data[:, 0]-res.data[:, i]).max() for i in + range(1, res.shape[1])]).max() < 1e-10) + self.assertArrayAlmostEqual(res.data[:, 0], self.src.data.reshape(-1)) + + def test_linear_rotated_regional(self): + sx_coord = self.src.coord(axis='x') + sy_coord = self.src.coord(axis='y') + grid_crs = iris.coord_systems.RotatedGeogCS( + 37.5, 177.5, ellipsoid=iris.coord_systems.GeogCS(6371229.0)) + grid_x = sx_coord.copy(np.linspace(350, 370, 100)) + grid_x.circular = False + grid_x.coord_system = grid_crs + grid_y = sy_coord.copy(np.linspace(-10, 10, 100)) + grid_y.coord_system = grid_crs + grid = iris.cube.Cube( + np.zeros([grid_y.points.size, grid_x.points.size])) + grid.add_dim_coord(grid_y, 0) + grid.add_dim_coord(grid_x, 1) + + res = self.src.regrid(grid, iris.analysis.Linear()) + + res_no_extrapolation = self.src.regrid( + grid, iris.analysis.Linear(extrapolation_mode='nan')) + + self.src.coord(axis='x').circular = False + res_non_circular = self.src.regrid(grid, iris.analysis.Linear()) + + res_non_circular_no_extrapolation = self.src.regrid( + grid, iris.analysis.Linear(extrapolation_mode='nan')) + + def munge(src_cube): + # Munge the source (duplicate source latitudes) so that we can + # utilise linear regridding in the ordinary case (that is, as + # though it was not a zonal mean). + src_cube2 = src_cube.copy() + src_cube2.coord(axis='x').points = -90 + src_cube2.coord(axis='x').bounds = [-180, 0] + src_cube.coord(axis='x').points = 90 + src_cube.coord(axis='x').bounds = [0, 180] + src_cubes = iris.cube.CubeList([src_cube, src_cube2]) + return src_cubes.concatenate_cube() + tar = munge(self.src).regrid(grid, iris.analysis.Linear()) + + self.assertArrayAlmostEqual(res.data, tar.data) + self.assertArrayAlmostEqual(res_no_extrapolation.data, tar.data) + self.assertArrayAlmostEqual(res_non_circular.data, tar.data) + + # Let's capture the fact that no extrapolation and a non circular + # source will not lead to the same results (given that Linear is a + # points based approach). + self.assertFalse(np.allclose(res_non_circular_no_extrapolation.data, + tar.data)) + + if __name__ == "__main__": tests.main() From 5b3a55b72a68a26366601fcce7c3cd7a719bbe0b Mon Sep 17 00:00:00 2001 From: cpelley Date: Fri, 13 Jul 2018 09:35:37 +0100 Subject: [PATCH 2/3] MAINT: Refactor of zonal mean testing --- lib/iris/tests/integration/test_regridding.py | 93 ++++++++++--------- 1 file changed, 51 insertions(+), 42 deletions(-) diff --git a/lib/iris/tests/integration/test_regridding.py b/lib/iris/tests/integration/test_regridding.py index dd4aecdabd..fe72f0b6b6 100644 --- a/lib/iris/tests/integration/test_regridding.py +++ b/lib/iris/tests/integration/test_regridding.py @@ -109,17 +109,7 @@ def test_nearest(self): self.assertArrayShapeStats(res, (1, 6, 3, 4), 315.890808, 11.000724) -def plot_cube(source): - import iris.plot as iplt - iplt.plt.switch_backend('TkAgg') - iplt.pcolormesh(source) - iplt.plt.gca().coastlines() - iplt.plt.colorbar(orientation='horizontal') - iplt.show() - - -@tests.skip_data -class TestZonalMean(tests.IrisTest): +class TestZonalMean_global(tests.IrisTest): def setUp(self): np.random.seed(0) self.src = iris.cube.Cube(np.random.random_integers(0, 10, (140, 1))) @@ -156,7 +146,13 @@ def test_linear_same_crs_global(self): range(1, res.shape[1])]).max() < 1e-10) self.assertArrayAlmostEqual(res.data[:, 0], self.src.data.reshape(-1)) - def test_linear_rotated_regional(self): + +class TestZonalMean_regional(TestZonalMean_global, tests.IrisTest): + def setUp(self): + super(TestZonalMean_regional, self).setUp() + + # Define a target grid and a target result (what we expect the + # regridder to return). sx_coord = self.src.coord(axis='x') sy_coord = self.src.coord(axis='y') grid_crs = iris.coord_systems.RotatedGeogCS( @@ -171,39 +167,52 @@ def test_linear_rotated_regional(self): grid.add_dim_coord(grid_y, 0) grid.add_dim_coord(grid_x, 1) - res = self.src.regrid(grid, iris.analysis.Linear()) + self.tar = self.zonal_mean_as_multi_column(self.src).regrid( + grid, iris.analysis.Linear()) + self.grid = grid - res_no_extrapolation = self.src.regrid( - grid, iris.analysis.Linear(extrapolation_mode='nan')) + def zonal_mean_as_multi_column(self, src_cube): + # Munge the source (duplicate source latitudes) so that we can + # utilise linear regridding as a conventional problem (that is, to + # duplicate columns so that it is no longer a zonal mean problem). + src_cube2 = src_cube.copy() + src_cube2.coord(axis='x').points = -90 + src_cube2.coord(axis='x').bounds = [-180, 0] + src_cube.coord(axis='x').points = 90 + src_cube.coord(axis='x').bounds = [0, 180] + src_cubes = iris.cube.CubeList([src_cube, src_cube2]) + return src_cubes.concatenate_cube() + def test_linear_rotated_regional(self): + # Ensure that zonal mean source data is linearly interpolated onto a + # high resolution target. We turn a zonal mean source into a multi + # column in order to derive the target values we expect the regridder + # to return (munge function below). + regridder = iris.analysis.Linear() + res = self.src.regrid(self.grid, regridder) + self.assertArrayAlmostEqual(res.data, self.tar.data) + + def test_linear_rotated_regional_no_extrapolation(self): + regridder = iris.analysis.Linear(extrapolation_mode='nan') + res = self.src.regrid(self.grid, regridder) + self.assertArrayAlmostEqual(res.data, self.tar.data) + + def test_linear_rotated_regional_not_circular(self): + regridder = iris.analysis.Linear() + self.src.coord(axis='x').circular = False + res = self.src.regrid(self.grid, regridder) + self.assertArrayAlmostEqual(res.data, self.tar.data) + + def test_linear_rotated_regional_no_extrapolation_not_circular(self): + # This technically is a test that confirm how zonal mean actually works + # in so far as, that extrapolation and circular source handling is the + # means by which these usecases are supported. + # In the case where the source is neither using extrapolation and is + # not circular, then 'nan' values will result. + regridder = iris.analysis.Linear(extrapolation_mode='nan') self.src.coord(axis='x').circular = False - res_non_circular = self.src.regrid(grid, iris.analysis.Linear()) - - res_non_circular_no_extrapolation = self.src.regrid( - grid, iris.analysis.Linear(extrapolation_mode='nan')) - - def munge(src_cube): - # Munge the source (duplicate source latitudes) so that we can - # utilise linear regridding in the ordinary case (that is, as - # though it was not a zonal mean). - src_cube2 = src_cube.copy() - src_cube2.coord(axis='x').points = -90 - src_cube2.coord(axis='x').bounds = [-180, 0] - src_cube.coord(axis='x').points = 90 - src_cube.coord(axis='x').bounds = [0, 180] - src_cubes = iris.cube.CubeList([src_cube, src_cube2]) - return src_cubes.concatenate_cube() - tar = munge(self.src).regrid(grid, iris.analysis.Linear()) - - self.assertArrayAlmostEqual(res.data, tar.data) - self.assertArrayAlmostEqual(res_no_extrapolation.data, tar.data) - self.assertArrayAlmostEqual(res_non_circular.data, tar.data) - - # Let's capture the fact that no extrapolation and a non circular - # source will not lead to the same results (given that Linear is a - # points based approach). - self.assertFalse(np.allclose(res_non_circular_no_extrapolation.data, - tar.data)) + res = self.src.regrid(self.grid, regridder) + self.assertTrue(np.isnan(res.data).all()) if __name__ == "__main__": From bc0d389c11977e4b48d9312b5e6e6ffb4025ba42 Mon Sep 17 00:00:00 2001 From: cpelley Date: Fri, 13 Jul 2018 11:46:20 +0100 Subject: [PATCH 3/3] MAINT: Documentation changes from review --- lib/iris/tests/integration/test_regridding.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/lib/iris/tests/integration/test_regridding.py b/lib/iris/tests/integration/test_regridding.py index fe72f0b6b6..47d155f3f2 100644 --- a/lib/iris/tests/integration/test_regridding.py +++ b/lib/iris/tests/integration/test_regridding.py @@ -167,6 +167,9 @@ def setUp(self): grid.add_dim_coord(grid_y, 0) grid.add_dim_coord(grid_x, 1) + # The target result is derived by regridding a multi-column version of + # the source to the target (i.e. turning a zonal mean regrid into a + # conventional regrid). self.tar = self.zonal_mean_as_multi_column(self.src).regrid( grid, iris.analysis.Linear()) self.grid = grid @@ -185,30 +188,32 @@ def zonal_mean_as_multi_column(self, src_cube): def test_linear_rotated_regional(self): # Ensure that zonal mean source data is linearly interpolated onto a - # high resolution target. We turn a zonal mean source into a multi - # column in order to derive the target values we expect the regridder - # to return (munge function below). + # high resolution target. regridder = iris.analysis.Linear() res = self.src.regrid(self.grid, regridder) self.assertArrayAlmostEqual(res.data, self.tar.data) def test_linear_rotated_regional_no_extrapolation(self): + # Capture the case where our source remains circular but we don't use + # extrapolation. regridder = iris.analysis.Linear(extrapolation_mode='nan') res = self.src.regrid(self.grid, regridder) self.assertArrayAlmostEqual(res.data, self.tar.data) def test_linear_rotated_regional_not_circular(self): + # Capture the case where our source is not circular but we utilise + # extrapolation. regridder = iris.analysis.Linear() self.src.coord(axis='x').circular = False res = self.src.regrid(self.grid, regridder) self.assertArrayAlmostEqual(res.data, self.tar.data) def test_linear_rotated_regional_no_extrapolation_not_circular(self): - # This technically is a test that confirm how zonal mean actually works - # in so far as, that extrapolation and circular source handling is the - # means by which these usecases are supported. + # Confirm how zonal mean actually works in so far as, that + # extrapolation and circular source handling is the means by which + # these usecases are supported. # In the case where the source is neither using extrapolation and is - # not circular, then 'nan' values will result. + # not circular, then 'nan' values will result (as we would expect). regridder = iris.analysis.Linear(extrapolation_mode='nan') self.src.coord(axis='x').circular = False res = self.src.regrid(self.grid, regridder)