Repository: dcs4cop/xcube (python)

Issue:
xcube server to verify valid cubes
**Is your feature request related to a problem? Please describe.**

xcube server must verify that configured cubes are valid, so that problems such as #103 are easier to understand address. 

**Describe the solution you'd like**

Call `xcube.api.assert_cube(dataset)` after opening datasets.



Pull request:
Forman 107 server to verify cubes
Closes #107.

@AliceBalfanz, to verify this PR, please test `xcube serve` with the original (illegal) cube that cause the actual problem. The exception message should now be much cleaner and point you to the actual issue.

Patch:
diff --git a/CHANGES.md b/CHANGES.md
index 415165d9..bda18153 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -2,6 +2,7 @@
 
 ### Enhancements
 
+* `xcube serve` now verifies that a configured cube is valid once it is opened. (#107)
 * Added new CLI command `xcube verify` performing data cube verification. (#19)
 * Reworked `xcube extract` to be finally useful and effective for point data extraction. (#102) 
 * `xcube server`can now filter datasets by point coordinate, e.g. `/datasets?point=12.5,52.8`. (#50) 
diff --git a/xcube/api/verify.py b/xcube/api/verify.py
index 34f2f7fc..e59352aa 100644
--- a/xcube/api/verify.py
+++ b/xcube/api/verify.py
@@ -25,7 +25,7 @@ import numpy as np
 import xarray as xr
 
 
-def assert_cube(dataset: xr.Dataset, name=None):
+def assert_cube(dataset: xr.Dataset, name=None) -> xr.Dataset:
     """
     Assert that the given *dataset* is a valid data cube.
 
@@ -40,6 +40,8 @@ def assert_cube(dataset: xr.Dataset, name=None):
         message += "- " + ";\n- ".join(report) + "."
         raise ValueError(message)
 
+    return dataset
+
 
 def verify_cube(dataset: xr.Dataset) -> List[str]:
     """
@@ -66,6 +68,7 @@ def verify_cube(dataset: xr.Dataset) -> List[str]:
     _check_time(dataset, "time", report)
     _check_lon_or_lat(dataset, "lat", -90, 90, report)
     _check_lon_or_lat(dataset, "lon", -180, 180, report)
+    # TODO (forman): verify bounds coordinate variables
     _check_data_variables(dataset, report)
     return report
 
@@ -122,6 +125,24 @@ def _check_coord_var(dataset, name, report):
         report.append(f"coordinate variable {name!r} must not be empty")
         return None
 
+    bnds_name = var.attrs.get('bounds', f'{name}_bnds')
+    if bnds_name in dataset.coords:
+        bnds_var = dataset.coords[bnds_name]
+        expected_shape = var.size, 2
+        expected_dtype = var.dtype
+        if len(bnds_var.dims) != 2 or bnds_var.dims[0] != name:
+            report.append(f"bounds coordinate variable {bnds_name!r}"
+                          f" must have dimensions ({name!r}, <bounds_dim>)")
+        if bnds_var.shape != expected_shape:
+            report.append(
+                f"shape of bounds coordinate variable {bnds_name!r}"
+                f" must be {expected_shape!r} but was {bnds_var.shape!r}")
+        if bnds_var.dtype != expected_dtype:
+            report.append(
+                f"type of bounds coordinate variable {bnds_name!r}"
+                f" must be {expected_dtype!r} but was {bnds_var.dtype!r}")
+        return None
+
     return var
 
 
@@ -137,7 +158,7 @@ def _check_lon_or_lat(dataset, name, min_value, max_value, report):
         report.append(f"values of coordinate variable {name!r}"
                       f" must be in the range {min_value} to {max_value}")
 
-    # TODO (forman): the following check is not valid for "lat" because we currently use wrong lat-order
+    # TODO (forman): the following check is not valid for "lat" because we currently use 'wrong' lat-order
     # TODO (forman): the following check is not valid for "lon" if a cube covers the antimeridian
     # if not np.all(np.diff(var.astype(np.float64)) > 0):
     #    report.append(f"values of coordinate variable {name!r} must be monotonic increasing")
diff --git a/xcube/webapi/context.py b/xcube/webapi/context.py
index 6bfc41bc..023e2c3c 100644
--- a/xcube/webapi/context.py
+++ b/xcube/webapi/context.py
@@ -32,6 +32,7 @@ import s3fs
 import xarray as xr
 import zarr
 
+from xcube.api import assert_cube
 from .cache import MemoryCacheStore, Cache
 from .defaults import DEFAULT_CMAP_CBAR, DEFAULT_CMAP_VMIN, DEFAULT_CMAP_VMAX, DEFAULT_TRACE_PERF
 from .errors import ServiceConfigError, ServiceError, ServiceBadRequestError, ServiceResourceNotFoundError
@@ -475,12 +476,12 @@ def open_ml_dataset_from_local_fs(ctx: ServiceContext, dataset_descriptor: Datas
 
     if data_format == FORMAT_NAME_NETCDF4:
         with measure_time(tag=f"opened local NetCDF dataset {path}"):
-            ds = xr.open_dataset(path)
+            ds = assert_cube(xr.open_dataset(path))
             return BaseMultiLevelDataset(ds)
 
     if data_format == FORMAT_NAME_ZARR:
         with measure_time(tag=f"opened local zarr dataset {path}"):
-            ds = xr.open_zarr(path)
+            ds = assert_cube(xr.open_zarr(path))
             return BaseMultiLevelDataset(ds)
 
     if data_format == FORMAT_NAME_LEVELS:
diff --git a/xcube/webapi/mldataset.py b/xcube/webapi/mldataset.py
index 48da453f..b7a99924 100644
--- a/xcube/webapi/mldataset.py
+++ b/xcube/webapi/mldataset.py
@@ -8,6 +8,7 @@ import xarray as xr
 import zarr
 
 from .im import TileGrid
+from ..api import assert_cube
 from ..util.geom import get_dataset_bounds
 from ..util.perf import measure_time
 
@@ -187,7 +188,7 @@ class FileStorageMultiLevelDataset(LazyMultiLevelDataset):
                     base_dir = os.path.dirname(self._dir_path)
                     level_path = os.path.join(base_dir, level_path)
         with measure_time(tag=f"opened local dataset {level_path} for level {index}"):
-            return xr.open_zarr(level_path, **zarr_kwargs)
+            return assert_cube(xr.open_zarr(level_path, **zarr_kwargs), name=level_path)
 
     def _get_tile_grid_lazily(self):
         """
@@ -258,7 +259,7 @@ class ObjectStorageMultiLevelDataset(LazyMultiLevelDataset):
         store = s3fs.S3Map(root=level_path, s3=self._obs_file_system, check=False)
         cached_store = zarr.LRUStoreCache(store, max_size=2 ** 28)
         with measure_time(tag=f"opened remote dataset {level_path} for level {index}"):
-            return xr.open_zarr(cached_store, **zarr_kwargs)
+            return assert_cube(xr.open_zarr(cached_store, **zarr_kwargs), name=level_path)
 
     def _get_tile_grid_lazily(self):
         """
@@ -379,7 +380,7 @@ class ComputedMultiLevelDataset(LazyMultiLevelDataset):
             raise self._exception_type(f"Failed to compute in-memory dataset {self._ds_id!r} at level {index} "
                                        f"from function {self._callable_name!r}: "
                                        f"expected an xarray.Dataset but got {type(computed_value)}")
-        return computed_value
+        return assert_cube(computed_value, name=self._ds_id)
 
 
 def _get_dataset_tile_grid(dataset: xr.Dataset, num_levels: int = None):


Test patch:
diff --git a/test/api/test_verify.py b/test/api/test_verify.py
index 546b59f2..f169c208 100644
--- a/test/api/test_verify.py
+++ b/test/api/test_verify.py
@@ -9,9 +9,50 @@ from xcube.api.verify import verify_cube, assert_cube
 
 class AssertAndVerifyCubeTest(unittest.TestCase):
 
-    def test_assert_cube(self):
+    def test_assert_cube_ok(self):
         cube = new_cube(variables=dict(precipitation=0.5))
-        cube["chl"] = xr.DataArray(np.random.rand(cube.dims["lat"], cube.dims["lon"]),
+        self.assertIs(cube, assert_cube(cube))
+
+    def test_assert_cube_without_bounds(self):
+        cube = new_cube(variables=dict(precipitation=0.5), drop_bounds=True)
+        self.assertIs(cube, assert_cube(cube))
+
+    def test_assert_cube_illegal_coord_var(self):
+        cube = new_cube(variables=dict(precipitation=0.5))
+        cube = cube.assign_coords(lat=xr.DataArray(np.outer(cube.lat, np.ones(cube.lon.size)),
+                                                   dims=("y", "x")),
+                                  lon=xr.DataArray(np.outer(np.ones(cube.lat.size), cube.lon),
+                                                   dims=("y", "x")))
+        with self.assertRaises(ValueError) as cm:
+            assert_cube(cube)
+        self.assertEqual("Dataset is not a valid data cube, because:\n"
+                         "- coordinate variable 'lat' must have a single dimension 'lat';\n"
+                         "- coordinate variable 'lon' must have a single dimension 'lon'.",
+                         f"{cm.exception}")
+
+    def test_assert_cube_illegal_coord_bounds_var(self):
+        cube = new_cube(variables=dict(precipitation=0.5))
+        lat_bnds = np.zeros((cube.time.size, cube.lat.size, 2))
+        lon_bnds = np.zeros((cube.time.size, cube.lon.size, 2), dtype=np.float16)
+        lat_bnds[:, :, :] = cube.lat_bnds
+        lon_bnds[:, :, :] = cube.lon_bnds
+        cube = cube.assign_coords(lat_bnds=xr.DataArray(lat_bnds, dims=("time", "lat", "bnds")),
+                                  lon_bnds=xr.DataArray(lon_bnds, dims=("time", "lon", "bnds")))
+        with self.assertRaises(ValueError) as cm:
+            assert_cube(cube)
+        self.assertEqual("Dataset is not a valid data cube, because:\n"
+                         "- bounds coordinate variable 'lat_bnds' must have dimensions ('lat', <bounds_dim>);\n"
+                         "- shape of bounds coordinate variable 'lat_bnds' must be (180, 2) but was (5, 180, 2);\n"
+                         "- bounds coordinate variable 'lon_bnds' must have dimensions ('lon', <bounds_dim>);\n"
+                         "- shape of bounds coordinate variable 'lon_bnds' must be (360, 2) but was (5, 360, 2);\n"
+                         "- type of bounds coordinate variable 'lon_bnds' must be dtype('float64')"
+                         " but was dtype('float16').",
+                         f"{cm.exception}")
+
+    def test_assert_cube_illegal_data_var(self):
+        cube = new_cube(variables=dict(precipitation=0.5))
+        shape = cube.dims["lat"], cube.dims["lon"]
+        cube["chl"] = xr.DataArray(np.random.rand(*shape),
                                    dims=("lat", "lon"),
                                    coords=dict(lat=cube.lat, lon=cube.lon))
         with self.assertRaises(ValueError) as cm: