Repository: dcs4cop/xcube
Pull request #39 (2019-03-18 13:42:39)

Forman pyramidize
When you accept, 

* xcube will have a new CLI command `levels` that can write multi-level pyramids from NetCDF and ZARR inputs
* xcube will have a new API functions `compute_levels`, `write_levels`, `read_levels` from module `xcube.api,levels`
* Any `xarray.Dataset` will have an `xcube` property that has a new `levels()` method delegating to API function `compute_levels`



Patch:
diff --git a/README.md b/README.md
index 74590651..2044ec61 100644
--- a/README.md
+++ b/README.md
@@ -46,6 +46,7 @@ with [coverage report](https://pytest-cov.readthedocs.io/en/latest/reporting.htm
 
 
 # Tools
+
 ## `xcube chunk`
 
     $ xcube chunk --help
@@ -275,5 +276,36 @@ Example:
     $ xcube gen -a -s 2000,1000 -r 0,50,5,52.5 -v conc_chl,conc_tsm,kd489,c2rcc_flags,quality_flags -n hiroc-cube -t snap-c2rcc D:\OneDrive\BC\EOData\HIGHROC\2017\01\*.nc
 
 
+## `xcube level`
+
+    $ xcube level --help
 
+    Usage: xcube level [OPTIONS] <input>
+    
+      Transform the given dataset by <input> into the levels of a multi-level
+      pyramid with spatial resolution decreasing by a factor of two in both
+      spatial dimensions and write the result to directory <output>.
+    
+    Options:
+      -o, --output <output>           Output directory. If omitted,
+                                      "<input>.levels" will be used.
+      -l, --link                      Link the <input> instead of converting it to
+                                      a level zero dataset. Use with care, as the
+                                      <input>'s internal spatial chunk sizes may
+                                      be inappropriate for imaging purposes.
+      -t, --tile-size <tile-size>     Tile size, given as single integer number or
+                                      as <tile-width>,<tile-height>. If omitted,
+                                      the tile size will be derived from the
+                                      <input>'s internal spatial chunk sizes. If
+                                      the <input> is not chunked, tile size will
+                                      be 512.
+      -n, --num-levels-max <num-levels-max>
+                                      Maximum number of levels to generate. If not
+                                      given, the number of levels will be derived
+                                      from spatial dimension and tile sizes.
+      --help                          Show this message and exit.
+
+    
+Example:
 
+    $ xcube level -l -t 720 data/cubes/test-cube.zarr
diff --git a/test/genl3/__init__.py b/test/genl3/__init__.py
deleted file mode 100644
index e69de29b..00000000
diff --git a/xcube/api/api.py b/xcube/api/api.py
index 4f211c39..e426f731 100644
--- a/xcube/api/api.py
+++ b/xcube/api/api.py
@@ -4,13 +4,14 @@ from typing import Dict, List, Mapping, Any, Union, Sequence
 import numpy as np
 import pandas as pd
 import xarray as xr
+from xcube.api.levels import compute_levels
 
 from .vars_to_dim import vars_to_dim
 from .chunk import chunk_dataset
 from .dump import dump_dataset
-from .new import new_cube
 from .extract import get_cube_values_for_points, get_cube_point_indexes, get_cube_values_for_indexes, \
     get_dataset_indexes, INDEX_NAME_PATTERN
+from .new import new_cube
 from .readwrite import read_cube, open_cube, write_cube
 from .select import select_vars
 from .verify import verify_cube
@@ -234,17 +235,17 @@ class XCubeAPI:
         return vars_to_dim(self._dataset, dim_name)
 
     def dump(self,
-             variable_names=None,
+             var_names=None,
              show_var_encoding=False) -> str:
         """
         Dump this dataset or its variables into a text string.
 
-        :param variable_names: names of variables to be dumped
+        :param var_names: names of variables to be dumped
         :param show_var_encoding: also dump variable encodings?
         :return: the dataset dump
         """
         return dump_dataset(self._dataset,
-                            variable_names=variable_names,
+                            var_names=var_names,
                             show_var_encoding=show_var_encoding)
 
     def verify(self) -> List[str]:
@@ -265,3 +266,22 @@ class XCubeAPI:
         :return: A new dataset. It is empty, if *var_names* is empty. It is *dataset*, if *var_names* is None.
         """
         return select_vars(self._dataset, var_names)
+
+    def levels(self, **kwargs) -> List[xr.Dataset]:
+        """
+        Transform this dataset into the levels of a multi-level pyramid with spatial resolution
+        decreasing by a factor of two in both spatial dimensions.
+
+        It is assumed that the spatial dimensions of each variable are the inner-most, that is, the last two elements
+        of a variable's shape provide the spatial dimension sizes.
+
+        :param spatial_dims: If given, only variables are considered whose last to dimension elements match the given *spatial_dims*.
+        :param spatial_shape: If given, only variables are considered whose last to shape elements match the given *spatial_shape*.
+        :param spatial_tile_shape: If given, chunking will match the provided *spatial_tile_shape*.
+        :param var_names: Variables to consider. If None, all variables with at least two dimensions are considered.
+        :param max_num_levels: If given, the maximum number of pyramid levels.
+        :param post_process_level: If given, the function will be called for each level and must return a dataset.
+        :param progress_monitor: If given, the function will be called for each level.
+        :return: A list of dataset instances representing the multi-level pyramid.
+        """
+        return compute_levels(self._dataset, **kwargs)
diff --git a/xcube/api/dump.py b/xcube/api/dump.py
index cabf0679..bf779cb7 100644
--- a/xcube/api/dump.py
+++ b/xcube/api/dump.py
@@ -2,18 +2,18 @@ import xarray as xr
 
 
 def dump_dataset(dataset: xr.Dataset,
-                 variable_names=None,
+                 var_names=None,
                  show_var_encoding=False) -> str:
     """
     Dumps a dataset or variables into a text string.
 
     :param dataset: input dataset
-    :param variable_names: names of variables to be dumped
+    :param var_names: names of variables to be dumped
     :param show_var_encoding: also dump variable encodings?
     :return: the dataset dump
     """
     lines = []
-    if not variable_names:
+    if not var_names:
         lines.append(str(dataset))
         if show_var_encoding:
             for var_name, var in dataset.coords.items():
@@ -23,7 +23,7 @@ def dump_dataset(dataset: xr.Dataset,
                 if var.encoding:
                     lines.append(dump_var_encoding(var, header=f"Encoding for data variable {var_name!r}:"))
     else:
-        for var_name in variable_names:
+        for var_name in var_names:
             var = dataset[var_name]
             lines.append(str(var))
             if show_var_encoding and var.encoding:
diff --git a/xcube/api/levels.py b/xcube/api/levels.py
new file mode 100644
index 00000000..d19b8f1c
--- /dev/null
+++ b/xcube/api/levels.py
@@ -0,0 +1,281 @@
+import os
+from typing import List, Callable, Sequence, Optional, Tuple
+
+import xarray as xr
+
+from .readwrite import read_dataset
+
+PyramidLevelCallback = Callable[[xr.Dataset, int, int], Optional[xr.Dataset]]
+
+
+def compute_levels(dataset: xr.Dataset,
+                   spatial_dims: Tuple[str, str] = None,
+                   spatial_shape: Tuple[int, int] = None,
+                   spatial_tile_shape: Tuple[int, int] = None,
+                   var_names: Sequence[str] = None,
+                   num_levels_max: int = None,
+                   post_process_level: PyramidLevelCallback = None,
+                   progress_monitor: PyramidLevelCallback = None) -> List[xr.Dataset]:
+    """
+    Transform the given *dataset* into the levels of a multi-level pyramid with spatial resolution
+    decreasing by a factor of two in both spatial dimensions.
+
+    It is assumed that the spatial dimensions of each variable are the inner-most, that is, the last two elements
+    of a variable's shape provide the spatial dimension sizes.
+
+    :param dataset: The input dataset to be turned into a multi-level pyramid.
+    :param spatial_dims: If given, only variables are considered whose last to dimension elements match the given *spatial_dims*.
+    :param spatial_shape: If given, only variables are considered whose last to shape elements match the given *spatial_shape*.
+    :param spatial_tile_shape: If given, chunking will match the provided *spatial_tile_shape*.
+    :param var_names: Variables to consider. If None, all variables with at least two dimensions are considered.
+    :param num_levels_max: If given, the maximum number of pyramid levels.
+    :param post_process_level: If given, the function will be called for each level and must return a dataset.
+    :param progress_monitor: If given, the function will be called for each level.
+    :return: A list of dataset instances representing the multi-level pyramid.
+    """
+    dropped_vars, spatial_shape, spatial_tile_shape = _filter_level_source_dataset(dataset, var_names, spatial_dims,
+                                                                                   spatial_shape, spatial_tile_shape)
+    if dropped_vars:
+        dataset = dataset.drop(dropped_vars)
+
+    if not tuple(dataset.data_vars):
+        raise ValueError("cannot compute pyramid levels because no suitable data variables were found")
+
+    if spatial_tile_shape is None:
+        spatial_tile_shape = min(spatial_shape[0], 512), min(spatial_shape[1], 512)
+
+    # Count num_levels
+    level_shapes = _compute_level_shapes(spatial_shape, spatial_tile_shape, num_levels_max=num_levels_max)
+    num_levels = len(level_shapes)
+
+    # Compute levels
+    level_dataset = dataset
+    level_datasets = []
+    for level in range(num_levels):
+        if level > 0:
+            # Down-sample levels
+            downsampled_vars = {}
+            for var_name in level_dataset.data_vars:
+                var = level_dataset.data_vars[var_name]
+                # For time being, we use the simplest and likely fastest downsampling I can think of
+                downsampled_var = var[..., ::2, ::2]
+                if downsampled_var.shape[-2:] != level_shapes[level]:
+                    import warnings
+                    warnings.warn(f"unexpected spatial shape for down-sampled variable {var_name!r}:"
+                                  f" expected {level_shapes[level]}, but found {downsampled_var.shape[-2:]}")
+                downsampled_vars[var_name] = downsampled_var
+            level_dataset = xr.Dataset(downsampled_vars, attrs=level_dataset.attrs)
+
+        level_dataset = _tile_level_dataset(level_dataset, spatial_tile_shape)
+
+        # Apply post processor, if any
+        if post_process_level is not None:
+            level_dataset = post_process_level(level_dataset, len(level_datasets), num_levels)
+
+        # Inform progress monitor, if any
+        if progress_monitor is not None:
+            progress_monitor(level_dataset, len(level_datasets), num_levels)
+
+        # Collect level dataset
+        level_datasets.append(level_dataset)
+
+    return level_datasets
+
+
+def write_levels(output_path: str,
+                 dataset: xr.Dataset = None,
+                 input_path: str = None,
+                 link_input: bool = False,
+                 progress_monitor: PyramidLevelCallback = None,
+                 **kwargs) -> List[xr.Dataset]:
+    """
+    Transform the given dataset given by a *dataset* instance or *input_path* string
+    into the levels of a multi-level pyramid with spatial resolution
+    decreasing by a factor of two in both spatial dimensions and write them to *output_path*.
+
+    One of *dataset* and *input_path* must be given.
+
+    :param output_path: Output path
+    :param dataset: Dataset to be converted and written as levels.
+    :param input_path: Input path to a dataset to be transformed and written as levels.
+    :param link_input: Just link the dataset at level zero instead of writing it.
+    :param progress_monitor: An optional progress monitor.
+    :param kwargs: Keyword-arguments accepted by the ``compute_levels()`` function.
+    :return: A list of dataset instances representing the multi-level pyramid.
+    """
+    if dataset is None and input_path is None:
+        raise ValueError("at least one of dataset or input_path must be given")
+
+    if link_input and input_path is None:
+        raise ValueError("input_path must be provided to link input")
+
+    _post_process_level = kwargs.pop("post_process_level") if "post_process_level" in kwargs else None
+
+    def post_process_level(level_dataset, index, num_levels):
+        if _post_process_level is not None:
+            level_dataset = _post_process_level(level_dataset, index, num_levels)
+
+        if index == 0 and link_input:
+            with open(os.path.join(output_path, f"{index}.link"), "w") as fp:
+                fp.write(input_path)
+        else:
+            path = os.path.join(output_path, f"{index}.zarr")
+            level_dataset.to_zarr(path)
+            level_dataset.close()
+            level_dataset = xr.open_zarr(path)
+
+        if progress_monitor is not None:
+            progress_monitor(level_dataset, index, num_levels)
+
+        return level_dataset
+
+    if not os.path.exists(output_path):
+        os.makedirs(output_path)
+
+    if dataset is None:
+        dataset = read_dataset(input_path)
+
+    return compute_levels(dataset, post_process_level=post_process_level, **kwargs)
+
+
+def read_levels(dir_path: str,
+                progress_monitor: PyramidLevelCallback = None) -> List[xr.Dataset]:
+    """
+    Read the of a multi-level pyramid with spatial resolution
+    decreasing by a factor of two in both spatial dimensions.
+
+    :param dir_path: The directory path.
+    :param progress_monitor: An optional progress monitor.
+    :return: A list of dataset instances representing the multi-level pyramid.
+    """
+    file_paths = os.listdir(dir_path)
+    level_paths = {}
+    num_levels = -1
+    for filename in file_paths:
+        file_path = os.path.join(dir_path, filename)
+        basename, ext = os.path.splitext(filename)
+        if basename.isdigit():
+            index = int(basename)
+            num_levels = max(num_levels, index + 1)
+            if os.path.isfile(file_path) and ext == ".link":
+                level_paths[index] = (ext, file_path)
+            elif os.path.isdir(file_path) and ext == ".zarr":
+                level_paths[index] = (ext, file_path)
+
+    if num_levels != len(level_paths):
+        raise ValueError(f"Inconsistent pyramid directory:"
+                         f" expected {num_levels} but found {len(level_paths)} entries:"
+                         f" {dir_path}")
+
+    levels = []
+    for index in range(num_levels):
+        ext, file_path = level_paths[index]
+        if ext == ".link":
+            with open(file_path, "r") as fp:
+                link_file_path = fp.read()
+            dataset = xr.open_zarr(link_file_path)
+        else:
+            dataset = xr.open_zarr(file_path)
+        if progress_monitor is not None:
+            progress_monitor(dataset, index, num_levels)
+        levels.append(dataset)
+    return levels
+
+
+def _tile_chunk(size, tile_size):
+    last_tile_size = size % tile_size
+    if last_tile_size != 0:
+        return (tile_size,) * (size // tile_size) + (last_tile_size,)
+    return tile_size
+
+
+def _tile_level_dataset(level_dataset, spatial_tile_shape):
+    tile_height, tile_width = spatial_tile_shape
+
+    # Chunk variables in level dataset according to spatial_tile_shape
+    chunked_vars = {}
+
+    # Chunk data variables according to tile size
+    for var_name in level_dataset.data_vars:
+        var = level_dataset.data_vars[var_name]
+        height, width = var.shape[-2:]
+        zarr_chunks = (1,) * (var.ndim - 2) + (tile_height,) + (tile_width,)
+        dask_chunks = (1,) * (var.ndim - 2) + (_tile_chunk(height, tile_height),) + (_tile_chunk(width, tile_width),)
+        dask_chunks = {var.dims[i]: dask_chunks[i] for i in range(var.ndim)}
+        chunked_var = var.chunk(chunks=dask_chunks)
+        chunked_var.encoding.update(chunks=zarr_chunks)
+        chunked_vars[var_name] = chunked_var
+
+    # Make coordinate variable chunks equal to their shape
+    # TODO (forman): find out if chunking the spatial coordinates according to tile size improves performance
+    for var_name in level_dataset.coords:
+        var = level_dataset.coords[var_name]
+        zarr_chunks = var.shape
+        dask_chunks = {var.dims[i]: var.shape[i] for i in range(var.ndim)}
+        chunked_var = var.chunk(chunks=dask_chunks)
+        chunked_var.encoding.update(chunks=zarr_chunks)
+        chunked_vars[var_name] = chunked_var
+
+    return level_dataset.assign(variables=chunked_vars)
+
+
+def _compute_level_shapes(spatial_shape, spatial_tile_shape, num_levels_max=None) -> List[Tuple[int, int]]:
+    height, width = spatial_shape
+    tile_height, tile_width = spatial_tile_shape
+    num_levels_max = num_levels_max or -1
+    level_shapes = [(height, width)]
+    while True:
+        width = (width + 1) // 2
+        height = (height + 1) // 2
+        if width < tile_width or height < tile_height or num_levels_max == len(level_shapes):
+            break
+        level_shapes.append((height, width))
+    return level_shapes
+
+
+def _filter_level_source_dataset(dataset,
+                                 var_names=None,
+                                 spatial_dims=None,
+                                 spatial_shape=None,
+                                 spatial_tile_shape=None):
+    if var_names:
+        var_names = set(var_names)
+        dropped_vars = list(set(dataset.data_vars).difference(var_names))
+    else:
+        var_names = set(dataset.data_vars)
+        dropped_vars = []
+
+    # Collect data variables to be dropped, derive missing information from spatial data variables
+    for var_name in var_names:
+
+        if var_name not in dataset.data_vars:
+            raise ValueError(f"variable {var_name} not found")
+
+        var = dataset[var_name]
+
+        if var.ndim < 2:
+            # Must have at least the two spatial dimensions
+            dropped_vars.append(var_name)
+            continue
+
+        if spatial_dims is None:
+            spatial_dims = var.dims[-2:]
+        elif spatial_dims != var.dims[-2:]:
+            # Spatial dimensions don't fit
+            dropped_vars.append(var_name)
+            continue
+
+        if spatial_shape is None:
+            spatial_shape = var.shape[-2:]
+        elif spatial_shape != var.shape[-2:]:
+            # Spatial dimension sizes don't fit
+            dropped_vars.append(var_name)
+            continue
+
+        if spatial_tile_shape is None and var.chunks is not None:
+            def chunk_to_int(chunk):
+                return chunk if isinstance(chunk, int) else max(chunk)
+
+            spatial_tile_shape = tuple(map(chunk_to_int, var.chunks[-2:]))
+
+    return dropped_vars, spatial_shape, spatial_tile_shape
diff --git a/xcube/api/vars_to_dim.py b/xcube/api/vars_to_dim.py
index 57f293b3..753d6a81 100644
--- a/xcube/api/vars_to_dim.py
+++ b/xcube/api/vars_to_dim.py
@@ -21,20 +21,36 @@
 
 import xarray as xr
 
-def vars_to_dim(cube: xr.Dataset, dim_name: str = 'newdim'):
+
+def vars_to_dim(cube: xr.Dataset,
+                dim_name: str = 'var',
+                var_name='data',
+                cube_asserted: bool = False):
     """
-    Convert data variables into a dimension
+    Convert data variables into a dimension.
 
     :param cube: The data cube.
-    :param dim_name: The name of the new dimension ['vars']
-    :return: A new data cube with the new dimension.
+    :param dim_name: The name of the new dimension and coordinate variable. Defaults to 'var'.
+    :param var_name: The name of the new, single data variable. Defaults to 'data'.
+    :param cube_asserted: If False, *cube* will be verified, otherwise it is expected to be a valid cube.
+    :return: A new data cube with data variables turned into a new dimension.
     """
 
-    if len(cube.variables) == 0:
-        raise ValueError("ERROR: Cube has no variables. Exiting.")
+    if not cube_asserted:
+        assert_cube(cube)
+
+    if var_name == dim_name:
+        raise ValueError("var_name must be different from dim_name")
+
+    data_var_names = [data_var_name for data_var_name in cube.data_vars]
+    if not data_var_names:
+        raise ValueError("cube must not be empty")
+
+    da = xr.concat([cube[data_var_name] for data_var_name in data_var_names], dim_name)
+    new_coord_var = xr.DataArray(data_var_names, dims=[dim_name])
+    da = da.assign_coords(**{dim_name: new_coord_var})
+
+    return xr.Dataset(dict(**{var_name: da}))
 
-    da = xr.concat([cube[var_name] for var_name in cube.data_vars], dim_name)
-    var_coords = xr.DataArray([var_name for var_name in cube.data_vars], dims=[dim_name])
-    da = da.assign_coords(**{dim_name: var_coords})
 
-    return xr.Dataset(dict(**{dim_name + '_vars': da}))
+from .verify import assert_cube
diff --git a/xcube/cli/cli.py b/xcube/cli/cli.py
index 7fd6df3f..7a597ee7 100644
--- a/xcube/cli/cli.py
+++ b/xcube/cli/cli.py
@@ -6,6 +6,8 @@ from xcube.cli.gen import gen
 from xcube.cli.grid import grid
 from xcube.version import version
 
+DEFAULT_TILE_SIZE = 512
+
 
 def _parse_kwargs(value: str, metavar: str = None) -> Dict[str, Any]:
     if value:
@@ -102,6 +104,76 @@ def chunk(input, output, format=None, params=None, chunks=None):
         write_dataset(chunked_dataset, output_path=output, format_name=format_name, **write_kwargs)
 
 
+@click.command(name="level")
+@click.argument('input', metavar='<input>')
+@click.option('--output', '-o', metavar='<output>',
+              help='Output directory. If omitted, "<input>.levels" will be used.')
+@click.option('--link', '-l', is_flag=True, flag_value=True,
+              help='Link the <input> instead of converting it to a level zero dataset. '
+                   'Use with care, as the <input>\'s internal spatial chunk sizes may be inappropriate '
+                   'for imaging purposes.')
+@click.option('--tile-size', '-t', metavar='<tile-size>',
+              help=f'Tile size, given as single integer number or as <tile-width>,<tile-height>. '
+              f'If omitted, the tile size will be derived from the <input>\'s '
+              f'internal spatial chunk sizes. '
+              f'If the <input> is not chunked, tile size will be {DEFAULT_TILE_SIZE}.')
+@click.option('--num-levels-max', '-n', metavar='<num-levels-max>', type=int,
+              help=f'Maximum number of levels to generate. '
+              f'If not given, the number of levels will be derived from '
+              f'spatial dimension and tile sizes.')
+def level(input, output, link, tile_size, num_levels_max):
+    """
+    Transform the given dataset by <input> into the levels of a multi-level pyramid with spatial resolution
+    decreasing by a factor of two in both spatial dimensions and write the result to directory <output>.
+    """
+    import time
+    import os
+    from xcube.api.levels import write_levels
+
+    input_path = input
+    output_path = output
+    link_input = link
+
+    if num_levels_max is not None and num_levels_max < 1:
+        raise click.ClickException(f"<num-levels-max> must be a positive integer")
+
+    if not output_path:
+        basename, ext = os.path.splitext(input_path)
+        output_path = os.path.join(os.path.dirname(input_path), basename + ".levels")
+
+    if os.path.exists(output_path):
+        raise click.ClickException(f"output {output_path!r} already exists")
+
+    spatial_tile_shape = None
+    if tile_size is not None:
+        try:
+            tile_size = int(tile_size)
+            tile_size = tile_size, tile_size
+        except ValueError:
+            tile_size = map(int, tile_size.split(","))
+            if tile_size != 2:
+                raise click.ClickException("Expected a pair of positive integers <tile-width>,<tile-height>")
+        if tile_size[0] < 1 or tile_size[1] < 1:
+            raise click.ClickException("<tile-size> must comprise positive integers")
+        spatial_tile_shape = tile_size[1], tile_size[0]
+
+    start_time = t0 = time.perf_counter()
+
+    # noinspection PyUnusedLocal
+    def progress_monitor(dataset, index, num_levels):
+        nonlocal t0
+        print(f"level {index + 1} of {num_levels} written after {time.perf_counter() - t0} seconds")
+        t0 = time.perf_counter()
+
+    levels = write_levels(output_path,
+                          input_path=input_path,
+                          link_input=link_input,
+                          progress_monitor=progress_monitor,
+                          spatial_tile_shape=spatial_tile_shape,
+                          num_levels_max=num_levels_max)
+    print(f"{len(levels)} level(s) written into {output_path} after {time.perf_counter() - start_time} seconds")
+
+
 @click.command(name="dump")
 @click.argument('dataset', metavar='<path>')
 @click.option('--variable', '-v', metavar='<variable>', multiple=True,
@@ -114,43 +186,42 @@ def dump(dataset, variable, encoding):
     """
     from xcube.api import open_dataset, dump_dataset
     with open_dataset(dataset) as ds:
-        text = dump_dataset(ds, variable_names=variable, show_var_encoding=encoding)
+        text = dump_dataset(ds, var_names=variable, show_var_encoding=encoding)
         print(text)
 
 
 # noinspection PyShadowingBuiltins
 @click.command(name="vars2dim")
 @click.argument('cube', metavar='<cube>')
-@click.option('--dim_name', '-d', metavar='<dimension name>', help="Name of the new dimension ['var']")
+@click.option('--dim_name', '-d', metavar='<dim-name>',
+              default='var',
+              help='Name of the new dimension into variables. Defaults to "var".')
+@click.option('--var_name', '-v', metavar='<var-name>',
+              default='data',
+              help='Name of the new variable that includes all variables. Defaults to "data".')
 @click.option('--output', '-o', metavar='<output>',
               help="Output file.")
 @click.option('--format', '-f', metavar='<format>', type=click.Choice(['zarr', 'netcdf']),
               help="Format of the output. If not given, guessed from <output>.")
-def vars2dim(cube, dim_name, output=None, format=None):
+def vars2dim(cube, var_name, dim_name, output=None, format=None):
     """
     Write a new dataset with identical data and compression but with new chunk sizes.
     """
 
     from xcube.util.dsio import guess_dataset_format
+    from xcube.api import open_dataset, vars_to_dim, write_dataset
     import os
 
-    if not dim_name:
-        dim_name = 'newdim'
-
     if not output:
-        basename = os.path.splitext(cube)[0]
-        ext = os.path.splitext(cube)[1]
-
-        output = basename + '_converted' + ext
+        dirname = os.path.dirname(cube)
+        basename = os.path.basename(cube)
+        basename, ext = os.path.splitext(basename)
+        output = os.path.join(dirname, basename + '-vars2dim' + ext)
 
     format_name = format if format else guess_dataset_format(output)
 
-    from xcube.api import open_dataset, vars_to_dim, write_dataset
-
     with open_dataset(input_path=cube) as ds:
-
-        converted_dataset = vars_to_dim(ds, dim_name=dim_name)
-
+        converted_dataset = vars_to_dim(ds, dim_name=dim_name, var_name=var_name)
         write_dataset(converted_dataset, output_path=output, format_name=format_name)
 
 
@@ -169,6 +240,7 @@ cli.add_command(extract)
 cli.add_command(grid)
 cli.add_command(vars2dim)
 cli.add_command(gen)
+cli.add_command(level)
 
 
 def main(args=None):
diff --git a/xcube/util/dsio.py b/xcube/util/dsio.py
index fb9bde65..f5c020be 100644
--- a/xcube/util/dsio.py
+++ b/xcube/util/dsio.py
@@ -21,6 +21,8 @@
 
 import glob
 import os
+import shutil
+import warnings
 from abc import abstractmethod, ABCMeta
 from typing import Set, Callable, List, Optional, Dict, Iterable, Tuple
 
@@ -28,7 +30,6 @@ import pandas as pd
 import s3fs
 import xarray as xr
 import zarr
-
 from xcube.util.objreg import get_obj_registry
 
 FORMAT_NAME_EXCEL = "excel"
@@ -485,12 +486,14 @@ def rimraf(path):
 
     :param path:  directory or single file
     """
-    import os
     if os.path.isdir(path):
-        import shutil
-        shutil.rmtree(path, ignore_errors=True)
+        try:
+            shutil.rmtree(path, ignore_errors=False)
+        except OSError:
+            warnings.warn(f"failed to remove file {path}")
     elif os.path.isfile(path):
         try:
             os.remove(path)
         except OSError:
+            warnings.warn(f"failed to remove file {path}")
             pass


Test patch:
diff --git a/test/api/test_api.py b/test/api/test_api.py
index 79608880..d1cd0cc7 100644
--- a/test/api/test_api.py
+++ b/test/api/test_api.py
@@ -30,3 +30,19 @@ class XCubeAPITest(unittest.TestCase):
                                    precipitation=0.4, temperature=275.2)
 
         self.assertIsInstance(dataset.xcube.vars_to_dim(), xr.Dataset)
+
+    def test_levels(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)
+        levels = dataset.xcube.levels(spatial_tile_shape=(45, 45))
+        self.assertIsInstance(levels, list)
+        self.assertEqual(3, len(levels))
+        self.assertTrue(all(isinstance(level, xr.Dataset) for level in levels))
+        self.assertTrue(all("precipitation" in level for level in levels))
+        self.assertTrue(all("temperature" in level for level in levels))
+        self.assertEqual([(5, 180, 360), (5, 90, 180), (5, 45, 90)],
+                         [level.precipitation.shape for level in levels])
+        self.assertEqual([((1, 1, 1, 1, 1), (45, 45, 45, 45), (45, 45, 45, 45, 45, 45, 45, 45)),
+                          ((1, 1, 1, 1, 1), (45, 45), (45, 45, 45, 45)),
+                          ((1, 1, 1, 1, 1), (45,), (45, 45))],
+                         [level.precipitation.chunks for level in levels])
diff --git a/test/api/test_levels.py b/test/api/test_levels.py
new file mode 100644
index 00000000..585b7053
--- /dev/null
+++ b/test/api/test_levels.py
@@ -0,0 +1,154 @@
+import os
+import time
+import unittest
+
+import numpy as np
+import xarray as xr
+
+from xcube.api.levels import compute_levels, write_levels, read_levels
+from xcube.util.dsio import rimraf
+
+
+def get_path(filename):
+    return os.path.join(os.path.dirname(__file__), filename)
+
+
+class PyramidTest(unittest.TestCase):
+    @classmethod
+    def create_test_dataset(cls, shape, chunks=None):
+        size = int(np.prod(shape))
+        dims = ["time", "y", "x"]
+        a_data = np.linspace(0, 1, size, dtype=np.float64).reshape(shape)
+        a = xr.DataArray(a_data, dims=dims)
+        b_data = np.linspace(-1, 0, size, dtype=np.float64).reshape(shape)
+        b = xr.DataArray(b_data, dims=dims)
+        if chunks:
+            a.encoding.update(chunks=chunks, chunksizes=chunks)
+            b.encoding.update(chunks=chunks, chunksizes=chunks)
+            return xr.Dataset(dict(a=a, b=b)).chunk(chunks={dims[i]: chunks[i] for i in range(len(dims))})
+        else:
+            return xr.Dataset(dict(a=a, b=b))
+
+    def test_compute_levels(self):
+        dataset = self.create_test_dataset(shape=(5, 200, 400), chunks=(1, 25, 25))
+        levels = compute_levels(dataset)
+        self._assert_levels_ok(levels,
+                               expected_num_levels=4,
+                               expected_shapes=[
+                                   (5, 200, 400),
+                                   (5, 100, 200),
+                                   (5, 50, 100),
+                                   (5, 25, 50),
+                               ],
+                               expected_chunks=[
+                                   ((1,) * 5, (25,) * 8, (25,) * 16),
+                                   ((1,) * 5, (25,) * 4, (25,) * 8),
+                                   ((1,) * 5, (25,) * 2, (25,) * 4),
+                                   ((1,) * 5, (25,) * 1, (25,) * 2),
+                               ])
+
+    def test_compute_levels_with_max_levels(self):
+        dataset = self.create_test_dataset(shape=(5, 200, 400), chunks=(1, 25, 25))
+        levels = compute_levels(dataset, num_levels_max=3)
+        self._assert_levels_ok(levels,
+                               expected_num_levels=3,
+                               expected_shapes=[
+                                   (5, 200, 400),
+                                   (5, 100, 200),
+                                   (5, 50, 100),
+                               ],
+                               expected_chunks=[
+                                   ((1,) * 5, (25,) * 8, (25,) * 16),
+                                   ((1,) * 5, (25,) * 4, (25,) * 8),
+                                   ((1,) * 5, (25,) * 2, (25,) * 4),
+                               ])
+
+    def test_write_read_levels_with_even_sizes(self):
+        shape = (5, 200, 400)
+        tile_shape = (25, 25)
+        expected_num_levels = 4
+        expected_shapes = [
+            (5, 200, 400),
+            (5, 100, 200),
+            (5, 50, 100),
+            (5, 25, 50),
+        ]
+        expected_chunks = [
+            ((1,) * 5, (25,) * 8, (25,) * 16),
+            ((1,) * 5, (25,) * 4, (25,) * 8),
+            ((1,) * 5, (25,) * 2, (25,) * 4),
+            ((1,) * 5, (25,) * 1, (25,) * 2),
+        ]
+        self._assert_io_ok(shape, tile_shape, expected_num_levels, expected_shapes, expected_chunks)
+
+    def test_write_read_levels_with_odd_sizes(self):
+        shape = (5, 203, 405)
+        tile_shape = (37, 38)
+        expected_num_levels = 3
+        expected_shapes = [
+            (5, 203, 405),
+            (5, 102, 203),
+            (5, 51, 102),
+        ]
+        expected_chunks = [
+            ((1,) * 5, (37, 37, 37, 37, 37, 18), (38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 25)),
+            ((1,) * 5, (37, 37, 28), (38, 38, 38, 38, 38, 13)),
+            ((1,) * 5, (37, 14), (38, 38, 26)),
+        ]
+        self._assert_io_ok(shape, tile_shape, expected_num_levels, expected_shapes, expected_chunks)
+
+    def _assert_io_ok(self, shape, tile_shape, expected_num_levels, expected_shapes, expected_chunks):
+
+        input_path = get_path("pyramid-input.nc")
+        output_path = get_path("pyramid-output")
+
+        rimraf(input_path)
+        rimraf(output_path)
+
+        try:
+            dataset = self.create_test_dataset(shape, chunks=None)
+            dataset.to_netcdf(input_path)
+
+            t0 = time.perf_counter()
+
+            levels = write_levels(output_path,
+                                  dataset=dataset,
+                                  spatial_tile_shape=tile_shape,
+                                  input_path=input_path)
+
+            print(f"write time total: ", time.perf_counter() - t0)
+
+            self._assert_levels_ok(levels, expected_num_levels, expected_shapes, expected_chunks)
+
+            t0 = time.perf_counter()
+
+            levels = read_levels(output_path)
+
+            print(f"read time total: ", time.perf_counter() - t0)
+
+            self._assert_levels_ok(levels, expected_num_levels, expected_shapes, expected_chunks)
+
+        finally:
+            rimraf(input_path)
+            rimraf(output_path)
+
+    def _assert_levels_ok(self, levels, expected_num_levels, expected_shapes, expected_chunks):
+        self.assertIsInstance(levels, list)
+        self.assertEqual(expected_num_levels, len(levels))
+
+        for i in range(expected_num_levels):
+            msg = f"at index {i}"
+            self.assertIsInstance(levels[i], xr.Dataset, msg=msg)
+            self.assertEqual(2, len(levels[i].data_vars), msg=msg)
+            self.assertIn("a", levels[i], msg=msg)
+            self.assertIn("b", levels[i], msg=msg)
+
+        for i in range(expected_num_levels):
+            msg = f"at index {i}"
+            self.assertEqual(expected_shapes[i], levels[i]["a"].shape, msg=msg)
+            self.assertEqual(expected_shapes[i], levels[i]["b"].shape, msg=msg)
+
+        for i in range(expected_num_levels):
+            msg = f"at index {i}"
+            self.assertEqual(expected_chunks[i], levels[i]["a"].chunks, msg=msg)
+            self.assertEqual(expected_chunks[i], levels[i]["b"].chunks, msg=msg)
diff --git a/test/api/test_vars_to_dim.py b/test/api/test_vars_to_dim.py
index 5966beb6..c9b92a22 100644
--- a/test/api/test_vars_to_dim.py
+++ b/test/api/test_vars_to_dim.py
@@ -1,22 +1,22 @@
 import unittest
-import xarray as xr
 
-from test.sampledata import new_test_dataset
+from xcube.api import new_cube
 from xcube.api.vars_to_dim import vars_to_dim
 
 
 class VarsToDimTest(unittest.TestCase):
     def test_vars_to_dim(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)
-
-        print(dataset.dims)
+        dataset = new_cube(variables=dict(precipitation=0.4, temperature=275.2))
 
         ds = vars_to_dim(dataset)
 
-        self.assertIn("newdim_vars", ds.variables)
-        newdim_vars = ds["newdim_vars"]
-        self.assertTrue(hasattr(newdim_vars, "encoding"))
-        self.assertEquals(len(dataset.dims), 3)
-
-        self.assertRaises(ValueError, vars_to_dim, xr.Dataset())
+        self.assertIn("var", ds.dims)
+        self.assertEqual(2, ds.dims["var"])
+        self.assertIn("var", ds.coords)
+        self.assertIn("data", ds.data_vars)
+        var_names = ds["var"]
+        self.assertEqual(("var",), var_names.dims)
+        self.assertTrue(hasattr(var_names, "encoding"))
+        self.assertEqual(2, len(var_names))
+        self.assertIn("precipitation", str(var_names[0]))
+        self.assertIn("temperature", str(var_names[1]))
diff --git a/test/cli/test_cli.py b/test/cli/test_cli.py
index 6f890c3d..a5ed343b 100644
--- a/test/cli/test_cli.py
+++ b/test/cli/test_cli.py
@@ -1,5 +1,4 @@
 import os
-import shutil
 import unittest
 from abc import ABCMeta
 from typing import List
@@ -7,10 +6,9 @@ from typing import List
 import click
 import click.testing
 import xarray as xr
-
-from test.sampledata import new_test_dataset
 from xcube.api.new import new_cube
 from xcube.cli.cli import cli, _parse_kwargs
+from xcube.util.dsio import rimraf
 
 TEST_NC_FILE = "test.nc"
 TEST_ZARR_DIR = "test.zarr"
@@ -22,17 +20,21 @@ class CliTest(unittest.TestCase, metaclass=ABCMeta):
         self.runner = click.testing.CliRunner()
         return self.runner.invoke(cli, args, catch_exceptions=False)
 
+    def outputs(self) -> List[str]:
+        return []
+
     def setUp(self):
-        super().setUp()
+        self._rm_outputs()
         dataset = new_cube(variables=dict(precipitation=0.4, temperature=275.2))
         dataset.to_netcdf(TEST_NC_FILE, mode="w")
         dataset.to_zarr(TEST_ZARR_DIR, mode="w")
 
     def tearDown(self):
-        if os.path.isdir(TEST_ZARR_DIR):
-            shutil.rmtree(TEST_ZARR_DIR, ignore_errors=True)
-        os.remove(TEST_NC_FILE)
-        super().tearDown()
+        self._rm_outputs()
+
+    def _rm_outputs(self):
+        for path in [TEST_ZARR_DIR, TEST_NC_FILE] + self.outputs():
+            rimraf(path)
 
 
 class DumpTest(CliTest):
@@ -69,32 +71,40 @@ class DumpTest(CliTest):
 
 
 class Vars2DimTest(CliTest):
+    TEST_OUTPUT = "test-vars2dim.zarr"
+
+    def outputs(self):
+        return [self.TEST_OUTPUT]
 
     def test_vars2dim(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)
-        dataset.to_zarr('tt.zarr')
-        output_path = "tt_converted.zarr"
-        result = self.invoke_cli(["vars2dim",
-                                  'tt.zarr'])
+        result = self.invoke_cli(["vars2dim", TEST_ZARR_DIR])
 
+        output_path = self.TEST_OUTPUT
         self.assertEqual(0, result.exit_code)
         self.assertTrue(os.path.isdir(output_path))
-        try:
-            ds = xr.open_zarr(output_path)
-            self.assertIn("newdim_vars", ds.variables)
-            newdim_vars = ds["newdim_vars"]
-            self.assertTrue(hasattr(newdim_vars, "encoding"))
-        finally:
-            print('hello')
-            shutil.rmtree(output_path, ignore_errors=True)
-            shutil.rmtree('tt.zarr', ignore_errors=True)
+
+        ds = xr.open_zarr(output_path)
+        self.assertIn("var", ds.dims)
+        self.assertEqual(2, ds.dims["var"])
+        self.assertIn("var", ds.coords)
+        self.assertIn("data", ds.data_vars)
+        var_names = ds["var"]
+        self.assertEqual(("var",), var_names.dims)
+        self.assertTrue(hasattr(var_names, "encoding"))
+        self.assertEqual("object", str(var_names.dtype))
+        self.assertEqual(2, len(var_names))
+        self.assertIn("precipitation", str(var_names[0]))
+        self.assertIn("temperature", str(var_names[1]))
 
 
 class ChunkTest(CliTest):
+    TEST_OUTPUT = "test-chunked.zarr"
+
+    def outputs(self):
+        return [self.TEST_OUTPUT]
 
     def test_chunk_zarr(self):
-        output_path = "test-chunked.zarr"
+        output_path = ChunkTest.TEST_OUTPUT
         result = self.invoke_cli(["chunk",
                                   TEST_ZARR_DIR,
                                   output_path,
@@ -102,15 +112,13 @@ class ChunkTest(CliTest):
         self.assertEqual("", result.output)
         self.assertEqual(0, result.exit_code)
         self.assertTrue(os.path.isdir(output_path))
-        try:
-            ds = xr.open_zarr(output_path)
-            self.assertIn("precipitation", ds)
-            precipitation = ds["precipitation"]
-            self.assertTrue(hasattr(precipitation, "encoding"))
-            self.assertIn("chunks", precipitation.encoding)
-            self.assertEqual(precipitation.encoding["chunks"], (1, 20, 40))
-        finally:
-            shutil.rmtree(output_path, ignore_errors=True)
+
+        ds = xr.open_zarr(output_path)
+        self.assertIn("precipitation", ds)
+        precipitation = ds["precipitation"]
+        self.assertTrue(hasattr(precipitation, "encoding"))
+        self.assertIn("chunks", precipitation.encoding)
+        self.assertEqual(precipitation.encoding["chunks"], (1, 20, 40))
 
     # TODO (forman): this test fails
     # netCDF4\_netCDF4.pyx:2437: in netCDF4._netCDF4.Dataset.createVariable
@@ -168,6 +176,90 @@ class ChunkTest(CliTest):
         self.assertEqual(1, result.exit_code)
 
 
+class LevelTest(CliTest):
+    TEST_OUTPUT = "test.levels"
+
+    def outputs(self):
+        return [LevelTest.TEST_OUTPUT]
+
+    def test_level_with_nc(self):
+        result = self.invoke_cli(["level",
+                                  "-t", "45",
+                                  "-o", LevelTest.TEST_OUTPUT,
+                                  TEST_NC_FILE,
+                                  ])
+        self.assertEqual(0, result.exit_code)
+        self.assertTrue(os.path.isdir(LevelTest.TEST_OUTPUT))
+        self.assertEqual({'0.zarr', '1.zarr', '2.zarr'}, set(os.listdir(LevelTest.TEST_OUTPUT)))
+
+    def test_level_with_zarr(self):
+        result = self.invoke_cli(["level",
+                                  "-t", "45",
+                                  "-o", LevelTest.TEST_OUTPUT,
+                                  TEST_ZARR_DIR,
+                                  ])
+        self.assertEqual(0, result.exit_code)
+        self.assertTrue(os.path.isdir(LevelTest.TEST_OUTPUT))
+        self.assertEqual({'0.zarr', '1.zarr', '2.zarr'}, set(os.listdir(LevelTest.TEST_OUTPUT)))
+
+    def test_level_with_zarr_link(self):
+        result = self.invoke_cli(["level",
+                                  "-l",
+                                  "-t", "45",
+                                  "-o", LevelTest.TEST_OUTPUT,
+                                  TEST_ZARR_DIR,
+                                  ])
+        self.assertEqual(0, result.exit_code)
+        self.assertTrue(os.path.isdir(LevelTest.TEST_OUTPUT))
+        self.assertEqual({'0.link', '1.zarr', '2.zarr'}, set(os.listdir(LevelTest.TEST_OUTPUT)))
+
+    def test_level_with_zarr_num_levels_max(self):
+        result = self.invoke_cli(["level",
+                                  "-t", "45",
+                                  "-n", "2",
+                                  "-o", LevelTest.TEST_OUTPUT,
+                                  TEST_ZARR_DIR,
+                                  ])
+        self.assertEqual(0, result.exit_code)
+        self.assertTrue(os.path.isdir(LevelTest.TEST_OUTPUT))
+        self.assertEqual({'0.zarr', '1.zarr'}, set(os.listdir(LevelTest.TEST_OUTPUT)))
+
+    def test_invalid_inputs(self):
+        result = self.invoke_cli(["level",
+                                  "-t", "a45",
+                                  "-o", LevelTest.TEST_OUTPUT,
+                                  TEST_NC_FILE,
+                                  ])
+        self.assertEqual("Error: Expected a pair of positive integers <tile-width>,<tile-height>\n",
+                         result.output)
+        self.assertEqual(1, result.exit_code)
+
+        result = self.invoke_cli(["level",
+                                  "-t", "-3",
+                                  "-o", LevelTest.TEST_OUTPUT,
+                                  TEST_NC_FILE,
+                                  ])
+        self.assertEqual("Error: <tile-size> must comprise positive integers\n",
+                         result.output)
+        self.assertEqual(1, result.exit_code)
+
+        result = self.invoke_cli(["level",
+                                  "-t", "45,45,45",
+                                  "-o", LevelTest.TEST_OUTPUT,
+                                  TEST_NC_FILE])
+        self.assertEqual("Error: Expected a pair of positive integers <tile-width>,<tile-height>\n",
+                         result.output)
+        self.assertEqual(1, result.exit_code)
+
+        result = self.invoke_cli(["level",
+                                  "-n", "0",
+                                  "-o", LevelTest.TEST_OUTPUT,
+                                  TEST_NC_FILE])
+        self.assertEqual("Error: <num-levels-max> must be a positive integer\n",
+                         result.output)
+        self.assertEqual(1, result.exit_code)
+
+
 class ParseTest(unittest.TestCase):
 
     def test_parse_kwargs(self):