Repository: dcs4cop/xcube

Issue:
Enhance multi-level dataset format specification
**Is your feature request related to a problem? Please describe.**

In order to allow for updating multi-level (multi-resolution) datasets (see #797), levels to be appended or updated must be generated according to the parameterisation of the original multi-level dataset.

We are currently lacking this information. 

**Describe the solution you'd like**

- [x] Enhance the [specification for multi-level datasets](https://github.com/dcs4cop/xcube/blob/master/docs/source/mldatasets.md) so the parameterisation information is persisted.
- [x] Enhance filesystem data store code to conform to the specification enhancement in [mldataset.py](https://github.com/dcs4cop/xcube/blob/master/xcube/core/store/fs/impl/mldataset.py) 


Discussion:

Patch:
diff --git a/CHANGES.md b/CHANGES.md
index de453dfc..f99ad53f 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -7,6 +7,14 @@
   It serves a single collection named "datasets" whose items are the
   datasets published by the service.
 
+* The xcube multi-level dataset specification has been enhanced. (#802)
+  - When writing multi-level datasets (`*.levels/`) we now create a new 
+    JSON file `.zlevels` that contains the parameters used to create the 
+    dataset.
+  - A new class `xcube.core.mldataset.FsMultiLevelDataset` that represents
+    a multi-level dataset persisted to some filesystem, like 
+    "file", "s3", "memory". It can also write datasets to the filesystem. 
+
 ### Fixes
 
 * The xcube server configuration parameters `url_prefix` and 
diff --git a/docs/source/mldatasets.md b/docs/source/mldatasets.md
index dd9c2c73..41d73859 100644
--- a/docs/source/mldatasets.md
+++ b/docs/source/mldatasets.md
@@ -23,8 +23,8 @@ In xcube, multi-resolution datasets are represented by the abstract class
 `xcube.core.mldataset.MultiLevelDataset`. The xcube data store framework
 refers to this datatype using the alias `mldataset`. The corresponding
 default data format is the xcube _Levels_ format, named `levels`.
-It is planned to also support Cloud Optimized GeoTIFF (COG) as format 
-for multi-resolution datasets in xcube.
+xcube also supports the Cloud Optimized GeoTIFF (COG) format 
+for reading multi-resolution datasets.
 
 The xcube Levels Format
 -----------------------
@@ -41,10 +41,12 @@ TODO (forman): link to xcube dataset convention
 
 The following is a multi-resolution dataset with three levels:
 
-    - test_pyramid.levels/
-        - 0.zarr/
-        - 1.zarr/
-        - 2.zarr/
+```text
+- test_pyramid.levels/
+    - 0.zarr/
+    - 1.zarr/
+    - 2.zarr/
+```
 
 An important use case is generating image pyramids from existing large 
 datasets without the need to create a copy of level zero.
@@ -55,10 +57,76 @@ The link file contains the path to the actual Zarr dataset
 to be used as level zero as a plain text string. It may be an absolute 
 path or a path relative to the top-level dataset.
 
-    - test_pyramid.levels/
-        - 0.link    # --> link to actual level zero dataset
-        - 1.zarr/
-        - 2.zarr/
+```text
+- test_pyramid.levels/
+    - 0.link    # --> link to actual level zero dataset
+    - 1.zarr/
+    - 2.zarr/
+```
+
+Starting with xcube 0.13.1, an additional, optional file `.zlevels` 
+has been made part of the levels format:
+
+```text
+- test_pyramid.levels/
+    - .zlevels
+    - 0.zarr/
+    - 1.zarr/
+    - 2.zarr/
+```
+
+If present, it is a text file comprising a JSON object with the following 
+properties:
+
+| Name               | Type                 | Description                                                   |
+|--------------------|----------------------|---------------------------------------------------------------|
+| `version`          | `"1.0"`              | Levels format version.                                        |
+| `num_levels`       | integer              | Number of levels in this dataset                              |
+| `use_saved_levels` | boolean              | If a next level shall be computed from the predecessor level. |
+| `tile_size`        | \[integer, integer\] | Tile size width and height in pixels.                         |
+| `agg_methods`      | object               | Mapping from variable name to aggregation method.             |
+
+Only `version` and `num_levels` are required.
+
+The properties of the `agg_methods` objects are the names of data variables.
+The values are aggregation methods. Valid values are
+
+| Value    | Description                                                  |
+|----------|--------------------------------------------------------------|
+| `first`  | Select the first pixel at (0,0) of a window of N x N pixels. | 
+| `min`    | Minimum value of a window of N x N pixels.                   | 
+| `max`    | Minimum value of a window of N x N pixels.                   | 
+| `mean`   | Mean value of a window of N x N pixels.                      | 
+| `median` | Median value of a window of N x N pixels.                    | 
+
+The following is an example of the `.zlevels` file for a dataset with the 
+data variables `CHL` (chlorophyll) if type `float32` and a variable 
+`qflags` of type `uint16`:
+
+```json
+{
+  "version": "1.0",
+  "num_levels": 8,
+  "use_saved_levels": true,
+  "tile_size": [2048, 2048],
+  "agg_methods": {
+    "CHL": "median",
+    "qflags": "first"
+  }
+}
+```
+
+---
+**xcube implementation note**: 
+When writing datasets as multi-level datasets and the `agg_methods` 
+parameter is missing, or a data variable's name is not contained in
+given `agg_methods` then `first` is used for variables that have 
+an integer data type and `median` for a floating point data type.
+In xcube Server, when opening datasets and converting them into 
+multi-level datasets on-the-fly, `agg_methods` is `first` for all 
+data variables for best performance. 
+---
+
 
 Related reads
 -------------
@@ -69,20 +137,15 @@ Related reads
   in zarr-developers / zarr-specs on GitHub.
 * [Package ndpyramid](https://github.com/carbonplan/ndpyramid)
 
+
 To be discussed
 ---------------
 
 * Allow links for all levels?
-* Make top-level directory a Zarr group (`.zgroup`)
-  and encode level metadata (e.g. `num_levels` and level links) in `.zattrs`, or 
-  even better `.zlevels`?
-* Allow a Zarr levels sub-group of the level zero Zarr dataset.
-  It would contain all levels without level zero, hence avoiding the need to 
-  link level zeros. 
-
-To do
------
-
-* Currently, the FS data stores treat relative link paths as relative
-  to the data store's `root`. See https://github.com/dcs4cop/xcube/pull/637
+* Do not write `0.link` file. Instead, provide in `.zlevels` where to find 
+  each level.
+* No longer use `.zarr` extension for levels. Just use the index as name.
+* Make top-level directory a Zarr group (`.zgroup`), so the multi-level 
+  dataset can be opened as a group using the `zarr` package.
+
 
diff --git a/xcube/core/mldataset/__init__.py b/xcube/core/mldataset/__init__.py
index 2e2b9d40..577d7ce6 100644
--- a/xcube/core/mldataset/__init__.py
+++ b/xcube/core/mldataset/__init__.py
@@ -1,5 +1,5 @@
 # The MIT License (MIT)
-# Copyright (c) 2022 by the xcube development team and contributors
+# Copyright (c) 2023 by the xcube team and contributors
 #
 # Permission is hereby granted, free of charge, to any person obtaining a
 # copy of this software and associated documentation files (the "Software"),
@@ -19,12 +19,15 @@
 # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 # DEALINGS IN THE SOFTWARE.
 
+
 from .abc import MultiLevelDataset
 from .base import BaseMultiLevelDataset
 from .combined import CombinedMultiLevelDataset
 from .computed import ComputedMultiLevelDataset
 from .computed import augment_ml_dataset
 from .computed import open_ml_dataset_from_python_code
+from .fs import FsMultiLevelDataset
+from .fs import FsMultiLevelDatasetError
 from .identity import IdentityMultiLevelDataset
 from .lazy import LazyMultiLevelDataset
 from .mapped import MappedMultiLevelDataset
diff --git a/xcube/core/mldataset/base.py b/xcube/core/mldataset/base.py
index 31d898fc..6cf90855 100644
--- a/xcube/core/mldataset/base.py
+++ b/xcube/core/mldataset/base.py
@@ -26,7 +26,7 @@ import xarray as xr
 from xcube.core.gridmapping import GridMapping
 from xcube.core.schema import rechunk_cube
 from xcube.core.subsampling import AggMethods
-from xcube.core.subsampling import assert_valid_agg_methods
+from xcube.core.subsampling import get_dataset_agg_methods
 from xcube.core.subsampling import subsample_dataset
 from xcube.core.tilingscheme import get_num_levels
 from xcube.util.assertions import assert_instance
@@ -54,23 +54,32 @@ class BaseMultiLevelDataset(LazyMultiLevelDataset):
                  base_dataset: xr.Dataset,
                  grid_mapping: Optional[GridMapping] = None,
                  num_levels: Optional[int] = None,
-                 agg_methods: AggMethods = 'first',
+                 agg_methods: Optional[AggMethods] = 'first',
                  ds_id: Optional[str] = None):
         assert_instance(base_dataset, xr.Dataset, name='base_dataset')
+        if grid_mapping is not None:
+            assert_instance(grid_mapping, GridMapping, name='grid_mapping')
 
         if grid_mapping is None:
             # TODO (forman): why not computing it lazily?
             grid_mapping = GridMapping.from_dataset(base_dataset,
                                                     tolerance=1e-4)
 
-        assert_valid_agg_methods(agg_methods)
-        self._agg_methods = agg_methods
+        self._agg_methods = get_dataset_agg_methods(
+            base_dataset,
+            xy_dim_names=grid_mapping.xy_dim_names,
+            agg_methods=agg_methods
+        )
 
         self._base_dataset = base_dataset
         super().__init__(grid_mapping=grid_mapping,
                          num_levels=num_levels,
                          ds_id=ds_id)
 
+    @property
+    def agg_methods(self):
+        return self._agg_methods
+
     def _get_num_levels_lazily(self) -> int:
         gm = self.grid_mapping
         return get_num_levels(gm.size, gm.tile_size)
diff --git a/xcube/core/mldataset/fs.py b/xcube/core/mldataset/fs.py
new file mode 100644
index 00000000..c9c22ad0
--- /dev/null
+++ b/xcube/core/mldataset/fs.py
@@ -0,0 +1,371 @@
+# The MIT License (MIT)
+# Copyright (c) 2023 by the xcube team and contributors
+#
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the "Software"),
+# to deal in the Software without restriction, including without limitation
+# the rights to use, copy, modify, merge, publish, distribute, sublicense,
+# and/or sell copies of the Software, and to permit persons to whom the
+# Software is furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+# DEALINGS IN THE SOFTWARE.
+import collections
+import json
+import math
+import pathlib
+import warnings
+from typing import Any, Optional, List, Union, Mapping, Dict, Tuple, Sequence
+
+import fsspec
+import fsspec.core
+import numpy as np
+import xarray as xr
+import zarr
+from functools import cached_property
+
+# noinspection PyUnresolvedReferences
+import xcube.core.zarrstore
+from xcube.core.gridmapping import GridMapping
+from xcube.core.subsampling import AggMethods, AggMethod
+from xcube.util.assertions import assert_instance
+from xcube.util.fspath import get_fs_path_class
+from xcube.util.fspath import resolve_path
+from xcube.util.types import ScalarOrPair
+from xcube.util.types import normalize_scalar_or_pair
+from .abc import MultiLevelDataset
+from .base import BaseMultiLevelDataset
+from .lazy import LazyMultiLevelDataset
+
+LEVELS_FORMAT_VERSION = "1.0"
+
+
+class FsMultiLevelDataset(LazyMultiLevelDataset):
+    _MIN_CACHE_SIZE = 1024 * 1024  # 1 MiB
+
+    def __init__(self,
+                 path: str,
+                 fs: Optional[fsspec.AbstractFileSystem] = None,
+                 fs_root: Optional[str] = None,
+                 fs_kwargs: Optional[Mapping[str, Any]] = None,
+                 cache_size: Optional[int] = None,
+                 consolidate: Optional[bool] = None,
+                 **zarr_kwargs):
+        if fs is None:
+            fs, path = fsspec.core.url_to_fs(path, **(fs_kwargs or {}))
+        assert_instance(fs, fsspec.AbstractFileSystem, name="fs")
+        assert_instance(path, str, name="data_id")
+        super().__init__(ds_id=path)
+        self._path = path
+        self._fs = fs
+        self._fs_root = fs_root
+        self._cache_size = cache_size
+        self._consolidate = consolidate
+        self._zarr_kwargs = zarr_kwargs
+        self._path_class = get_fs_path_class(fs)
+
+    @property
+    def path(self) -> str:
+        return self._path
+
+    @property
+    def fs(self) -> fsspec.AbstractFileSystem:
+        return self._fs
+
+    @property
+    def cache_size(self) -> Optional[int]:
+        return self._cache_size
+
+    @cached_property
+    def size_weights(self) -> np.ndarray:
+        """Size weights are used to distribute the cache size
+        over the levels.
+        """
+        return self.compute_size_weights(self.num_levels)
+
+    @cached_property
+    def tile_size(self) -> Optional[Sequence[int]]:
+        spec = self._levels_spec
+        return spec.get("tile_size")
+
+    @cached_property
+    def use_saved_levels(self) -> Optional[bool]:
+        spec = self._levels_spec
+        return spec.get("use_saved_levels")
+
+    @cached_property
+    def base_dataset_path(self) -> Optional[str]:
+        spec = self._levels_spec
+        return spec.get("base_dataset_path")
+
+    @cached_property
+    def agg_methods(self) -> Optional[Mapping[str, AggMethod]]:
+        spec = self._levels_spec
+        return spec.get("agg_methods")
+
+    @cached_property
+    def _levels_spec(self) -> Mapping[str, Any]:
+        path = f"{self._path}/.zlevels"
+        spec = {}
+        if self.fs.exists(path):
+            with self.fs.open(f"{self._path}/.zlevels") as fp:
+                spec = json.load(fp)
+            if not isinstance(spec, collections.Mapping):
+                raise TypeError('Unexpected .zlevels file.'
+                                ' Must be a JSON object.')
+        # TODO (forman): validate JSON object
+        return spec
+
+    def _get_dataset_lazily(self, index: int, parameters) \
+            -> xr.Dataset:
+
+        cache_size = self._cache_size
+
+        fs = self._fs
+
+        ds_path = self._get_path(self._path)
+        link_path = ds_path / f'{index}.link'
+        if fs.isfile(str(link_path)):
+            # If file "{index}.link" exists, we have a link to
+            # a level Zarr and open this instead,
+            with fs.open(str(link_path), 'r') as fp:
+                level_path = self._get_path(fp.read())
+            if not level_path.is_absolute() \
+                    and not self._is_path_relative_to_path(level_path,
+                                                           ds_path):
+                level_path = resolve_path(ds_path / level_path)
+        else:
+            # Nominal "{index}.zarr" must exist
+            level_path = ds_path / f'{index}.zarr'
+
+        level_zarr_store = fs.get_mapper(str(level_path))
+
+        consolidated = self._consolidate if self._consolidate is not None \
+            else ('.zmetadata' in level_zarr_store)
+
+        if isinstance(cache_size, int) \
+                and cache_size >= self._MIN_CACHE_SIZE:
+            # compute cache size for level weighted by
+            # size in pixels for each level
+            cache_size = math.ceil(self.size_weights[index] * cache_size)
+            if cache_size >= self._MIN_CACHE_SIZE:
+                level_zarr_store = zarr.LRUStoreCache(level_zarr_store,
+                                                      max_size=cache_size)
+
+        try:
+            level_dataset = xr.open_zarr(level_zarr_store,
+                                         consolidated=consolidated,
+                                         **self._zarr_kwargs)
+        except ValueError as e:
+            raise FsMultiLevelDatasetError(
+                f'Failed to open'
+                f' dataset {level_path!r}:'
+                f' {e}'
+            ) from e
+
+        level_dataset.zarr_store.set(level_zarr_store)
+        return level_dataset
+
+    @staticmethod
+    def _is_path_relative_to_path(level_path, ds_path):
+        if hasattr(level_path, 'is_relative_to'):
+            # Python >=3.9
+            return level_path.is_relative_to(ds_path)
+        try:
+            # Python <3.9
+            level_path.relative_to(ds_path)
+            return True
+        except ValueError:
+            return False
+
+    @classmethod
+    def compute_size_weights(cls, num_levels: int) -> np.ndarray:
+        weights = (2 ** np.arange(0, num_levels, dtype=np.float64)) ** 2
+        return weights[::-1] / np.sum(weights)
+
+    def _get_num_levels_lazily(self) -> int:
+        spec = self._levels_spec
+        num_levels = spec.get("num_levels")
+        levels = self._get_levels()
+        if num_levels is None:
+            num_levels = len(levels)
+        expected_levels = list(range(num_levels))
+        for level in expected_levels:
+            if level != levels[level]:
+                raise FsMultiLevelDatasetError(
+                    f'Inconsistent'
+                    f' multi-level dataset {self.ds_id!r},'
+                    f' expected levels {expected_levels!r}'
+                    f' found {levels!r}'
+                )
+        return num_levels
+
+    def _get_levels(self) -> List[int]:
+        levels = []
+        paths = [self._get_path(entry['name'])
+                 for entry in self._fs.listdir(self._path, detail=True)]
+        for path in paths:
+            # No ext, i.e. dir_name = "<level>", is proposed by
+            # https://github.com/zarr-developers/zarr-specs/issues/50.
+            # xcube already selected dir_name = "<level>.zarr".
+            basename = path.stem
+            if path.stem and path.suffix in ('', '.zarr', '.link'):
+                try:
+                    level = int(basename)
+                except ValueError:
+                    continue
+                levels.append(level)
+        levels = sorted(levels)
+        return levels
+
+    def _get_path(self, *args) -> pathlib.PurePath:
+        return self._path_class(*args)
+
+    @classmethod
+    def write_dataset(cls,
+                      dataset: Union[xr.Dataset, MultiLevelDataset],
+                      path: str,
+                      fs: Optional[fsspec.AbstractFileSystem] = None,
+                      fs_root: Optional[str] = None,
+                      fs_kwargs: Optional[Mapping[str, Any]] = None,
+                      replace: bool = False,
+                      num_levels: Optional[int] = None,
+                      consolidated: bool = True,
+                      tile_size: Optional[ScalarOrPair[int]] = None,
+                      use_saved_levels: bool = False,
+                      base_dataset_path: Optional[str] = None,
+                      agg_methods: Optional[AggMethods] = None,
+                      **zarr_kwargs) -> str:
+
+        assert_instance(dataset, (xr.Dataset, MultiLevelDataset),
+                        name='dataset')
+
+        if fs is None:
+            fs, path = fsspec.core.url_to_fs(path, **(fs_kwargs or {}))
+        if tile_size is not None:
+            tile_size = normalize_scalar_or_pair(tile_size,
+                                                 item_type=int,
+                                                 name='tile_size')
+
+        assert_instance(path, str, name='path')
+        assert_instance(fs, fsspec.AbstractFileSystem, name='fs')
+
+        if isinstance(dataset, MultiLevelDataset):
+            ml_dataset = dataset
+            if tile_size:
+                warnings.warn(
+                    'tile_size is ignored for multi-level datasets'
+                )
+            if agg_methods:
+                warnings.warn(
+                    'agg_methods is ignored for multi-level datasets'
+                )
+        else:
+            base_dataset: xr.Dataset = dataset
+            grid_mapping = None
+            if tile_size is not None:
+                grid_mapping = GridMapping.from_dataset(base_dataset)
+                x_name, y_name = grid_mapping.xy_dim_names
+                # noinspection PyTypeChecker
+                base_dataset = base_dataset.chunk({x_name: tile_size[0],
+                                                   y_name: tile_size[1]})
+                # noinspection PyTypeChecker
+                grid_mapping = grid_mapping.derive(tile_size=tile_size)
+            ml_dataset = BaseMultiLevelDataset(base_dataset,
+                                               grid_mapping=grid_mapping,
+                                               num_levels=num_levels,
+                                               agg_methods=agg_methods)
+
+        if use_saved_levels:
+            ml_dataset = BaseMultiLevelDataset(
+                ml_dataset.get_dataset(0),
+                grid_mapping=ml_dataset.grid_mapping,
+                agg_methods=agg_methods
+            )
+
+        path_class = get_fs_path_class(fs)
+        data_path = path_class(path)
+        fs.mkdirs(str(data_path), exist_ok=replace)
+
+        if num_levels is None or num_levels <= 0:
+            num_levels_max = ml_dataset.num_levels
+        else:
+            num_levels_max = min(num_levels, ml_dataset.num_levels)
+
+        with fs.open(str(data_path / '.zlevels'), mode='w') as fp:
+            levels_data: Dict[str, Any] = dict(
+                version=LEVELS_FORMAT_VERSION,
+                num_levels=num_levels_max
+            )
+            if use_saved_levels is not None:
+                levels_data.update(use_saved_levels=bool(use_saved_levels))
+            if base_dataset_path:
+                levels_data.update(base_dataset_path=base_dataset_path)
+            if tile_size is not None:
+                levels_data.update(tile_size=list(tile_size))
+            if hasattr(ml_dataset, 'agg_methods'):
+                levels_data.update(agg_methods=dict(ml_dataset.agg_methods))
+            json.dump(levels_data, fp, indent=2)
+
+        for index in range(num_levels_max):
+            level_dataset = ml_dataset.get_dataset(index)
+            if base_dataset_path and index == 0:
+                assert_instance(fs_root, str, name="fs_root")
+
+                # Write file "0.link" instead of copying
+                # level zero dataset to "0.zarr".
+
+                # Compute a relative base dataset path first
+                base_dataset_path = path_class(fs_root, base_dataset_path)
+                data_parent_path = data_path.parent
+                try:
+                    base_dataset_path = base_dataset_path.relative_to(
+                        data_parent_path
+                    )
+                except ValueError as e:
+                    raise FsMultiLevelDatasetError(
+                        f'Invalid base_dataset_id: {base_dataset_path}'
+                    ) from e
+                base_dataset_path = '..' / base_dataset_path
+
+                # Then write relative base dataset path into link file
+                link_path = data_path / f'{index}.link'
+                with fs.open(str(link_path), mode='w') as fp:
+                    fp.write(base_dataset_path.as_posix())
+            else:
+                # Write level "{index}.zarr"
+                level_path = data_path / f'{index}.zarr'
+                level_zarr_store = fs.get_mapper(str(level_path), create=True)
+                try:
+                    level_dataset.to_zarr(
+                        level_zarr_store,
+                        mode='w' if replace else None,
+                        consolidated=consolidated,
+                        **zarr_kwargs
+                    )
+                except ValueError as e:
+                    # TODO: remove already written data!
+                    raise FsMultiLevelDatasetError(
+                        f'Failed to write dataset {path}: {e}'
+                    ) from e
+                if use_saved_levels:
+                    level_dataset = xr.open_zarr(level_zarr_store,
+                                                 consolidated=consolidated)
+                    level_dataset.zarr_store.set(level_zarr_store)
+                    ml_dataset.set_dataset(index, level_dataset)
+
+        return path
+
+
+class FsMultiLevelDatasetError(ValueError):
+
+    def __init__(self, message: str):
+        super().__init__(message)
diff --git a/xcube/core/store/fs/impl/dataset.py b/xcube/core/store/fs/impl/dataset.py
index c09e6fd3..df8d5b1a 100644
--- a/xcube/core/store/fs/impl/dataset.py
+++ b/xcube/core/store/fs/impl/dataset.py
@@ -1,23 +1,23 @@
 # The MIT License (MIT)
-# Copyright (c) 2021 by the xcube development team and contributors
+# Copyright (c) 2023 by the xcube team and contributors
 #
-# Permission is hereby granted, free of charge, to any person obtaining a copy of
-# this software and associated documentation files (the "Software"), to deal in
-# the Software without restriction, including without limitation the rights to
-# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
-# of the Software, and to permit persons to whom the Software is furnished to do
-# so, subject to the following conditions:
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the "Software"),
+# to deal in the Software without restriction, including without limitation
+# the rights to use, copy, modify, merge, publish, distribute, sublicense,
+# and/or sell copies of the Software, and to permit persons to whom the
+# Software is furnished to do so, subject to the following conditions:
 #
-# The above copyright notice and this permission notice shall be included in all
-# copies or substantial portions of the Software.
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
 #
 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-# SOFTWARE.
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+# DEALINGS IN THE SOFTWARE.
 
 from abc import ABC
 from typing import Tuple, Optional
@@ -36,6 +36,7 @@ from xcube.core.zarrstore import LoggingZarrStore
 from xcube.core.zarrstore import ZarrStoreHolder
 from xcube.util.assertions import assert_instance
 from xcube.util.assertions import assert_true
+from xcube.util.fspath import is_local_fs
 from xcube.util.jsonencoder import to_json_value
 from xcube.util.jsonschema import JsonArraySchema
 from xcube.util.jsonschema import JsonBooleanSchema
@@ -45,7 +46,6 @@ from xcube.util.jsonschema import JsonObjectSchema
 from xcube.util.jsonschema import JsonStringSchema
 from xcube.util.temp import new_temp_file
 from ..accessor import FsDataAccessor
-from ..helpers import is_local_fs
 from ...datatype import DATASET_TYPE
 from ...datatype import DataType
 from ...error import DataStoreError
diff --git a/xcube/core/store/fs/impl/geodataframe.py b/xcube/core/store/fs/impl/geodataframe.py
index e2df3495..dba36f7a 100644
--- a/xcube/core/store/fs/impl/geodataframe.py
+++ b/xcube/core/store/fs/impl/geodataframe.py
@@ -1,35 +1,34 @@
 # The MIT License (MIT)
-# Copyright (c) 2021 by the xcube development team and contributors
+# Copyright (c) 2023 by the xcube team and contributors
 #
-# Permission is hereby granted, free of charge, to any person obtaining a copy of
-# this software and associated documentation files (the "Software"), to deal in
-# the Software without restriction, including without limitation the rights to
-# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
-# of the Software, and to permit persons to whom the Software is furnished to do
-# so, subject to the following conditions:
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the "Software"),
+# to deal in the Software without restriction, including without limitation
+# the rights to use, copy, modify, merge, publish, distribute, sublicense,
+# and/or sell copies of the Software, and to permit persons to whom the
+# Software is furnished to do so, subject to the following conditions:
 #
-# The above copyright notice and this permission notice shall be included in all
-# copies or substantial portions of the Software.
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
 #
 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-# SOFTWARE.
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+# DEALINGS IN THE SOFTWARE.
 
 from abc import abstractmethod
-from typing import Tuple, Optional
 
 import geopandas as gpd
 import pandas as pd
 
 from xcube.util.assertions import assert_instance
+from xcube.util.fspath import is_local_fs
 from xcube.util.jsonschema import JsonObjectSchema
 from xcube.util.temp import new_temp_file
 from ..accessor import FsDataAccessor
-from ..helpers import is_local_fs
 from ...datatype import DataType
 from ...datatype import GEO_DATA_FRAME_TYPE
 
diff --git a/xcube/core/store/fs/impl/mldataset.py b/xcube/core/store/fs/impl/mldataset.py
index b9f53d1a..4a8b8c02 100644
--- a/xcube/core/store/fs/impl/mldataset.py
+++ b/xcube/core/store/fs/impl/mldataset.py
@@ -1,37 +1,30 @@
 # The MIT License (MIT)
-# Copyright (c) 2022 by the xcube development team and contributors
+# Copyright (c) 2023 by the xcube team and contributors
 #
-# Permission is hereby granted, free of charge, to any person obtaining a copy of
-# this software and associated documentation files (the "Software"), to deal in
-# the Software without restriction, including without limitation the rights to
-# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
-# of the Software, and to permit persons to whom the Software is furnished to do
-# so, subject to the following conditions:
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the "Software"),
+# to deal in the Software without restriction, including without limitation
+# the rights to use, copy, modify, merge, publish, distribute, sublicense,
+# and/or sell copies of the Software, and to permit persons to whom the
+# Software is furnished to do so, subject to the following conditions:
 #
-# The above copyright notice and this permission notice shall be included in all
-# copies or substantial portions of the Software.
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
 #
 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-# SOFTWARE.
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+# DEALINGS IN THE SOFTWARE.
 
-import math
-import pathlib
-import warnings
-from typing import Dict, Any, List, Union, Optional
+from typing import Union
 
-import fsspec
-import numpy as np
 import xarray as xr
-import zarr
 
-from xcube.core.gridmapping import GridMapping
-from xcube.core.mldataset import BaseMultiLevelDataset
-from xcube.core.mldataset import LazyMultiLevelDataset
+from xcube.core.mldataset import FsMultiLevelDataset
+from xcube.core.mldataset import FsMultiLevelDatasetError
 from xcube.core.mldataset import MultiLevelDataset
 from xcube.core.subsampling import AGG_METHODS
 # Note, we need the following reference to register the
@@ -46,14 +39,11 @@ from xcube.util.jsonschema import JsonIntegerSchema
 from xcube.util.jsonschema import JsonNullSchema
 from xcube.util.jsonschema import JsonObjectSchema
 from xcube.util.jsonschema import JsonStringSchema
-from xcube.util.types import normalize_scalar_or_pair
 from .dataset import DatasetZarrFsDataAccessor
-from ..helpers import get_fs_path_class
-from ..helpers import resolve_path
+from ... import DataStoreError
 from ...datatype import DATASET_TYPE
 from ...datatype import DataType
 from ...datatype import MULTI_LEVEL_DATASET_TYPE
-from ...error import DataStoreError
 
 
 class MultiLevelDatasetLevelsFsDataAccessor(DatasetZarrFsDataAccessor):
@@ -71,8 +61,14 @@ class MultiLevelDatasetLevelsFsDataAccessor(DatasetZarrFsDataAccessor):
 
     def open_data(self, data_id: str, **open_params) -> MultiLevelDataset:
         assert_instance(data_id, str, name='data_id')
-        fs, root, open_params = self.load_fs(open_params)
-        return FsMultiLevelDataset(fs, root, data_id, **open_params)
+        fs, fs_root, open_params = self.load_fs(open_params)
+        try:
+            return FsMultiLevelDataset(data_id,
+                                       fs=fs,
+                                       fs_root=fs_root,
+                                       **open_params)
+        except FsMultiLevelDatasetError as e:
+            raise DataStoreError(f"{e}") from e
 
     def get_write_data_params_schema(self) -> JsonObjectSchema:
         schema = super().get_write_data_params_schema()  # creates deep copy
@@ -131,114 +127,31 @@ class MultiLevelDatasetLevelsFsDataAccessor(DatasetZarrFsDataAccessor):
                    replace: bool = False,
                    **write_params) -> str:
         assert_instance(data, MultiLevelDataset, name='data')
-        return self.write_generic_data(data,
+        return self.write_data_generic(data,
                                        data_id,
                                        replace=replace,
                                        **write_params)
 
-    def write_generic_data(self,
+    def write_data_generic(self,
                            data: Union[xr.Dataset, MultiLevelDataset],
                            data_id: str,
                            replace: bool = False,
                            **write_params) -> str:
         assert_instance(data, (xr.Dataset, MultiLevelDataset), name='data')
-        assert_instance(data_id, str, name='data_id')
-        tile_size = write_params.pop('tile_size', None)
-        agg_methods = write_params.pop('agg_methods', None)
-        if isinstance(data, MultiLevelDataset):
-            ml_dataset = data
-            if tile_size:
-                warnings.warn(
-                    'tile_size is ignored for multi-level datasets'
-                )
-            if agg_methods:
-                warnings.warn(
-                    'agg_methods is ignored for multi-level datasets'
-                )
-        else:
-            base_dataset: xr.Dataset = data
-            grid_mapping = None
-            if tile_size is not None:
-                tile_size = normalize_scalar_or_pair(
-                    tile_size, item_type=int, name='tile_size'
-                )
-                grid_mapping = GridMapping.from_dataset(base_dataset)
-                x_name, y_name = grid_mapping.xy_dim_names
-                # noinspection PyTypeChecker
-                base_dataset = base_dataset.chunk({x_name: tile_size[0],
-                                                   y_name: tile_size[1]})
-                # noinspection PyTypeChecker
-                grid_mapping = grid_mapping.derive(tile_size=tile_size)
-            ml_dataset = BaseMultiLevelDataset(base_dataset,
-                                               grid_mapping=grid_mapping,
-                                               agg_methods=agg_methods)
-        fs, root, write_params = self.load_fs(write_params)
-        consolidated = write_params.pop('consolidated', True)
-        use_saved_levels = write_params.pop('use_saved_levels', False)
-        base_dataset_id = write_params.pop('base_dataset_id', None)
-        num_levels = write_params.pop('num_levels', None)
-
-        if use_saved_levels:
-            ml_dataset = BaseMultiLevelDataset(
-                ml_dataset.get_dataset(0),
-                grid_mapping=ml_dataset.grid_mapping,
-                agg_methods=agg_methods
+        fs, fs_root, write_params = self.load_fs(write_params)
+        base_dataset_id = write_params.pop("base_dataset_id", None)
+        try:
+            return FsMultiLevelDataset.write_dataset(
+                data,
+                data_id,
+                fs=fs,
+                fs_root=fs_root,
+                replace=replace,
+                base_dataset_path=base_dataset_id,
+                **write_params
             )
-
-        path_class = get_fs_path_class(fs)
-        data_path = path_class(data_id)
-        fs.mkdirs(str(data_path), exist_ok=replace)
-
-        if num_levels is None or num_levels <= 0:
-            num_levels_max = ml_dataset.num_levels
-        else:
-            num_levels_max = min(num_levels, ml_dataset.num_levels)
-
-        for index in range(num_levels_max):
-            level_dataset = ml_dataset.get_dataset(index)
-            if base_dataset_id and index == 0:
-                # Write file "0.link" instead of copying
-                # level zero dataset to "0.zarr".
-
-                # Compute a relative base dataset path first
-                base_dataset_path = path_class(root, base_dataset_id)
-                data_parent_path = data_path.parent
-                try:
-                    base_dataset_path = base_dataset_path.relative_to(
-                        data_parent_path
-                    )
-                except ValueError as e:
-                    raise DataStoreError(
-                        f'invalid base_dataset_id: {base_dataset_id}'
-                    ) from e
-                base_dataset_path = '..' / base_dataset_path
-
-                # Then write relative base dataset path into link file
-                link_path = data_path / f'{index}.link'
-                with fs.open(str(link_path), mode='w') as fp:
-                    fp.write(base_dataset_path.as_posix())
-            else:
-                # Write level "{index}.zarr"
-                level_path = data_path / f'{index}.zarr'
-                level_zarr_store = fs.get_mapper(str(level_path), create=True)
-                try:
-                    level_dataset.to_zarr(
-                        level_zarr_store,
-                        mode='w' if replace else None,
-                        consolidated=consolidated,
-                        **write_params
-                    )
-                except ValueError as e:
-                    # TODO: remove already written data!
-                    raise DataStoreError(f'Failed to write'
-                                         f' dataset {data_id}: {e}') from e
-                if use_saved_levels:
-                    level_dataset = xr.open_zarr(level_zarr_store,
-                                                 consolidated=consolidated)
-                    level_dataset.zarr_store.set(level_zarr_store)
-                    ml_dataset.set_dataset(index, level_dataset)
-
-        return data_id
+        except FsMultiLevelDatasetError as e:
+            raise DataStoreError(f"{e}") from e
 
 
 class DatasetLevelsFsDataAccessor(MultiLevelDatasetLevelsFsDataAccessor):
@@ -260,133 +173,7 @@ class DatasetLevelsFsDataAccessor(MultiLevelDatasetLevelsFsDataAccessor):
                    replace: bool = False,
                    **write_params) -> str:
         assert_instance(data, xr.Dataset, name='data')
-        return self.write_generic_data(data,
+        return self.write_data_generic(data,
                                        data_id,
                                        replace=replace,
                                        **write_params)
-
-
-class FsMultiLevelDataset(LazyMultiLevelDataset):
-    _MIN_CACHE_SIZE = 1024 * 1024  # 1 MiB
-
-    def __init__(self,
-                 fs: fsspec.AbstractFileSystem,
-                 root: Optional[str],
-                 data_id: str,
-                 **open_params: Dict[str, Any]):
-        super().__init__(ds_id=data_id)
-        self._fs = fs
-        self._root = root
-        self._path = data_id
-        self._open_params = open_params
-        self._path_class = get_fs_path_class(fs)
-        self._size_weights: Optional[np.ndarray] = None
-
-    @property
-    def size_weights(self) -> np.ndarray:
-        if self._size_weights is None:
-            with self.lock:
-                self._size_weights = self.compute_size_weights(
-                    self.num_levels
-                )
-        return self._size_weights
-
-    def _get_dataset_lazily(self, index: int, parameters) \
-            -> xr.Dataset:
-
-        open_params = dict(self._open_params)
-
-        cache_size = open_params.pop('cache_size', None)
-
-        fs = self._fs
-
-        ds_path = self._get_path(self._path)
-        link_path = ds_path / f'{index}.link'
-        if fs.isfile(str(link_path)):
-            # If file "{index}.link" exists, we have a link to
-            # a level Zarr and open this instead,
-            with fs.open(str(link_path), 'r') as fp:
-                level_path = self._get_path(fp.read())
-            if not level_path.is_absolute() \
-                    and not self._is_path_relative_to_path(level_path,
-                                                           ds_path):
-                level_path = resolve_path(ds_path / level_path)
-        else:
-            # Nominal "{index}.zarr" must exist
-            level_path = ds_path / f'{index}.zarr'
-
-        level_zarr_store = fs.get_mapper(str(level_path))
-
-        consolidated = open_params.pop('consolidated',
-                                       '.zmetadata' in level_zarr_store)
-
-        if isinstance(cache_size, int) \
-                and cache_size >= self._MIN_CACHE_SIZE:
-            # compute cache size for level weighted by
-            # size in pixels for each level
-            cache_size = math.ceil(self.size_weights[index] * cache_size)
-            if cache_size >= self._MIN_CACHE_SIZE:
-                level_zarr_store = zarr.LRUStoreCache(level_zarr_store,
-                                                      max_size=cache_size)
-
-        try:
-            level_dataset = xr.open_zarr(level_zarr_store,
-                                         consolidated=consolidated,
-                                         **open_params)
-        except ValueError as e:
-            raise DataStoreError(f'Failed to open'
-                                 f' dataset {level_path!r}:'
-                                 f' {e}') from e
-
-        level_dataset.zarr_store.set(level_zarr_store)
-        return level_dataset
-
-    @staticmethod
-    def _is_path_relative_to_path(level_path, ds_path):
-        if hasattr(level_path, 'is_relative_to'):
-            # Python >=3.9
-            return level_path.is_relative_to(ds_path)
-        try:
-            # Python <3.9
-            level_path.relative_to(ds_path)
-            return True
-        except ValueError:
-            return False
-
-    @staticmethod
-    def compute_size_weights(num_levels: int) -> np.ndarray:
-        weights = (2 ** np.arange(0, num_levels, dtype=np.float64)) ** 2
-        return weights[::-1] / np.sum(weights)
-
-    def _get_num_levels_lazily(self) -> int:
-        levels = self._get_levels()
-        num_levels = len(levels)
-        expected_levels = list(range(num_levels))
-        for level in expected_levels:
-            if level != levels[level]:
-                raise DataStoreError(f'Inconsistent'
-                                     f' multi-level dataset {self.ds_id!r},'
-                                     f' expected levels {expected_levels!r}'
-                                     f' found {levels!r}')
-        return num_levels
-
-    def _get_levels(self) -> List[int]:
-        levels = []
-        paths = [self._get_path(entry['name'])
-                 for entry in self._fs.listdir(self._path, detail=True)]
-        for path in paths:
-            # No ext, i.e. dir_name = "<level>", is proposed by
-            # https://github.com/zarr-developers/zarr-specs/issues/50.
-            # xcube already selected dir_name = "<level>.zarr".
-            basename = path.stem
-            if path.stem and path.suffix in ('', '.zarr', '.link'):
-                try:
-                    level = int(basename)
-                except ValueError:
-                    continue
-                levels.append(level)
-        levels = sorted(levels)
-        return levels
-
-    def _get_path(self, *args) -> pathlib.PurePath:
-        return self._path_class(*args)
diff --git a/xcube/core/store/fs/store.py b/xcube/core/store/fs/store.py
index ce5f5db6..a6ef4bf6 100644
--- a/xcube/core/store/fs/store.py
+++ b/xcube/core/store/fs/store.py
@@ -1,23 +1,23 @@
 # The MIT License (MIT)
-# Copyright (c) 2022 by the xcube development team and contributors
+# Copyright (c) 2023 by the xcube team and contributors
 #
-# Permission is hereby granted, free of charge, to any person obtaining a copy of
-# this software and associated documentation files (the "Software"), to deal in
-# the Software without restriction, including without limitation the rights to
-# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
-# of the Software, and to permit persons to whom the Software is furnished to do
-# so, subject to the following conditions:
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the "Software"),
+# to deal in the Software without restriction, including without limitation
+# the rights to use, copy, modify, merge, publish, distribute, sublicense,
+# and/or sell copies of the Software, and to permit persons to whom the
+# Software is furnished to do so, subject to the following conditions:
 #
-# The above copyright notice and this permission notice shall be included in all
-# copies or substantial portions of the Software.
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
 #
 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-# SOFTWARE.
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+# DEALINGS IN THE SOFTWARE.
 
 import fnmatch
 import os.path
@@ -38,6 +38,7 @@ from xcube.util.assertions import assert_in
 from xcube.util.assertions import assert_instance
 from xcube.util.assertions import assert_true
 from xcube.util.extension import Extension
+from xcube.util.fspath import is_local_fs
 from xcube.util.jsonschema import JsonArraySchema
 from xcube.util.jsonschema import JsonBooleanSchema
 from xcube.util.jsonschema import JsonComplexSchema
@@ -47,7 +48,6 @@ from xcube.util.jsonschema import JsonObjectSchema
 from xcube.util.jsonschema import JsonStringSchema
 from .accessor import FsAccessor
 from .accessor import STORAGE_OPTIONS_PARAM_NAME
-from .helpers import is_local_fs
 from ..accessor import DataOpener
 from ..accessor import DataWriter
 from ..accessor import find_data_opener_extensions
diff --git a/xcube/core/subsampling.py b/xcube/core/subsampling.py
index 9caf0d69..c578f0f0 100644
--- a/xcube/core/subsampling.py
+++ b/xcube/core/subsampling.py
@@ -41,7 +41,7 @@ def subsample_dataset(
         dataset: xr.Dataset,
         step: int,
         xy_dim_names: Optional[Tuple[str, str]] = None,
-        agg_methods: AggMethods = None
+        agg_methods: Optional[AggMethods] = None
 ) -> xr.Dataset:
     """
     Subsample *dataset* with given integer subsampling *step*.
@@ -60,6 +60,8 @@ def subsample_dataset(
         "auto", "first", "min", "max", "mean", "median".
         If "auto", the default, "first" is used for integer variables
         and "mean" for floating point variables.
+    :return: The subsampled dataset or a tuple comprising the
+        subsampled dataset and the effective aggregation methods.
     """
     assert_instance(dataset, xr.Dataset, name='dataset')
     assert_instance(step, int, name='step')
@@ -67,11 +69,15 @@ def subsample_dataset(
 
     x_name, y_name = xy_dim_names or ('y', 'x')
 
+    agg_methods = get_dataset_agg_methods(dataset,
+                                          xy_dim_names=xy_dim_names,
+                                          agg_methods=agg_methods)
+
     new_data_vars = dict()
     new_coords = None  # used to collect coordinates from coarsen
     for var_name, var in dataset.data_vars.items():
-        if x_name in var.dims or y_name in var.dims:
-            agg_method = find_agg_method(agg_methods, var_name, var.dtype)
+        if var_name in agg_methods:
+            agg_method = agg_methods[var_name]
             if agg_method == 'first':
                 slices = get_variable_subsampling_slices(
                     var, step, xy_dim_names=xy_dim_names
@@ -125,7 +131,25 @@ def subsample_dataset(
                       attrs=dataset.attrs)
 
 
-def assert_valid_agg_methods(agg_methods: AggMethods):
+def get_dataset_agg_methods(dataset: xr.Dataset,
+                            xy_dim_names: Optional[Tuple[str, str]] = None,
+                            agg_methods: Optional[AggMethods] = None):
+    assert_instance(dataset, xr.Dataset, name='dataset')
+    assert_valid_agg_methods(agg_methods)
+
+    x_name, y_name = xy_dim_names or ('y', 'x')
+
+    dataset_agg_methods = dict()
+
+    for var_name, var in dataset.data_vars.items():
+        if x_name in var.dims or y_name in var.dims:
+            agg_method = find_agg_method(agg_methods, var_name, var.dtype)
+            dataset_agg_methods[str(var_name)] = agg_method
+
+    return dataset_agg_methods
+
+
+def assert_valid_agg_methods(agg_methods: Optional[AggMethods]):
     """Assert that the given *agg_methods* are valid."""
     assert_instance(agg_methods,
                     (type(None), str, collections.abc.Mapping),
diff --git a/xcube/core/store/fs/helpers.py b/xcube/util/fspath.py
similarity index 76%
rename from xcube/core/store/fs/helpers.py
rename to xcube/util/fspath.py
index 50a86a99..9af67776 100644
--- a/xcube/core/store/fs/helpers.py
+++ b/xcube/util/fspath.py
@@ -1,23 +1,23 @@
 # The MIT License (MIT)
-# Copyright (c) 2022 by the xcube development team and contributors
+# Copyright (c) 2023 by the xcube team and contributors
 #
-# Permission is hereby granted, free of charge, to any person obtaining a copy of
-# this software and associated documentation files (the "Software"), to deal in
-# the Software without restriction, including without limitation the rights to
-# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
-# of the Software, and to permit persons to whom the Software is furnished to do
-# so, subject to the following conditions:
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the "Software"),
+# to deal in the Software without restriction, including without limitation
+# the rights to use, copy, modify, merge, publish, distribute, sublicense,
+# and/or sell copies of the Software, and to permit persons to whom the
+# Software is furnished to do so, subject to the following conditions:
 #
-# The above copyright notice and this permission notice shall be included in all
-# copies or substantial portions of the Software.
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
 #
 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-# SOFTWARE.
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+# DEALINGS IN THE SOFTWARE.
 
 import pathlib
 from typing import Type, Iterator


Test patch:
diff --git a/test/cli/test_level.py b/test/cli/test_level.py
index 376f953b..630fee16 100644
--- a/test/cli/test_level.py
+++ b/test/cli/test_level.py
@@ -48,7 +48,8 @@ class LevelDataTest(CliDataTest):
         self._assert_result_ok(
             result, [((1, 1, 1, 1, 1), (45, 45, 45, 45), (90, 90, 90, 90)),
                      ((1, 1, 1, 1, 1), (45, 45), (90, 90)),
-                     ((1, 1, 1, 1, 1), (45,), (90,))],
+                     ((1, 1, 1, 1, 1), (45,), (90,)),
+                     ((1, 1, 1, 1, 1), (23,), (45,))],
             LevelDataTest.TEST_OUTPUT,
             'Multi-level dataset written to /test.levels after .*\n'
         )
@@ -89,7 +90,7 @@ class LevelDataTest(CliDataTest):
                                   ])
         self.assertEqual(0, result.exit_code)
         self.assertTrue(os.path.isdir(LevelDataTest.TEST_OUTPUT))
-        self.assertEqual({'0.zarr', '1.zarr', '2.zarr'},
+        self.assertEqual({'.zlevels', '0.zarr', '1.zarr', '2.zarr'},
                          set(os.listdir(LevelDataTest.TEST_OUTPUT)))
 
     def test_level_with_zarr(self):
@@ -100,7 +101,7 @@ class LevelDataTest(CliDataTest):
                                   ])
         self.assertEqual(0, result.exit_code)
         self.assertTrue(os.path.isdir(LevelDataTest.TEST_OUTPUT))
-        self.assertEqual({'0.zarr', '1.zarr', '2.zarr'},
+        self.assertEqual({'.zlevels', '0.zarr', '1.zarr', '2.zarr'},
                          set(os.listdir(LevelDataTest.TEST_OUTPUT)))
 
     def test_level_with_zarr_agg_method(self):
@@ -112,7 +113,7 @@ class LevelDataTest(CliDataTest):
                                   ])
         self.assertEqual(0, result.exit_code)
         self.assertTrue(os.path.isdir(LevelDataTest.TEST_OUTPUT))
-        self.assertEqual({'0.zarr', '1.zarr', '2.zarr'},
+        self.assertEqual({'.zlevels', '0.zarr', '1.zarr', '2.zarr'},
                          set(os.listdir(LevelDataTest.TEST_OUTPUT)))
 
     def test_level_with_zarr_agg_methods(self):
@@ -125,7 +126,7 @@ class LevelDataTest(CliDataTest):
         ])
         self.assertEqual(0, result.exit_code)
         self.assertTrue(os.path.isdir(LevelDataTest.TEST_OUTPUT))
-        self.assertEqual({'0.zarr', '1.zarr', '2.zarr'},
+        self.assertEqual({'.zlevels', '0.zarr', '1.zarr', '2.zarr'},
                          set(os.listdir(LevelDataTest.TEST_OUTPUT)))
 
     def test_level_with_zarr_link(self):
@@ -137,7 +138,7 @@ class LevelDataTest(CliDataTest):
                                   ])
         self.assertEqual(0, result.exit_code)
         self.assertTrue(os.path.isdir(LevelDataTest.TEST_OUTPUT))
-        self.assertEqual({'0.link', '1.zarr', '2.zarr'},
+        self.assertEqual({'.zlevels', '0.link', '1.zarr', '2.zarr'},
                          set(os.listdir(LevelDataTest.TEST_OUTPUT)))
 
     def test_level_with_zarr_num_levels_max(self):
@@ -149,7 +150,7 @@ class LevelDataTest(CliDataTest):
                                   ])
         self.assertEqual(0, result.exit_code)
         self.assertTrue(os.path.isdir(LevelDataTest.TEST_OUTPUT))
-        self.assertEqual({'0.zarr', '1.zarr'},
+        self.assertEqual({'.zlevels', '0.zarr', '1.zarr'},
                          set(os.listdir(LevelDataTest.TEST_OUTPUT)))
 
     def test_invalid_inputs(self):
diff --git a/test/core/mldataset/test_base.py b/test/core/mldataset/test_base.py
index c3d8a3c9..050d3bce 100644
--- a/test/core/mldataset/test_base.py
+++ b/test/core/mldataset/test_base.py
@@ -15,6 +15,7 @@ class BaseMultiLevelDatasetTest(unittest.TestCase):
         self.assertIsInstance(ml_ds.grid_mapping, GridMapping)
         self.assertIsNotNone(ml_ds.lock)
         self.assertEqual(3, ml_ds.num_levels)
+        self.assertEqual({'noise': 'first'}, ml_ds.agg_methods)
 
     def test_resolutions(self):
         ds = get_test_dataset()
diff --git a/test/core/mldataset/test_fs.py b/test/core/mldataset/test_fs.py
new file mode 100644
index 00000000..5860d28c
--- /dev/null
+++ b/test/core/mldataset/test_fs.py
@@ -0,0 +1,168 @@
+# The MIT License (MIT)
+# Copyright (c) 2023 by the xcube team and contributors
+#
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the "Software"),
+# to deal in the Software without restriction, including without limitation
+# the rights to use, copy, modify, merge, publish, distribute, sublicense,
+# and/or sell copies of the Software, and to permit persons to whom the
+# Software is furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+# DEALINGS IN THE SOFTWARE.
+
+
+import math
+import unittest
+from typing import Optional, List, Mapping
+
+import fsspec
+import fsspec.core
+import xarray as xr
+
+from xcube.core.mldataset import FsMultiLevelDataset
+from xcube.core.new import new_cube
+from xcube.core.subsampling import AggMethod
+
+
+class FsMultiLevelDatasetTest(unittest.TestCase):
+
+    def setUp(self) -> None:
+        self.dataset = new_cube(
+            width=512,
+            height=256,
+            x_res=360 / 512,
+            y_res=180 / 256,
+            time_periods=1,
+            variables=dict(CHL=0.8, qflags=1)
+        )
+
+        self.fs: fsspec.AbstractFileSystem = fsspec.core.get_filesystem_class(
+            "memory"
+        )()
+
+    def tearDown(self) -> None:
+        self.fs.rm("/", recursive=True)
+
+    def test_io_nl_4_ts_256(self):
+        self.assert_io_ok(
+            "test.levels",
+            num_levels=4,
+            agg_methods=None,
+            tile_size=256,
+            use_saved_levels=False,
+            base_dataset_path=None,
+            expected_files=[
+                ".zlevels",
+                "0.zarr",
+                "1.zarr",
+                "2.zarr",
+                "3.zarr"
+            ],
+            expected_num_levels=4,
+            expected_tile_size=[256, 256],
+            expected_agg_methods={'CHL': 'mean', 'qflags': 'first'},
+        )
+
+    def test_io_nl_4_ts_256_agg(self):
+        self.assert_io_ok(
+            "test.levels",
+            num_levels=4,
+            agg_methods={'CHL': 'median', 'qflags': 'max'},
+            tile_size=256,
+            use_saved_levels=False,
+            base_dataset_path=None,
+            expected_files=[
+                ".zlevels",
+                "0.zarr",
+                "1.zarr",
+                "2.zarr",
+                "3.zarr"
+            ],
+            expected_num_levels=4,
+            expected_tile_size=[256, 256],
+            expected_agg_methods={'CHL': 'median', 'qflags': 'max'},
+        )
+
+    def test_io_nl_4_ts_256_base(self):
+        self.dataset.to_zarr(self.fs.get_mapper("test.zarr"))
+        self.assert_io_ok(
+            "test.levels",
+            num_levels=4,
+            agg_methods=None,
+            tile_size=256,
+            use_saved_levels=False,
+            base_dataset_path="test.zarr",
+            expected_files=[
+                ".zlevels",
+                "0.link",
+                "1.zarr",
+                "2.zarr",
+                "3.zarr"
+            ],
+            expected_num_levels=4,
+            expected_tile_size=[256, 256],
+            expected_agg_methods={'CHL': 'mean', 'qflags': 'first'},
+        )
+
+    def assert_io_ok(
+            self,
+            path: str,
+            num_levels: Optional[int],
+            agg_methods: Optional[Mapping[str, AggMethod]],
+            tile_size: Optional[int],
+            use_saved_levels: bool,
+            base_dataset_path: Optional[str],
+            expected_files: List[str],
+            expected_num_levels: int,
+            expected_agg_methods: Optional[Mapping[str, AggMethod]],
+            expected_tile_size: List[int]
+    ):
+        fs = self.fs
+        FsMultiLevelDataset.write_dataset(self.dataset,
+                                          path,
+                                          fs=fs,
+                                          fs_root="",
+                                          replace=True,
+                                          num_levels=num_levels,
+                                          agg_methods=agg_methods,
+                                          tile_size=tile_size,
+                                          use_saved_levels=use_saved_levels,
+                                          base_dataset_path=base_dataset_path)
+        self.assertTrue(fs.isdir(path))
+        self.assertEqual(set([f"/{path}/{f}" for f in expected_files]),
+                         set(fs.listdir(path, detail=False)))
+
+        ml_dataset = FsMultiLevelDataset(path, fs=fs)
+        self.assertEquals(expected_num_levels, ml_dataset.num_levels)
+        self.assertEquals(expected_agg_methods, ml_dataset.agg_methods)
+        self.assertEquals(expected_tile_size, ml_dataset.tile_size)
+        self.assertEquals(use_saved_levels, ml_dataset.use_saved_levels)
+        self.assertEquals(base_dataset_path, ml_dataset.base_dataset_path)
+        self.assertEquals(None, ml_dataset.cache_size)
+        self.assertEquals(num_levels, len(ml_dataset.size_weights))
+        self.assertTrue(all(ml_dataset.size_weights > 0.0))
+        for i in range(num_levels):
+            self.assertIsInstance(ml_dataset.get_dataset(i), xr.Dataset)
+
+    def test_compute_size_weights(self):
+        size = 2 ** 28
+        weighted_sizes = list(map(
+            math.ceil, size * FsMultiLevelDataset.compute_size_weights(5)
+        ))
+        self.assertEqual(
+            [201523393,
+             50380849,
+             12595213,
+             3148804,
+             787201],
+            weighted_sizes
+        )
diff --git a/test/core/store/fs/impl/test_mldataset.py b/test/core/store/fs/impl/test_mldataset.py
deleted file mode 100644
index 13f7eac8..00000000
--- a/test/core/store/fs/impl/test_mldataset.py
+++ /dev/null
@@ -1,50 +0,0 @@
-#  The MIT License (MIT)
-#  Copyright (c) 2022 by the xcube development team and contributors
-#
-#  Permission is hereby granted, free of charge, to any person obtaining a
-#  copy of this software and associated documentation files (the "Software"),
-#  to deal in the Software without restriction, including without limitation
-#  the rights to use, copy, modify, merge, publish, distribute, sublicense,
-#  and/or sell copies of the Software, and to permit persons to whom the
-#  Software is furnished to do so, subject to the following conditions:
-#
-#  The above copyright notice and this permission notice shall be included in
-#  all copies or substantial portions of the Software.
-#
-#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-#  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-#  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
-#  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
-#  DEALINGS IN THE SOFTWARE.
-
-import math
-#
-#  Permission is hereby granted, free of charge, to any person obtaining a
-#  copy of this software and associated documentation files (the "Software"),
-#  to deal in the Software without restriction, including without limitation
-#  the rights to use, copy, modify, merge, publish, distribute, sublicense,
-#  and/or sell copies of the Software, and to permit persons to whom the
-#  Software is furnished to do so, subject to the following conditions:
-#
-#
-import unittest
-
-from xcube.core.store.fs.impl.mldataset import FsMultiLevelDataset
-
-
-class FsMultiLevelDatasetTest(unittest.TestCase):
-    def test_compute_size_weights(self):
-        size = 2 ** 28
-        weighted_sizes = list(map(
-            math.ceil, size * FsMultiLevelDataset.compute_size_weights(5)
-        ))
-        self.assertEqual(
-            [201523393,
-             50380849,
-             12595213,
-             3148804,
-             787201],
-            weighted_sizes
-        )
diff --git a/test/core/test_subsampling.py b/test/core/test_subsampling.py
index d23a0410..88b0ec16 100644
--- a/test/core/test_subsampling.py
+++ b/test/core/test_subsampling.py
@@ -26,6 +26,7 @@ import xarray as xr
 
 from xcube.core.new import new_cube
 from xcube.core.subsampling import find_agg_method
+from xcube.core.subsampling import get_dataset_agg_methods
 from xcube.core.subsampling import subsample_dataset
 
 
@@ -53,7 +54,6 @@ class SubsampleDatasetTest(unittest.TestCase):
                                 variables=dict(var_1=test_data_1,
                                                var_2=test_data_2))
 
-
     def test_subsample_dataset_none(self):
         subsampled_dataset = subsample_dataset(self.dataset,
                                                step=2,
@@ -82,6 +82,7 @@ class SubsampleDatasetTest(unittest.TestCase):
             np.array([2., 4., 6.]),
         )
 
+
     def test_subsample_dataset_first(self):
         subsampled_dataset = subsample_dataset(self.dataset,
                                                step=2,
@@ -216,6 +217,22 @@ class SubsampleDatasetTest(unittest.TestCase):
             expected_var_2
         )
 
+    def test_get_dataset_agg_methods(self):
+        agg_methods = get_dataset_agg_methods(self.dataset,
+                                              agg_methods=None)
+        self.assertEqual({'var_1': 'first', 'var_2': 'mean'},
+                         agg_methods)
+
+        agg_methods = get_dataset_agg_methods(self.dataset,
+                                              agg_methods='mean')
+        self.assertEqual({'var_1': 'mean', 'var_2': 'mean'},
+                         agg_methods)
+
+        agg_methods = get_dataset_agg_methods(self.dataset,
+                                              agg_methods='max')
+        self.assertEqual({'var_1': 'max', 'var_2': 'max'},
+                         agg_methods)
+
     # noinspection PyTypeChecker
     def test_find_agg_method(self):
         for m in ('first', 'min', 'max', 'mean', 'median'):
diff --git a/test/core/store/fs/test_helpers.py b/test/util/test_fspath.py
similarity index 93%
rename from test/core/store/fs/test_helpers.py
rename to test/util/test_fspath.py
index 751f8625..b89581d9 100644
--- a/test/core/store/fs/test_helpers.py
+++ b/test/util/test_fspath.py
@@ -4,11 +4,22 @@ import unittest
 import fsspec
 import fsspec.implementations.local
 
-from xcube.core.store.fs.helpers import get_fs_path_class
-from xcube.core.store.fs.helpers import resolve_path
+from xcube.util.fspath import get_fs_path_class
+from xcube.util.fspath import resolve_path
 
 
-class FsHelpersTest(unittest.TestCase):
+class FsPathTest(unittest.TestCase):
+    def test_get_fs_path_class(self):
+        self.assertIs(pathlib.PurePath,
+                      get_fs_path_class(
+                          fsspec.implementations.local.LocalFileSystem()
+                      ))
+        self.assertIs(pathlib.PurePath,
+                      get_fs_path_class(fsspec.filesystem("file")))
+        self.assertIs(pathlib.PurePosixPath,
+                      get_fs_path_class(fsspec.filesystem("s3")))
+        self.assertIs(pathlib.PurePosixPath,
+                      get_fs_path_class(fsspec.filesystem("memory")))
 
     def test_resolve_path_succeeds(self):
         path = pathlib.PurePosixPath("a", "b", "c")
@@ -48,14 +59,3 @@ class FsHelpersTest(unittest.TestCase):
         self.assertEqual('cannot resolve path, misplaced ".."',
                          f'{cm.exception}')
 
-    def test_get_fs_path_class(self):
-        self.assertIs(pathlib.PurePath,
-                      get_fs_path_class(
-                          fsspec.implementations.local.LocalFileSystem()
-                      ))
-        self.assertIs(pathlib.PurePath,
-                      get_fs_path_class(fsspec.filesystem("file")))
-        self.assertIs(pathlib.PurePosixPath,
-                      get_fs_path_class(fsspec.filesystem("s3")))
-        self.assertIs(pathlib.PurePosixPath,
-                      get_fs_path_class(fsspec.filesystem("memory")))