Repository: dcs4cop/xcube (python)

Issue:
Resample .zarr Cube
**Describe the bug**
Error during the xcube resample for Sentinel-2 dataset.

**To Reproduce**
1. Generate the Sentinel-2 MSI irregular time cube using the `xcube gen`.

To generate it I took the netCDF files (tiles 31UET, 31UES, 31UDT, 31UDS).

**I use the config parameters:**
![TICKET01_01](https://user-images.githubusercontent.com/32440152/67864720-58e14c80-fb26-11e9-91db-e0fcaa0c152f.png)

The **output_size** and **output_region** I took from `xcube grid abox`:
![TICKET01_02](https://user-images.githubusercontent.com/32440152/67864773-77474800-fb26-11e9-8526-7855d69ab2b5.png)

**The .zarr cube was created with no errors.**
![ticket01_06](https://user-images.githubusercontent.com/32440152/67864968-c7bea580-fb26-11e9-87da-cf69c515be50.png)

**2. Run the `xcube resample` with the command line:**
![ticket01_03](https://user-images.githubusercontent.com/32440152/67865040-e9b82800-fb26-11e9-990c-8d9881cb53bd.png)

**HERE IS THE ERROR MESSAGE:**
![xcube_ticket01_resample_issue_with_chunk](https://user-images.githubusercontent.com/32440152/67864136-71049c00-fb25-11e9-8e33-5c156fd7f554.png)

**Run again the `xcube resample` but now with output as netCDF4:**
![ticket01_04](https://user-images.githubusercontent.com/32440152/67865109-0f453180-fb27-11e9-95e4-bbf64fc17b46.png)

**With output as netCDF4 the resample successfully finished.**
![ticket01_05](https://user-images.githubusercontent.com/32440152/67865248-461b4780-fb27-11e9-8c0d-54539ca9e651.png)

**Expected behavior**
I was expected to resample using as output the .zarr format (netCDF4 is too big and too slow).

**Screenshots**
Added during the text.

**Desktop (please complete the following information):**
 - OS: Linux Ubuntu 18.04
 - Browser: Firefox (not necessary)
 - Version: concerning the xcube version I updated the xcube and xcube-gen-bc before start (30/10/2019)

**Additional context**
I used the same `xcube resample` before for the Sentinel-3 (full North Sea) and it worked perfect to generate the 1D resample cube in .zarr format. **I think my problem is due the chunk option.**
What should I do to edit the chunk size to a *ideal match*.

Thank you for the help!!


Pull request description:
Forman 212 resample bug
* `xcube resample` now correctly re-chunks its output. By default, chunking of the 
  `time` dimension is set to one. (#212)
* refactored out new function `xcube.core.update.update_dataset_chunk_encoding()`

Closes #212



Patch:
diff --git a/CHANGES.md b/CHANGES.md
index 30b3bf37..fddbf17b 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -16,6 +16,10 @@
 
 * CLI commands execute much faster now when invoked with the `--help` and `--info` options.
 
+### Fixes
+
+* `xcube resample` now correctly re-chunks its output. By default, chunking of the 
+  `time` dimension is set to one. (#212)
 
 ### Incompatible changes
 
diff --git a/xcube/cli/resample.py b/xcube/cli/resample.py
index 58ade0a3..f73c273e 100644
--- a/xcube/cli/resample.py
+++ b/xcube/cli/resample.py
@@ -167,7 +167,7 @@ def _resample_in_time(input_path: str = None,
                       methods: Sequence[str] = (DEFAULT_RESAMPLING_METHOD,),
                       frequency: str = DEFAULT_RESAMPLING_FREQUENCY,
                       offset: str = None,
-                      base: str = DEFAULT_RESAMPLING_BASE,
+                      base: int = DEFAULT_RESAMPLING_BASE,
                       interp_kind: str = DEFAULT_INTERPOLATION_KIND,
                       tolerance: str = None,
                       dry_run: bool = False,
@@ -176,6 +176,7 @@ def _resample_in_time(input_path: str = None,
     from xcube.core.dsio import open_cube
     from xcube.core.dsio import write_cube
     from xcube.core.resample import resample_in_time
+    from xcube.core.update import update_dataset_chunk_encoding
 
     if not output_format:
         output_format = guess_dataset_format(output_path)
@@ -191,9 +192,15 @@ def _resample_in_time(input_path: str = None,
                                   base=base,
                                   interp_kind=interp_kind,
                                   tolerance=tolerance,
+                                  time_chunk_size=1,
                                   var_names=variables,
                                   metadata=metadata)
 
+        agg_ds = update_dataset_chunk_encoding(agg_ds,
+                                               chunk_sizes={},
+                                               format_name=output_format,
+                                               in_place=True)
+
         monitor(f'Writing resampled cube to {output_path!r}...')
         if not dry_run:
             write_cube(agg_ds, output_path, output_format, cube_asserted=True)
diff --git a/xcube/core/chunk.py b/xcube/core/chunk.py
index 7d74d653..327b2847 100644
--- a/xcube/core/chunk.py
+++ b/xcube/core/chunk.py
@@ -4,41 +4,24 @@ from typing import Dict, Tuple, Iterable
 import numpy as np
 import xarray as xr
 
-from xcube.constants import FORMAT_NAME_ZARR, FORMAT_NAME_NETCDF4
+from xcube.core.update import update_dataset_chunk_encoding
 
 
 def chunk_dataset(dataset: xr.Dataset,
                   chunk_sizes: Dict[str, int] = None,
                   format_name: str = None) -> xr.Dataset:
     """
-    Chunk dataset and update encodings for given format.
+    Chunk *dataset* using *chunk_sizes* and optionally update encodings for given *format_name*.
 
     :param dataset: input dataset
     :param chunk_sizes: mapping from dimension name to new chunk size
-    :param format_name: format, e.g. "zarr" or "netcdf4"
-    :return: the re-chunked dataset
+    :param format_name: optional format, e.g. "zarr" or "netcdf4"
+    :return: the (re)chunked dataset
     """
-    chunked_ds = dataset.chunk(chunks=chunk_sizes)
-
-    # Update encoding so writing of chunked_ds recognizes new chunks
-    chunk_sizes_attr_name = None
-    if format_name == FORMAT_NAME_ZARR:
-        chunk_sizes_attr_name = "chunks"
-    if format_name == FORMAT_NAME_NETCDF4:
-        chunk_sizes_attr_name = "chunksizes"
-    if chunk_sizes_attr_name:
-        for var_name in chunked_ds.variables:
-            var = chunked_ds[var_name]
-            if chunk_sizes:
-                sizes = tuple(chunk_sizes[dim_name] if dim_name in chunk_sizes
-                              else var.shape[var.dims.index(dim_name)]
-                              for dim_name in var.dims)
-                var.encoding.update({chunk_sizes_attr_name: sizes})
-            elif chunk_sizes_attr_name in var.encoding:
-                # Remove any explicit and wrong specification so writing will use Dask chunks (TBC!)
-                del var.encoding[chunk_sizes_attr_name]
-
-    return chunked_ds
+    dataset = dataset.chunk(chunks=chunk_sizes)
+    if format_name:
+        dataset = update_dataset_chunk_encoding(dataset, chunk_sizes=chunk_sizes, format_name=format_name)
+    return dataset
 
 
 def get_empty_dataset_chunks(dataset: xr.Dataset) -> Dict[str, Tuple[Tuple[int, ...]]]:
diff --git a/xcube/core/resample.py b/xcube/core/resample.py
index 5b995407..9f353b1a 100644
--- a/xcube/core/resample.py
+++ b/xcube/core/resample.py
@@ -23,32 +23,44 @@ from typing import Dict, Any, Sequence, Union
 
 import xarray as xr
 
+from xcube.core.schema import CubeSchema
 from xcube.core.select import select_vars
+from xcube.core.verify import assert_cube
 
 
 def resample_in_time(cube: xr.Dataset,
                      frequency: str,
                      method: Union[str, Sequence[str]],
                      offset=None,
-                     base: str = 0,
+                     base: int = 0,
                      tolerance=None,
                      interp_kind=None,
+                     time_chunk_size=None,
                      var_names: Sequence[str] = None,
-                     metadata: Dict[str, Any] = None):
+                     metadata: Dict[str, Any] = None,
+                     cube_asserted: bool = False) -> xr.Dataset:
     """
     Resample a xcube dataset in the time dimension.
 
     :param cube: The xcube dataset.
-    :param frequency: Resampling frequency.
+    :param frequency: Temporal aggregation frequency. Use format "<count><offset>"
+        "where <offset> is one of 'H', 'D', 'W', 'M', 'Q', 'Y'.
     :param method: Resampling method or sequence of resampling methods.
-    :param offset: Offset used to adjust the resampled time labels. Some pandas date offset strings are supported.
-    :param base: Resampling method.
+    :param offset: Offset used to adjust the resampled time labels.
+        Uses same syntax as *frequency*.
+    :param base: For frequencies that evenly subdivide 1 day, the "origin" of the
+        aggregated intervals. For example, for '24H' frequency, base could range from 0 through 23.
+    :param time_chunk_size: If not None, the chunk size to be used for the "time" dimension.
     :param var_names: Variable names to include.
     :param tolerance: Time tolerance for selective upsampling methods. Defaults to *frequency*.
     :param interp_kind: Kind of interpolation if *method* is 'interpolation'.
     :param metadata: Output metadata.
+    :param cube_asserted: If False, *cube* will be verified, otherwise it is expected to be a valid cube.
     :return: A new xcube dataset resampled in time.
     """
+    if not cube_asserted:
+        assert_cube(cube)
+
     if var_names:
         cube = select_vars(cube, var_names)
 
@@ -70,7 +82,8 @@ def resample_in_time(cube: xr.Dataset,
         resampling_method = getattr(resampler, method)
         kwargs = get_method_kwargs(method, frequency, interp_kind, tolerance)
         resampled_cube = resampling_method(**kwargs)
-        resampled_cube = resampled_cube.rename({var_name: f'{var_name}_{method}' for var_name in resampled_cube.data_vars})
+        resampled_cube = resampled_cube.rename(
+            {var_name: f'{var_name}_{method}' for var_name in resampled_cube.data_vars})
         resampled_cubes.append(resampled_cube)
 
     if len(resampled_cubes) == 1:
@@ -86,7 +99,14 @@ def resample_in_time(cube: xr.Dataset,
     # TODO: add other time_coverage_ attributes
     resampled_cube.attrs.update(time_coverage_start=time_coverage_start,
                                 time_coverage_end=time_coverage_end)
-    return resampled_cube
+
+    schema = CubeSchema.new(cube)
+    chunk_sizes = {schema.dims[i]: schema.chunks[i] for i in range(schema.ndim)}
+
+    if isinstance(time_chunk_size, int) and time_chunk_size >= 0:
+        chunk_sizes['time'] = time_chunk_size
+
+    return resampled_cube.chunk(chunk_sizes)
 
 
 def get_method_kwargs(method, frequency, interp_kind, tolerance):
diff --git a/xcube/core/update.py b/xcube/core/update.py
index f1b18be2..28bdb9a6 100644
--- a/xcube/core/update.py
+++ b/xcube/core/update.py
@@ -24,6 +24,8 @@ from typing import Any, Dict
 
 import xarray as xr
 
+from xcube.constants import FORMAT_NAME_NETCDF4
+from xcube.constants import FORMAT_NAME_ZARR
 from xcube.util.config import NameDictPairList
 
 _LON_ATTRS_DATA = ('lon', 'lon_bnds', 'degrees_east',
@@ -192,3 +194,48 @@ def update_dataset_var_attrs(dataset: xr.Dataset,
             var.attrs.update(var_attrs)
 
     return dataset
+
+
+def update_dataset_chunk_encoding(dataset: xr.Dataset,
+                                  chunk_sizes: Dict[str, int] = None,
+                                  format_name: str = None,
+                                  in_place: bool = False) -> xr.Dataset:
+    """
+    Update each variable's encoding in *dataset* with respect to *chunk_sizes*
+    so *dataset* is written in chunks for given *format_name*.
+
+    :param dataset: input dataset.
+    :param chunk_sizes: the chunk sizes to be used for the encoding.
+        If None, any chunking encoding is removed.
+    :param format_name: format name, e.g. "zarr" or "netcdf4".
+    :param in_place: If ``True``, *dataset* will be modified in place and returned.
+    """
+    if format_name == FORMAT_NAME_ZARR:
+        chunk_sizes_attr_name = "chunks"
+    elif format_name == FORMAT_NAME_NETCDF4:
+        chunk_sizes_attr_name = "chunksizes"
+    else:
+        return dataset
+    if not in_place:
+        dataset = dataset.copy()
+    for var_name in dataset.variables:
+        var = dataset[var_name]
+        if chunk_sizes is not None:
+            def get_size(i):
+                dim_name = var.dims[i]
+                size = chunk_sizes.get(dim_name)
+                if isinstance(size, int):
+                    return size
+                if var.chunks:
+                    size = var.chunks[i]
+                    if isinstance(size, int):
+                        return size
+                    if len(size):
+                        return size[0]
+                return var.shape[i]
+
+            var.encoding.update({chunk_sizes_attr_name: tuple(map(get_size, range(var.ndim)))})
+        elif chunk_sizes_attr_name in var.encoding:
+            # Remove any explicit and possibly unintended specification
+            del var.encoding[chunk_sizes_attr_name]
+    return dataset


Test patch:
diff --git a/test/core/test_chunk.py b/test/core/test_chunk.py
index 9127a9b4..8e267e21 100644
--- a/test/core/test_chunk.py
+++ b/test/core/test_chunk.py
@@ -29,6 +29,26 @@ class ChunkDatasetTest(unittest.TestCase):
         self.assertEqual({}, chunked_dataset.precipitation.encoding)
         self.assertEqual({}, chunked_dataset.temperature.encoding)
 
+        dataset = dataset.chunk(dict(time=2, lat=10, lon=20))
+
+        chunked_dataset = chunk_dataset(dataset,
+                                        chunk_sizes=None,
+                                        format_name="zarr")
+        self.assertEqual({}, chunked_dataset.precipitation.encoding)
+        self.assertEqual({}, chunked_dataset.temperature.encoding)
+
+        chunked_dataset = chunk_dataset(dataset,
+                                        chunk_sizes={},
+                                        format_name="zarr")
+        self.assertEqual({'chunks': (2, 10, 20)}, chunked_dataset.precipitation.encoding)
+        self.assertEqual({'chunks': (2, 10, 20)}, chunked_dataset.temperature.encoding)
+
+        chunked_dataset = chunk_dataset(dataset,
+                                        chunk_sizes=dict(time=1),
+                                        format_name="zarr")
+        self.assertEqual({'chunks': (1, 10, 20)}, chunked_dataset.precipitation.encoding)
+        self.assertEqual({'chunks': (1, 10, 20)}, chunked_dataset.temperature.encoding)
+
     def test_unchunk_dataset(self):
         dataset = new_test_dataset(["2010-01-01", "2010-01-02", "2010-01-03", "2010-01-04", "2010-01-05"],
                                    precipitation=0.4, temperature=275.2)
diff --git a/test/core/test_resample.py b/test/core/test_resample.py
index ddbc343f..256e8daa 100644
--- a/test/core/test_resample.py
+++ b/test/core/test_resample.py
@@ -4,47 +4,71 @@ import numpy as np
 import pandas as pd
 
 from test.sampledata import new_test_dataset
+from xcube.core.chunk import chunk_dataset
 from xcube.core.resample import resample_in_time
+from xcube.core.schema import CubeSchema
 
 
-class GenerateL3CubeTest(unittest.TestCase):
+class ResampleInTimeTest(unittest.TestCase):
 
-    def test_resample(self):
+    def setUp(self) -> None:
         num_times = 30
 
-        time, temperature, precipitation = zip(*[(('2017-07-0%s' if i < 9 else '2017-07-%s') % (i + 1),
-                                                  272 + 0.1 * i,
-                                                  120 - 0.2 * i) for i in range(num_times)])
-
-        ds1 = new_test_dataset(time, temperature=temperature, precipitation=precipitation)
-
-        ds2 = resample_in_time(ds1, '3D', ['min', 'max'])
-        self.assertIsNot(ds2, ds1)
-        self.assertIn('time', ds2)
-        self.assertIn('temperature_min', ds2)
-        self.assertIn('temperature_max', ds2)
-        self.assertIn('precipitation_min', ds2)
-        self.assertIn('precipitation_max', ds2)
-        self.assertEqual(('time',), ds2.time.dims)
-        self.assertEqual(('time', 'lat', 'lon'), ds2.temperature_min.dims)
-        self.assertEqual(('time', 'lat', 'lon'), ds2.temperature_max.dims)
-        self.assertEqual(('time', 'lat', 'lon'), ds2.precipitation_min.dims)
-        self.assertEqual(('time', 'lat', 'lon'), ds2.precipitation_max.dims)
-        self.assertEqual((num_times / 3,), ds2.time.shape)
-        self.assertEqual((num_times / 3, 180, 360), ds2.temperature_min.shape)
-        self.assertEqual((num_times / 3, 180, 360), ds2.temperature_max.shape)
-        self.assertEqual((num_times / 3, 180, 360), ds2.precipitation_min.shape)
-        self.assertEqual((num_times / 3, 180, 360), ds2.precipitation_max.shape)
-        np.testing.assert_equal(ds2.time.values,
-                                np.array(pd.to_datetime(
-                                    ['2017-07-01', '2017-07-04', '2017-07-07', '2017-07-10',
-                                     '2017-07-13', '2017-07-16', '2017-07-19', '2017-07-22',
-                                     '2017-07-25', '2017-07-28'])))
-        np.testing.assert_allclose(ds2.temperature_min.values[..., 0, 0],
-                                   np.array([272., 272.3, 272.6, 272.9, 273.2, 273.5, 273.8, 274.1, 274.4, 274.7]))
-        np.testing.assert_allclose(ds2.temperature_max.values[..., 0, 0],
-                                   np.array([272.2, 272.5, 272.8, 273.1, 273.4, 273.7, 274., 274.3, 274.6, 274.9]))
-        np.testing.assert_allclose(ds2.precipitation_min.values[..., 0, 0],
-                                   np.array([119.6, 119., 118.4, 117.8, 117.2, 116.6, 116., 115.4, 114.8, 114.2]))
-        np.testing.assert_allclose(ds2.precipitation_max.values[..., 0, 0],
-                                   np.array([120., 119.4, 118.8, 118.2, 117.6, 117., 116.4, 115.8, 115.2, 114.6]))
+        time = []
+        periods = ['1D', '1D', '3D', '4D', '2D']
+        t = pd.to_datetime('2017-07-01T10:30:15Z', utc=True)
+        for i in range(num_times):
+            time.append(t.isoformat())
+            t += pd.to_timedelta(periods[i % len(periods)])
+
+        temperature, precipitation = zip(*[(272 + 0.1 * i, 120 - 0.2 * i) for i in range(num_times)])
+
+        input_cube = new_test_dataset(time, temperature=temperature, precipitation=precipitation)
+        input_cube = chunk_dataset(input_cube, chunk_sizes=dict(time=1, lat=90, lon=180))
+        self.input_cube = input_cube
+
+    def test_resample_in_time(self):
+        resampled_cube = resample_in_time(self.input_cube, '2W', ['min', 'max'])
+        self.assertIsNot(resampled_cube, self.input_cube)
+        self.assertIn('time', resampled_cube)
+        self.assertIn('temperature_min', resampled_cube)
+        self.assertIn('temperature_max', resampled_cube)
+        self.assertIn('precipitation_min', resampled_cube)
+        self.assertIn('precipitation_max', resampled_cube)
+        self.assertEqual(('time',), resampled_cube.time.dims)
+        self.assertEqual(('time', 'lat', 'lon'), resampled_cube.temperature_min.dims)
+        self.assertEqual(('time', 'lat', 'lon'), resampled_cube.temperature_max.dims)
+        self.assertEqual(('time', 'lat', 'lon'), resampled_cube.precipitation_min.dims)
+        self.assertEqual(('time', 'lat', 'lon'), resampled_cube.precipitation_max.dims)
+        self.assertEqual((6,), resampled_cube.time.shape)
+        self.assertEqual((6, 180, 360), resampled_cube.temperature_min.shape)
+        self.assertEqual((6, 180, 360), resampled_cube.temperature_max.shape)
+        self.assertEqual((6, 180, 360), resampled_cube.precipitation_min.shape)
+        self.assertEqual((6, 180, 360), resampled_cube.precipitation_max.shape)
+        np.testing.assert_equal(resampled_cube.time.values,
+                                np.array(
+                                    ['2017-06-25T00:00:00Z', '2017-07-09T00:00:00Z',
+                                     '2017-07-23T00:00:00Z', '2017-08-06T00:00:00Z',
+                                     '2017-08-20T00:00:00Z', '2017-09-03T00:00:00Z'], dtype=np.datetime64))
+        np.testing.assert_allclose(resampled_cube.temperature_min.values[..., 0, 0],
+                                   np.array([272.0, 272.4, 273.0, 273.8, 274.4, 274.9]))
+        np.testing.assert_allclose(resampled_cube.temperature_max.values[..., 0, 0],
+                                   np.array([272.3, 272.9, 273.7, 274.3, 274.8, 274.9]))
+        np.testing.assert_allclose(resampled_cube.precipitation_min.values[..., 0, 0],
+                                   np.array([119.4, 118.2, 116.6, 115.4, 114.4, 114.2]))
+        np.testing.assert_allclose(resampled_cube.precipitation_max.values[..., 0, 0],
+                                   np.array([120.0, 119.2, 118.0, 116.4, 115.2, 114.2]))
+
+        schema = CubeSchema.new(resampled_cube)
+        self.assertEqual(3, schema.ndim)
+        self.assertEqual(('time', 'lat', 'lon'), schema.dims)
+        self.assertEqual((6, 180, 360), schema.shape)
+        self.assertEqual((1, 90, 180), schema.chunks)
+
+    def test_resample_in_time_with_time_chunk_size(self):
+        resampled_cube = resample_in_time(self.input_cube, '2D', ['min', 'max'], time_chunk_size=5)
+        schema = CubeSchema.new(resampled_cube)
+        self.assertEqual(3, schema.ndim)
+        self.assertEqual(('time', 'lat', 'lon'), schema.dims)
+        self.assertEqual((33, 180, 360), schema.shape)
+        self.assertEqual((5, 90, 180), schema.chunks)
diff --git a/test/sampledata.py b/test/sampledata.py
index db010201..8e6eea49 100644
--- a/test/sampledata.py
+++ b/test/sampledata.py
@@ -29,8 +29,10 @@ def new_test_dataset(time, height=180, **indexers):
             raise ValueError()
         data_vars[name] = (['time', 'lat', 'lon'],
                            np.concatenate(tuple(np.full(shape, values[i]) for i in range(num_times))))
+    time = np.array(pd.to_datetime(time), dtype=np.datetime64)
     return xr.Dataset(data_vars,
-                      coords=dict(time=(['time'], pd.to_datetime(time)),
+                      coords=dict(time=(
+                      ['time'], time, dict(units='nanoseconds since 1970-01-01', calendar='proleptic_gregorian')),
                                   lat=(['lat'], np.linspace(-90 + res, +90 - res, height)),
                                   lon=(['lon'], np.linspace(-180 + res, +180 - res, width))))