ENH: Working zonal mean linear regridding for circular sources or with use of extrapolation

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:

lib/iris/tests/integration/test_regridding.py

@@ -109,5 +109,116 @@ def test_nearest(self):
         self.assertArrayShapeStats(res, (1, 6, 3, 4), 315.890808, 11.000724)
 
 
+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)))
+        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))
+
+
+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(
+            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)
+
+        # 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
+
+    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.
+        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):
+        # 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 (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)
+        self.assertTrue(np.isnan(res.data).all())
+
+
 if __name__ == "__main__":
     tests.main()