Repository: dcs4cop/xcube
Pull request #282 (2020-03-03 17:54:28)

Fix rectify_dataset()
In this PR:

- renamed `rectify_dataset()` kwarg `is_y_axis_inverted` into `is_y_reversed`
- fixed `rectify_dataset()` to correctly deal with `is_y_reversed` kwark if `tile_size` (= dask mode) is given too
- added new unit-tests for tiles + y axis reversed

@AliceBalfanz, please note, this PR may introduce new problems:

If we call `rectify_dataset(dataset, is_y_reversed=True, ...)` as done by `XYInputProcessor`, the chunking in the y-dimension changes, e.g. y-chunksizes (512, 512, 512, 512, 43) turns into (43, 512, 512, 512, 512). Zarr doesn't like this, you may have to call `chunk_dataset()` in `gen_cube()` which may lead to performance loss. In the future, we should only generate cubes with `is_y_reversed=False`!!!


Patch:
diff --git a/xcube/core/gen/iproc.py b/xcube/core/gen/iproc.py
index 977b8ea2..8be29445 100644
--- a/xcube/core/gen/iproc.py
+++ b/xcube/core/gen/iproc.py
@@ -25,6 +25,7 @@ from typing import Any, Collection, Dict, Mapping, Optional, Tuple, Union
 
 import numpy as np
 import xarray as xr
+
 from xcube.constants import CRS_WKT_EPSG_4326
 from xcube.constants import EXTENSION_POINT_INPUT_PROCESSORS
 from xcube.core.geocoding import GeoCoding
@@ -326,7 +327,7 @@ class XYInputProcessor(InputProcessor, metaclass=ABCMeta):
                                       compute_subset=False,
                                       geo_coding=geo_coding,
                                       output_geom=output_geom,
-                                      is_y_axis_inverted=True)
+                                      is_y_reversed=True)
 
             if dataset is not None and geo_coding.is_geo_crs and geo_coding.xy_names != ('lon', 'lat'):
                 dataset = dataset.rename({geo_coding.x_name: 'lon', geo_coding.y_name: 'lat'})
diff --git a/xcube/core/imgeom.py b/xcube/core/imgeom.py
index 39c60674..7e7d4f65 100644
--- a/xcube/core/imgeom.py
+++ b/xcube/core/imgeom.py
@@ -178,7 +178,7 @@ class ImageGeom:
     def coord_vars(self,
                    xy_names: Tuple[str, str],
                    is_lon_normalized: bool = False,
-                   is_y_axis_inverted: bool = False) -> Mapping[str, xr.DataArray]:
+                   is_y_reversed: bool = False) -> Mapping[str, xr.DataArray]:
         x_name, y_name = xy_names
         x_attrs, y_attrs = (_LON_ATTRS, _LAT_ATTRS) if self.is_geo_crs else ({}, {})
         w, h = self.size
@@ -199,7 +199,7 @@ class ImageGeom:
         y_bnds_0_data = np.linspace(y1, y2 - res, h)
         y_bnds_1_data = np.linspace(y1 + res, y2, h)
 
-        if is_y_axis_inverted:
+        if is_y_reversed:
             y_data = y_data[::-1]
             y_bnds_1_data, y_bnds_0_data = y_bnds_0_data[::-1], y_bnds_1_data[::-1]
 
diff --git a/xcube/core/rectify.py b/xcube/core/rectify.py
index f505a445..45de920d 100644
--- a/xcube/core/rectify.py
+++ b/xcube/core/rectify.py
@@ -37,7 +37,7 @@ def rectify_dataset(dataset: xr.Dataset,
                     geo_coding: GeoCoding = None,
                     xy_names: Tuple[str, str] = None,
                     output_geom: ImageGeom = None,
-                    is_y_axis_inverted: bool = False,
+                    is_y_reversed: bool = False,
                     tile_size: Union[int, Tuple[int, int]] = None,
                     output_ij_names: Tuple[str, str] = None,
                     load_xy: bool = False,
@@ -66,7 +66,7 @@ def rectify_dataset(dataset: xr.Dataset,
     :param xy_names: Optional tuple of the x- and y-coordinate variables in *dataset*. Ignored if *geo_coding* is given.
     :param output_geom: Optional output geometry. If not given, output geometry will be computed
         to spatially fit *dataset* and to retain its spatial resolution.
-    :param is_y_axis_inverted: Whether the y-axis labels in the output should be in inverse order.
+    :param is_y_reversed: Whether the y-axis labels in the output should be in reverse order.
     :param tile_size: Optional tile size for the output.
     :param output_ij_names: If given, a tuple of variable names in which to store the computed source pixel
         coordinates in the returned output.
@@ -120,13 +120,18 @@ def rectify_dataset(dataset: xr.Dataset,
 
     dst_src_ij_array = get_dst_src_ij_images(src_geo_coding,
                                              output_geom,
-                                             is_y_axis_inverted,
                                              uv_delta)
 
+    if is_y_reversed:
+        # Note that the following reverse operation may change any y-axis' chunking.
+        # This is the case if the y-chunksize does not integer-divide y-size, e.g.
+        # y-chunksizes (512, 512, 512, 273) will change into (273, 512, 512, 512).
+        dst_src_ij_array = dst_src_ij_array[:, ::-1]
+
     dst_dims = src_geo_coding.xy_names[::-1]
     dst_ds_coords = output_geom.coord_vars(xy_names=src_geo_coding.xy_names,
                                            is_lon_normalized=src_geo_coding.is_lon_normalized,
-                                           is_y_axis_inverted=is_y_axis_inverted)
+                                           is_y_reversed=is_y_reversed)
     dst_vars = dict()
     for src_var_name, src_var in src_vars.items():
         if load_vars:
@@ -194,7 +199,6 @@ def _is_2d_var(var: xr.DataArray, two_d_coord_var: xr.DataArray) -> bool:
 
 def _compute_ij_images_xarray_numpy(src_geo_coding: GeoCoding,
                                     output_geom: ImageGeom,
-                                    is_dst_y_axis_inverted: bool,
                                     uv_delta: float) -> np.ndarray:
     """Compute numpy.ndarray destination image with source pixel i,j coords from xarray.DataArray x,y sources """
     dst_width = output_geom.width
@@ -213,14 +217,11 @@ def _compute_ij_images_xarray_numpy(src_geo_coding: GeoCoding,
                                       dst_y_min,
                                       dst_xy_res,
                                       uv_delta)
-    if is_dst_y_axis_inverted:
-        dst_src_ij_images = dst_src_ij_images[:, ::-1, :]
     return dst_src_ij_images
 
 
 def _compute_ij_images_xarray_dask(src_geo_coding: GeoCoding,
                                    output_geom: ImageGeom,
-                                   is_dst_y_axis_inverted: bool,
                                    uv_delta: float) -> da.Array:
     """Compute dask.array.Array destination image with source pixel i,j coords from xarray.DataArray x,y sources """
     dst_width = output_geom.width
@@ -236,13 +237,13 @@ def _compute_ij_images_xarray_dask(src_geo_coding: GeoCoding,
     # Compute an empirical xy_border as a function of the number of tiles, because the more tiles we have
     # the smaller the destination xy-bboxes and the higher the risk to not find any source ij-bbox for
     # a given xy-bbox.
+    # xy_border will not be larger than half of the coverage of a tile.
     #
-    num_tiles_x = (dst_width + dst_tile_width - 1) / dst_tile_width
-    num_tiles_y = (dst_height + dst_tile_height - 1) / dst_tile_height
+    num_tiles_x = dst_width / dst_tile_width
+    num_tiles_y = dst_height / dst_tile_height
     max_num_tiles = max(num_tiles_x, num_tiles_y)
     xy_border = min(2 * max_num_tiles * dst_xy_res,
-                    min(0.5 * (dst_x_max - dst_x_min),
-                        0.5 * (dst_y_max - dst_y_min)))
+                    min(0.5 * (dst_x_max - dst_x_min), 0.5 * (dst_y_max - dst_y_min)))
 
     dst_xy_bboxes = output_geom.xy_bboxes
     src_ij_bboxes = src_geo_coding.ij_bboxes(dst_xy_bboxes, xy_border=xy_border, ij_border=1)
@@ -264,7 +265,6 @@ def _compute_ij_images_xarray_dask(src_geo_coding: GeoCoding,
                                        dst_x_min,
                                        dst_y_min,
                                        dst_xy_res,
-                                       is_dst_y_axis_inverted,
                                        uv_delta
                                    ),
                                    name='ij_pixels',
@@ -281,7 +281,6 @@ def _compute_ij_images_xarray_dask_block(dtype: np.dtype,
                                          dst_x_min: float,
                                          dst_y_min: float,
                                          dst_xy_res: float,
-                                         is_dst_y_axis_inverted: bool,
                                          uv_delta: float) -> np.ndarray:
     """Compute dask.array.Array destination block with source pixel i,j coords from xarray.DataArray x,y sources """
     dst_src_ij_block = np.full(block_shape, np.nan, dtype=dtype)
@@ -301,8 +300,6 @@ def _compute_ij_images_xarray_dask_block(dtype: np.dtype,
                                         dst_y_min + dst_y_slice_start * dst_xy_res,
                                         dst_xy_res,
                                         uv_delta)
-    if is_dst_y_axis_inverted:
-        dst_src_ij_block = dst_src_ij_block[:, ::-1, :]
     return dst_src_ij_block
 
 


Test patch:
diff --git a/test/core/test_imgeom.py b/test/core/test_imgeom.py
index e009671a..55db4cf5 100644
--- a/test/core/test_imgeom.py
+++ b/test/core/test_imgeom.py
@@ -160,10 +160,10 @@ class ImageGeomTest(SourceDatasetMixin, unittest.TestCase):
                                     (1250., 1260.),
                                 ))
 
-    def test_coord_vars_y_axis_inverted(self):
+    def test_coord_vars_y_reversed(self):
         image_geom = ImageGeom(size=(10, 6), x_min=-2600.0, y_min=1200.0, xy_res=10.0)
 
-        cv = image_geom.coord_vars(xy_names=('x', 'y'), is_y_axis_inverted=True)
+        cv = image_geom.coord_vars(xy_names=('x', 'y'), is_y_reversed=True)
         self._assert_coord_vars(cv,
                                 (10, 6),
                                 ('x', 'y'),
diff --git a/test/core/test_rectify.py b/test/core/test_rectify.py
index c0d7af23..d35e217b 100644
--- a/test/core/test_rectify.py
+++ b/test/core/test_rectify.py
@@ -97,12 +97,12 @@ class RectifyDatasetTest(SourceDatasetMixin, unittest.TestCase):
                                                 dtype=lat.dtype))
         np.testing.assert_almost_equal(rad.values, self.expected_rad_13x13(rad.dtype))
 
-    def test_rectify_2x2_to_13x13_inverse_y_axis(self):
+    def test_rectify_2x2_to_13x13_y_reversed(self):
         src_ds = self.new_source_dataset()
 
         output_geom = ImageGeom(size=(13, 13), x_min=-0.25, y_min=49.75, xy_res=0.5)
 
-        dst_ds = rectify_dataset(src_ds, output_geom=output_geom, is_y_axis_inverted=True)
+        dst_ds = rectify_dataset(src_ds, output_geom=output_geom, is_y_reversed=True)
         lon, lat, rad = self._assert_shape_and_dim(dst_ds, (13, 13))
         np.testing.assert_almost_equal(lon.values,
                                        np.array([0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5., 5.5, 6.],
@@ -113,6 +113,22 @@ class RectifyDatasetTest(SourceDatasetMixin, unittest.TestCase):
                                                 dtype=lat.dtype))
         np.testing.assert_almost_equal(rad.values, self.expected_rad_13x13(rad.dtype)[::-1])
 
+    def test_rectify_2x2_to_13x13_y_reversed_dask_5x5(self):
+        src_ds = self.new_source_dataset()
+
+        output_geom = ImageGeom(size=(13, 13), x_min=-0.25, y_min=49.75, xy_res=0.5, tile_size=5)
+
+        dst_ds = rectify_dataset(src_ds, output_geom=output_geom, is_y_reversed=True)
+        lon, lat, rad = self._assert_shape_and_dim(dst_ds, (13, 13), chunks=((3, 5, 5), (5, 5, 3)))
+        np.testing.assert_almost_equal(lon.values,
+                                       np.array([0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5., 5.5, 6.],
+                                                dtype=lon.dtype))
+        np.testing.assert_almost_equal(lat.values,
+                                       np.array([56., 55.5, 55., 54.5, 54.0, 53.5, 53.0, 52.5, 52.0, 51.5, 51.,
+                                                 50.5, 50.],
+                                                dtype=lat.dtype))
+        np.testing.assert_almost_equal(rad.values, self.expected_rad_13x13(rad.dtype)[::-1])
+
     def test_rectify_2x2_to_13x13_dask_7x7(self):
         src_ds = self.new_source_dataset()
         output_geom = ImageGeom(size=(13, 13), x_min=-0.25, y_min=49.75, xy_res=0.5, tile_size=7)
@@ -141,7 +157,6 @@ class RectifyDatasetTest(SourceDatasetMixin, unittest.TestCase):
                                        np.array([50., 50.5, 51.0, 51.5, 52.0, 52.5, 53.0, 53.5, 54.0, 54.5, 55.,
                                                  55.5, 56.],
                                                 dtype=lat.dtype))
-        print(rad.values)
         np.testing.assert_almost_equal(rad.values, self.expected_rad_13x13(rad.dtype))
 
     def test_rectify_2x2_to_13x13_dask_3x13(self):
@@ -384,11 +399,21 @@ class RectifySentinel2DatasetTest(SourceDatasetMixin, unittest.TestCase):
              [nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]])
 
         dst_ds = rectify_dataset(src_ds, tile_size=None)
+        self.assertEqual(None, dst_ds.rrs_665.chunks)
         np.testing.assert_almost_equal(dst_ds.rrs_665.values, expected_data)
 
         dst_ds = rectify_dataset(src_ds, tile_size=5)
+        self.assertEqual(((5, 2), (5, 5)), dst_ds.rrs_665.chunks)
         np.testing.assert_almost_equal(dst_ds.rrs_665.values, expected_data)
 
+        dst_ds = rectify_dataset(src_ds, tile_size=None, is_y_reversed=True)
+        self.assertEqual(None, dst_ds.rrs_665.chunks)
+        np.testing.assert_almost_equal(dst_ds.rrs_665.values, expected_data[::-1])
+
+        dst_ds = rectify_dataset(src_ds, tile_size=5, is_y_reversed=True)
+        self.assertEqual(((2, 5), (5, 5)), dst_ds.rrs_665.chunks)
+        np.testing.assert_almost_equal(dst_ds.rrs_665.values, expected_data[::-1])
+
 
 def create_s2plus_dataset():
     x = xr.DataArray([310005., 310015., 310025., 310035., 310045.], dims=["x"],