Repository: dcs4cop/xcube

Issue:
Implement OGC API - Coverages
**Is your feature request related to a problem? Please describe.**
We need an implementation of OGC API - Coverages for OGC Testbed-19.

**Describe the solution you'd like**
`xcube serve` should provide a standard-compliant OGC API - Coverages instance along with its other APIs.

**Describe alternatives you've considered**
 - Use the reference Python implementation in [pygeoapi](https://pygeoapi.io/) -- but this would be hard to integrate with xcube's existing server framework and pull in various additional dependencies.
 - Don't implement the API -- but this would largely exclude us from Testbed-19 participation, since OGC API - Coverages is a core component of the draft Geodatacube API.

**Additional context**
 - [Draft GDC API](https://m-mohr.github.io/geodatacube-api/)
 - [Testbed-19 call for participation(https://portal.ogc.org/files/?artifact_id=104098)


Discussion:

Patch:
diff --git a/CHANGES.md b/CHANGES.md
index eba84ccc..cdc4416d 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -2,6 +2,8 @@
 
 ### Enhancements
 
+* Added a basic implementation of the draft version of OGC API - Coverages
+
 ### Fixes
 
 ### Other changes
diff --git a/docs/source/webapi.rst b/docs/source/webapi.rst
index 1f4c6d28..02b340d5 100644
--- a/docs/source/webapi.rst
+++ b/docs/source/webapi.rst
@@ -18,7 +18,8 @@ API, clients can
 * Get places (GeoJSON features including vector data) that can be associated
   with xcube datasets;
 * Perform compute operations on datasets, with the results calculated on
-  demand and presented as new, dynamically generated datasets.
+  demand and presented as new, dynamically generated datasets;
+* Retrieve coverages backed by datasets using the OGC API - Coverages API.
 
 Later versions of API will also allow for xcube dataset management including
 generation, modification, and deletion of xcube datasets.
diff --git a/xcube/plugin.py b/xcube/plugin.py
index 349fa3db..083cd977 100644
--- a/xcube/plugin.py
+++ b/xcube/plugin.py
@@ -229,6 +229,7 @@ def _register_server_apis(ext_registry: extension.ExtensionRegistry):
         'tiles',
         'timeseries',
         'volumes',
+        'ows.coverages',
         'ows.stac',
         'ows.wmts',
         's3',
diff --git a/xcube/webapi/ows/coverages/__init__.py b/xcube/webapi/ows/coverages/__init__.py
new file mode 100644
index 00000000..d003a730
--- /dev/null
+++ b/xcube/webapi/ows/coverages/__init__.py
@@ -0,0 +1,22 @@
+# 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.
+
+from .routes import api
diff --git a/xcube/webapi/ows/coverages/api.py b/xcube/webapi/ows/coverages/api.py
new file mode 100644
index 00000000..170e6dce
--- /dev/null
+++ b/xcube/webapi/ows/coverages/api.py
@@ -0,0 +1,31 @@
+# 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.
+
+from xcube.server.api import Api
+from xcube.webapi.ows.coverages.context import CoveragesContext
+
+api = Api(
+    'ows.coverages',
+    description='OGC API - Coverages',
+    required_apis=['datasets'],
+    create_ctx=CoveragesContext,
+    config_schema=None,
+)
diff --git a/xcube/webapi/ows/coverages/context.py b/xcube/webapi/ows/coverages/context.py
new file mode 100644
index 00000000..6ea47818
--- /dev/null
+++ b/xcube/webapi/ows/coverages/context.py
@@ -0,0 +1,38 @@
+# 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.
+
+
+from xcube.server.api import Context
+from xcube.webapi.common.context import ResourcesContext
+
+from ...datasets.context import DatasetsContext
+
+
+class CoveragesContext(ResourcesContext):
+    def __init__(self, server_ctx: Context):
+        super().__init__(server_ctx)
+        self._datasets_ctx = server_ctx.get_api_ctx(
+            'datasets', cls=DatasetsContext
+        )
+
+    @property
+    def datasets_ctx(self) -> DatasetsContext:
+        return self._datasets_ctx
diff --git a/xcube/webapi/ows/coverages/controllers.py b/xcube/webapi/ows/coverages/controllers.py
new file mode 100644
index 00000000..39045f11
--- /dev/null
+++ b/xcube/webapi/ows/coverages/controllers.py
@@ -0,0 +1,314 @@
+# 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 os
+import tempfile
+from typing import Mapping, Sequence, Optional
+
+import numpy as np
+import pyproj
+import xarray as xr
+
+from xcube.webapi.datasets.context import DatasetsContext
+
+
+def get_coverage_as_json(ctx: DatasetsContext, collection_id: str):
+    """
+    Return a JSON representation of the specified coverage
+
+    Currently, the range set component is omitted.
+
+    :param ctx: a dataset context
+    :param collection_id: the ID of the collection providing the coverage
+    :return: a JSON representation of the coverage
+    """
+    return {
+        'id': collection_id,
+        'type': 'CoverageByDomainAndRange',
+        'envelope': get_collection_envelope(ctx, collection_id),
+        'domainSet': get_coverage_domainset(ctx, collection_id),
+        'rangeSet': {
+            'type': 'RangeSet',
+            # TODO: Wait for next update to API specification before
+            #  implementing the data block -- not clear yet whether this
+            #  is being deprecated along with the rangeSet endpoint.
+            'dataBlock': {'type': 'VDataBlock', 'values': ['TODO']},
+        },
+        'rangeType': get_coverage_rangetype(ctx, collection_id),
+        'metadata': get_collection_metadata(ctx, collection_id),
+    }
+
+
+def get_coverage_data(
+    ctx: DatasetsContext,
+    collection_id: str,
+    query: Mapping[str, Sequence[str]],
+    content_type: str,
+) -> Optional[bytes]:
+    """
+    Return coverage data from a dataset
+
+    This method currently returns coverage data from a dataset as either
+    TIFF or NetCDF. The bbox, datetime, and properties parameters are
+    currently handled.
+
+    :param ctx: a datasets context
+    :param collection_id: the dataset from which to return the coverage
+    :param query: the HTTP query parameters
+    :param content_type: the MIME type of the desired output format
+    :return: the coverage as bytes in the requested output format, or None
+             if the requested output format is not supported
+    """
+    ds = get_dataset(ctx, collection_id)
+    if 'bbox' in query:
+        bbox = list(map(float, query['bbox'][0].split(',')))
+        ds = ds.sel(lat=slice(bbox[0], bbox[2]), lon=slice(bbox[1], bbox[3]))
+    if 'datetime' in query:
+        timespec = query['datetime']
+        if '/' in timespec:
+            timefrom, timeto = timespec[0].split('/')
+            ds = ds.sel(time=slice(timefrom, timeto))
+        else:
+            ds = ds.sel(time=timespec, method='nearest').squeeze()
+    if 'properties' in query:
+        vars_to_keep = set(query['properties'][0].split(','))
+        data_vars = set(ds.data_vars)
+        vars_to_drop = list(data_vars - vars_to_keep)
+        ds = ds.drop_vars(vars_to_drop)
+    if content_type in ['image/tiff', 'application/x-geotiff']:
+        return dataset_to_tiff(ds)
+    if content_type in ['application/netcdf', 'application/x-netcdf']:
+        return dataset_to_netcdf(ds)
+    return None
+
+
+def dataset_to_tiff(ds: xr.Dataset) -> bytes:
+    """
+    Return an in-memory TIFF representing a dataset
+
+    :param ds: a dataset
+    :return: TIFF-formatted bytes representing the dataset
+    """
+    with tempfile.TemporaryDirectory() as tempdir:
+        path = os.path.join(tempdir, 'out.tiff')
+        ds.rio.to_raster(path)
+        with open(path, 'rb') as fh:
+            data = fh.read()
+    return data
+
+
+def dataset_to_netcdf(ds: xr.Dataset) -> bytes:
+    """
+    Return an in-memory NetCDF representing a dataset
+
+    :param ds: a dataset
+    :return: NetCDF-formatted bytes representing the dataset
+    """
+    with tempfile.TemporaryDirectory() as tempdir:
+        path = os.path.join(tempdir, 'out.nc')
+        ds.to_netcdf(path)
+        with open(path, 'rb') as fh:
+            data = fh.read()
+    return data
+
+
+def get_coverage_domainset(ctx: DatasetsContext, collection_id: str):
+    """
+    Return the domain set of a dataset-backed coverage
+
+    The domain set is the set of input parameters (e.g. geographical extent,
+    time span) for which a coverage is defined.
+
+    :param ctx: a datasets context
+    :param collection_id: the dataset for which to return the domain set
+    :return: a dictionary representing an OGC API - Coverages domain set
+    """
+    ds = get_dataset(ctx, collection_id)
+    grid_limits = dict(
+        type='GridLimits',
+        srsName=f'http://www.opengis.net/def/crs/OGC/0/Index{len(ds.dims)}D',
+        axisLabels=list(ds.dims),
+        axis=[_get_grid_limits_axis(ds, dim) for dim in ds.dims],
+    )
+    grid = dict(
+        type='GeneralGridCoverage',
+        srsName=get_crs_from_dataset(ds),
+        axisLabels=list(ds.dims.keys()),
+        axis=_get_axes_properties(ds),
+        gridLimits=grid_limits,
+    )
+    return dict(type='DomainSet', generalGrid=grid)
+
+
+def get_collection_metadata(ctx: DatasetsContext, collection_id: str):
+    """
+    Return a metadata dictionary for a dataset
+
+    The metadata is taken directly from the dataset attributes.
+
+    :param ctx: a datasets context
+    :param collection_id: the dataset for which to return the metadata
+    :return: a dictionary of metadata keys and values
+    """
+    ds = get_dataset(ctx, collection_id)
+    return ds.attrs
+
+
+def get_dataset(ctx: DatasetsContext, collection_id: str):
+    """
+    Get a dataset from a datasets context
+
+    :param ctx: a datasets context
+    :param collection_id: the ID of a dataset in the context
+    :return: the dataset
+    """
+    ml_dataset = ctx.get_ml_dataset(collection_id)
+    ds = ml_dataset.get_dataset(0)
+    assert isinstance(ds, xr.Dataset)
+    return ds
+
+
+def _get_axes_properties(ds: xr.Dataset) -> list[dict]:
+    return [_get_axis_properties(ds, dim) for dim in ds.dims]
+
+
+def _get_axis_properties(ds: xr.Dataset, dim: str) -> dict:
+    axis = ds.coords[dim]
+    return dict(
+        type='RegularAxis',
+        axisLabel=dim,
+        lowerBound=axis[0].item(),
+        upperBound=axis[-1].item(),
+        resolution=abs((axis[-1] - axis[0]).item() / len(axis)),
+        uomLabel=_get_units(ds, dim),
+    )
+
+
+def _get_grid_limits_axis(ds: xr.Dataset, dim: str):
+    return dict(
+        type='IndexAxis', axisLabel=dim, lowerBound=0, upperBound=len(ds[dim])
+    )
+
+
+def _get_units(ds: xr.Dataset, dim: str):
+    coord = ds.coords[dim]
+    if hasattr(coord, 'attrs') and 'units' in coord.attrs:
+        return coord.attrs['units']
+    else:
+        return 'degrees'
+
+
+def get_crs_from_dataset(ds: xr.Dataset) -> str:
+    """
+    Return the CRS of a dataset as a string. The CRS is taken from the
+    metadata of the crs or spatial_ref variables, if available.
+    "EPSG:4326" is used as a fallback.
+
+    :param ds: a dataset
+    :return: a string representation of the dataset's CRS, or "EPSG:4326"
+             if the CRS cannot be determined
+    """
+    for var_name in 'crs', 'spatial_ref':
+        if var_name in ds.variables:
+            var = ds[var_name]
+            if 'spatial_ref' in var.attrs:
+                crs_string = ds[var_name].attrs['spatial_ref']
+                return pyproj.crs.CRS(crs_string).to_string()
+    return 'EPSG:4326'
+
+
+def get_coverage_rangetype(ctx: DatasetsContext, collection_id: str) -> dict[str, list]:
+    """
+    Return the range type of a dataset
+
+    The range type describes the data types of the dataset's variables
+    """
+    ds = get_dataset(ctx, collection_id)
+    result = dict(type='DataRecord', field=[])
+    for var_name in ds.data_vars:
+        result['field'].append(
+            dict(
+                type='Quantity',
+                name=var_name,
+                description=get_dataarray_description(ds[var_name]),
+                encodingInfo=dict(
+                    dataType=dtype_to_opengis_datatype(ds[var_name].dtype)
+                ),
+            )
+        )
+    return result
+
+
+def dtype_to_opengis_datatype(dt: np.dtype) -> str:
+    """
+    Convert a NumPy dtype to an equivalent OpenGIS type identifier string.
+
+    :param dt: a NumPy dtype
+    :return: an equivalent OpenGIS type identifier string, or an empty string
+             if the dtype is not recognized
+    """
+    nbits = 8 * np.dtype(dt).itemsize
+    int_size_map = {8: 'Byte', 16: 'Short', 32: 'Int', 64: 'Long'}
+    prefix = 'http://www.opengis.net/def/dataType/OGC/0/'
+    if np.issubdtype(dt, np.floating):
+        opengis_type = f'{prefix}float{nbits}'
+    elif np.issubdtype(dt, np.signedinteger):
+        opengis_type = f'{prefix}signed{int_size_map[nbits]}'
+    elif np.issubdtype(dt, np.unsignedinteger):
+        opengis_type = f'{prefix}unsigned{int_size_map[nbits]}'
+    elif 'datetime64' in str(dt):
+        opengis_type = 'http://www.opengis.net/def/bipm/UTC'
+    else:
+        opengis_type = ''  # TODO decide what to do in this case
+    return opengis_type
+
+
+def get_dataarray_description(da: xr.DataArray) -> str:
+    """
+    Return a string describing a DataArray, either from an attribute or,
+    as a fallback, from its name attribute.
+
+    :param da: a DataArray
+    :return: a string describing the DataArray
+    """
+    if hasattr(da, 'attrs'):
+        for attr in ['description', 'long_name', 'standard_name', 'name']:
+            if attr in da.attrs:
+                return da.attrs[attr]
+    return str(da.name)
+
+
+def get_collection_envelope(ds_ctx, collection_id):
+    """
+    Return the OGC API - Coverages envelope of a dataset.
+
+    The envelope comprises the extents of all the dataset's dimensions.
+
+    :param ds_ctx: a datasets context
+    :param collection_id: a dataset ID within the given context
+    :return: the envelope of the specified dataset
+    """
+    ds = get_dataset(ds_ctx, collection_id)
+    return {
+        'type': 'EnvelopeByAxis',
+        'srsName': get_crs_from_dataset(ds),
+        'axisLabels': list(ds.dims.keys()),
+        'axis': _get_axes_properties(ds),
+    }
diff --git a/xcube/webapi/ows/coverages/routes.py b/xcube/webapi/ows/coverages/routes.py
new file mode 100644
index 00000000..a842a11a
--- /dev/null
+++ b/xcube/webapi/ows/coverages/routes.py
@@ -0,0 +1,205 @@
+# 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 json
+import re
+from typing import Collection, Optional
+import fnmatch
+from xcube.server.api import ApiHandler, ApiRequest
+from .api import api
+from .context import CoveragesContext
+from .controllers import (
+    get_coverage_as_json,
+    get_coverage_domainset,
+    get_coverage_rangetype,
+    get_collection_metadata,
+    get_coverage_data,
+)
+
+
+# noinspection PyAbstractClass,PyMethodMayBeStatic
+@api.route('/catalog/collections/{collectionId}/coverage')
+class CoveragesCoverageHandler(ApiHandler[CoveragesContext]):
+    """
+    Return coverage data
+
+    This is the main coverage endpoint: the one that returns the actual
+    data (as opposed to metadata), as TIFF or NetCDF. It can also provide
+    representations in HTML and JSON.
+    """
+    # noinspection PyPep8Naming
+    @api.operation(
+        operation_id='coveragesCoverage',
+        summary='A coverage in OGC API - Coverages',
+    )
+    async def get(self, collectionId: str):
+        ds_ctx = self.ctx.datasets_ctx
+        available_types = [
+            'image/tiff',
+            'application/x-geotiff',
+            'text/html',
+            'application/json',
+            'application/netcdf',
+            'application/x-netcdf',
+        ]
+        content_type = negotiate_content_type(self.request, available_types)
+        if content_type is None:
+            self.response.set_status(415)
+            return await self.response.finish(
+                f'Available media types: {", ".join(available_types)}\n'
+            )
+        elif content_type == 'text/html':
+            result = (
+                '<html><title>Collection</title><body><pre>\n' +
+                json.dumps(get_coverage_as_json(ds_ctx, collectionId), indent=2)
+                + '\n</pre></body></html>'
+            )
+        elif content_type == 'application/json':
+            result = get_coverage_as_json(ds_ctx, collectionId)
+        else:
+            result = get_coverage_data(
+                ds_ctx, collectionId, self.request.query, content_type
+            )
+        return await self.response.finish(result, content_type=content_type)
+
+
+# noinspection PyAbstractClass,PyMethodMayBeStatic
+@api.route('/catalog/collections/{collectionId}/coverage/domainset')
+class CoveragesDomainsetHandler(ApiHandler[CoveragesContext]):
+    """
+    Describe the domain set of a coverage
+
+    The domain set is the set of input parameters (e.g. geographical extent,
+    time span) for which the coverage is defined.
+    """
+    # noinspection PyPep8Naming
+    @api.operation(
+        operation_id='coveragesDomainSet',
+        summary='OGC API - Coverages - domain set',
+    )
+    async def get(self, collectionId: str):
+        domain_set = get_coverage_domainset(
+            self.ctx.datasets_ctx, collectionId
+        )
+        return self.response.finish(domain_set)
+
+
+# noinspection PyAbstractClass,PyMethodMayBeStatic
+@api.route('/catalog/collections/{collectionId}/coverage/rangetype')
+class CoveragesRangetypeHandler(ApiHandler[CoveragesContext]):
+    """
+    Describe the range type of a coverage
+
+    The range type describes the types of the data contained in the range
+    of this coverage. For a data cube, these would typically correspond to
+    the types of the variables or bands.
+    """
+    # noinspection PyPep8Naming
+    @api.operation(
+        operation_id='coveragesRangeType',
+        summary='OGC API - Coverages - range type',
+    )
+    async def get(self, collectionId: str):
+        range_type = get_coverage_rangetype(
+            self.ctx.datasets_ctx, collectionId
+        )
+        return self.response.finish(range_type)
+
+
+@api.route('/catalog/collections/{collectionId}/coverage/metadata')
+class CoveragesMetadataHandler(ApiHandler[CoveragesContext]):
+    """
+    Return coverage metadata
+
+    The metadata is taken from the source dataset's attributes
+    """
+    # noinspection PyPep8Naming
+    @api.operation(
+        operation_id='coveragesMetadata',
+        summary='OGC API - Coverages - metadata',
+    )
+    async def get(self, collectionId: str):
+        return self.response.finish(
+            get_collection_metadata(self.ctx.datasets_ctx, collectionId)
+        )
+
+
+@api.route('/catalog/collections/{collectionId}/coverage/rangeset')
+class CoveragesRangesetHandler(ApiHandler[CoveragesContext]):
+    """
+    Handle rangeset endpoint with a "not supported" response
+
+    This endpoint has been deprecated
+    (see https://github.com/m-mohr/geodatacube-api/pull/1 ),
+    but we handle it to make clear that its non-implementation is a deliberate
+    decision.
+    """
+    # noinspection PyPep8Naming
+    @api.operation(
+        operation_id='coveragesRangeset',
+        summary='OGC API - Coverages - rangeset',
+    )
+    async def get(self, collectionId: str):
+        self.response.set_status(501)
+        return self.response.finish(
+            'The rangeset endpoint has been deprecated and is not supported.'
+        )
+
+
+def negotiate_content_type(
+    request: ApiRequest, available: Collection[str]
+) -> Optional[str]:
+    """
+    Determine a MIME content type based on client and server capabilities
+
+    Client preferences are determined both from the standard HTTP Accept
+    header provided by the client, and by an optional `f` parameter which
+    can override the Accept header.
+
+    :param request: HTTP request with Accept header and/or f parameter
+    :param available: List of MIME types that the server can produce
+    :return: the MIME type that is most acceptable to both client and server,
+             or None if there is no MIME type acceptable to both
+    """
+    if 'f' in request.query:  # overrides headers
+        content_type = request.query['f'][0]
+        return (
+            content_type
+            if available is None or content_type in available
+            else None
+        )
+
+    accept = re.split(', *', request.headers.get('Accept'))
+
+    def parse_part(part: str) -> tuple[float, str]:
+        if ';q=' in part:
+            subparts = part.split(';q=')
+            return float(subparts[1]), subparts[0]
+        else:
+            return 1, part
+
+    type_specs = sorted([parse_part(part) for part in accept], reverse=True)
+    types = [ts[1] for ts in type_specs]
+    for allowed_type in types:
+        for available_type in available:
+            # We (ab)use fnmatch to match * wildcards from accept headers
+            if fnmatch.fnmatch(available_type, allowed_type):
+                return available_type
+    return None


Test patch:
diff --git a/test/webapi/ows/coverages/__init__.py b/test/webapi/ows/coverages/__init__.py
new file mode 100644
index 00000000..454a3bf5
--- /dev/null
+++ b/test/webapi/ows/coverages/__init__.py
@@ -0,0 +1,20 @@
+# 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.
diff --git a/test/webapi/ows/coverages/expected.json b/test/webapi/ows/coverages/expected.json
new file mode 100644
index 00000000..272c98ae
--- /dev/null
+++ b/test/webapi/ows/coverages/expected.json
@@ -0,0 +1,188 @@
+{
+  "id": "demo",
+  "type": "CoverageByDomainAndRange",
+  "envelope": {
+    "type": "EnvelopeByAxis",
+    "srsName": "EPSG:4326",
+    "axisLabels": [
+      "lat",
+      "bnds",
+      "lon",
+      "time"
+    ],
+    "axis": [
+      {
+        "type": "RegularAxis",
+        "axisLabel": "lat",
+        "lowerBound": 52.49875,
+        "upperBound": 50.00125,
+        "resolution": 0.0024975000000000023,
+        "uomLabel": "degrees_north"
+      },
+      {
+        "type": "RegularAxis",
+        "axisLabel": "bnds",
+        "lowerBound": 0,
+        "upperBound": 1,
+        "resolution": 0.5,
+        "uomLabel": "degrees"
+      },
+      {
+        "type": "RegularAxis",
+        "axisLabel": "lon",
+        "lowerBound": 0.00125,
+        "upperBound": 4.99875,
+        "resolution": 0.00249875,
+        "uomLabel": "degrees_east"
+      },
+      {
+        "type": "RegularAxis",
+        "axisLabel": "time",
+        "lowerBound": 1484561361834255872,
+        "upperBound": 1485773193836892416,
+        "resolution": 242366400527308.8,
+        "uomLabel": "degrees"
+      }
+    ]
+  },
+  "domainSet": {
+    "type": "DomainSet",
+    "generalGrid": {
+      "type": "GeneralGridCoverage",
+      "srsName": "EPSG:4326",
+      "axisLabels": [
+        "lat",
+        "bnds",
+        "lon",
+        "time"
+      ],
+      "axis": [
+        {
+          "type": "RegularAxis",
+          "axisLabel": "lat",
+          "lowerBound": 52.49875,
+          "upperBound": 50.00125,
+          "resolution": 0.0024975000000000023,
+          "uomLabel": "degrees_north"
+        },
+        {
+          "type": "RegularAxis",
+          "axisLabel": "bnds",
+          "lowerBound": 0,
+          "upperBound": 1,
+          "resolution": 0.5,
+          "uomLabel": "degrees"
+        },
+        {
+          "type": "RegularAxis",
+          "axisLabel": "lon",
+          "lowerBound": 0.00125,
+          "upperBound": 4.99875,
+          "resolution": 0.00249875,
+          "uomLabel": "degrees_east"
+        },
+        {
+          "type": "RegularAxis",
+          "axisLabel": "time",
+          "lowerBound": 1484561361834255872,
+          "upperBound": 1485773193836892416,
+          "resolution": 242366400527308.8,
+          "uomLabel": "degrees"
+        }
+      ],
+      "gridLimits": {
+        "type": "GridLimits",
+        "srsName": "http://www.opengis.net/def/crs/OGC/0/Index4D",
+        "axisLabels": [
+          "lat",
+          "bnds",
+          "lon",
+          "time"
+        ],
+        "axis": [
+          {
+            "type": "IndexAxis",
+            "axisLabel": "lat",
+            "lowerBound": 0,
+            "upperBound": 1000
+          },
+          {
+            "type": "IndexAxis",
+            "axisLabel": "bnds",
+            "lowerBound": 0,
+            "upperBound": 2
+          },
+          {
+            "type": "IndexAxis",
+            "axisLabel": "lon",
+            "lowerBound": 0,
+            "upperBound": 2000
+          },
+          {
+            "type": "IndexAxis",
+            "axisLabel": "time",
+            "lowerBound": 0,
+            "upperBound": 5
+          }
+        ]
+      }
+    }
+  },
+  "rangeSet": {
+    "type": "RangeSet",
+    "dataBlock": {
+      "type": "VDataBlock",
+      "values": [
+        "TODO"
+      ]
+    }
+  },
+  "rangeType": {
+    "type": "DataRecord",
+    "field": [
+      {
+        "type": "Quantity",
+        "name": "c2rcc_flags",
+        "description": "C2RCC quality flags",
+        "encodingInfo": {
+          "dataType": "http://www.opengis.net/def/dataType/OGC/0/float64"
+        }
+      },
+      {
+        "type": "Quantity",
+        "name": "conc_chl",
+        "description": "Chlorophyll concentration",
+        "encodingInfo": {
+          "dataType": "http://www.opengis.net/def/dataType/OGC/0/float64"
+        }
+      },
+      {
+        "type": "Quantity",
+        "name": "conc_tsm",
+        "description": "Total suspended matter dry weight concentration",
+        "encodingInfo": {
+          "dataType": "http://www.opengis.net/def/dataType/OGC/0/float64"
+        }
+      },
+      {
+        "type": "Quantity",
+        "name": "kd489",
+        "description": "Irradiance attenuation coefficient at 489 nm",
+        "encodingInfo": {
+          "dataType": "http://www.opengis.net/def/dataType/OGC/0/float64"
+        }
+      },
+      {
+        "type": "Quantity",
+        "name": "quality_flags",
+        "description": "Classification and quality flags",
+        "encodingInfo": {
+          "dataType": "http://www.opengis.net/def/dataType/OGC/0/float64"
+        }
+      }
+    ]
+  },
+  "metadata": {
+    "Conventions": "CF-1.7"
+  }
+}
\ No newline at end of file
diff --git a/test/webapi/ows/coverages/test_context.py b/test/webapi/ows/coverages/test_context.py
new file mode 100644
index 00000000..64d550c6
--- /dev/null
+++ b/test/webapi/ows/coverages/test_context.py
@@ -0,0 +1,42 @@
+# 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 unittest
+from typing import Optional, Union, Mapping, Any
+
+from test.webapi.helpers import get_api_ctx
+from xcube.server.api import Context
+from xcube.webapi.datasets.context import DatasetsContext
+from xcube.webapi.ows.coverages.context import CoveragesContext
+
+
+def get_coverages_ctx(
+    server_config: Optional[Union[str, Mapping[str, Any]]] = None
+) -> CoveragesContext:
+    return get_api_ctx("ows.coverages", CoveragesContext, server_config)
+
+
+class CoveragesContextTest(unittest.TestCase):
+    def test_ctx_ok(self):
+        ctx = get_coverages_ctx()
+        self.assertIsInstance(ctx.server_ctx, Context)
+        self.assertIsInstance(ctx.datasets_ctx, DatasetsContext)
diff --git a/test/webapi/ows/coverages/test_controllers.py b/test/webapi/ows/coverages/test_controllers.py
new file mode 100644
index 00000000..59728861
--- /dev/null
+++ b/test/webapi/ows/coverages/test_controllers.py
@@ -0,0 +1,126 @@
+# 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 json
+import os
+import tempfile
+import unittest
+from io import BytesIO
+from pathlib import Path
+import numpy as np
+import xarray as xr
+import rioxarray
+
+from test.webapi.ows.coverages.test_context import get_coverages_ctx
+from xcube.webapi.ows.coverages.controllers import (
+    get_coverage_as_json,
+    get_coverage_data,
+    get_crs_from_dataset,
+    dtype_to_opengis_datatype,
+    get_dataarray_description,
+)
+
+
+class CoveragesControllersTest(unittest.TestCase):
+    def test_get_coverage_as_json(self):
+        result = get_coverage_as_json(get_coverages_ctx().datasets_ctx, 'demo')
+        self.assertIsInstance(result, dict)
+        path = Path(__file__).parent / 'expected.json'
+        # with open(path, mode="w") as fp:
+        #    json.dump(result, fp, indent=2)
+        with open(path, mode="r") as fp:
+            expected_result = json.load(fp)
+        self.assertEqual(expected_result, result)
+
+    def test_get_coverage_data_tiff(self):
+        query = dict(
+            bbox=['52,1,51,2'],
+            datetime=['2017-01-25'],
+            properties=['conc_chl'],
+        )
+        result = get_coverage_data(
+            get_coverages_ctx().datasets_ctx, 'demo', query, 'image/tiff'
+        )
+        with BytesIO(result) as fh:
+            da = rioxarray.open_rasterio(fh)
+            self.assertIsInstance(da, xr.DataArray)
+            self.assertEqual(('band', 'y', 'x'), da.dims)
+            self.assertEqual('Chlorophyll concentration', da.long_name)
+            self.assertEqual((1, 400, 400), da.shape)
+
+    def test_get_coverage_data_netcdf(self):
+        query = dict(
+            bbox=['52,1,51,2'],
+            datetime=['2017-01-25'],
+            properties=['conc_chl,kd489'],
+        )
+        result = get_coverage_data(
+            get_coverages_ctx().datasets_ctx,
+            'demo',
+            query,
+            'application/netcdf',
+        )
+        # We can't read this directly from memory: the netcdf4 engine only
+        # reads from filesystem paths, the h5netcdf engine (which can read
+        # from memory) isn't an xcube dependency, and the scipy engine only
+        # handles NetCDF 3.
+
+        with tempfile.TemporaryDirectory() as tempdir:
+            path = os.path.join(tempdir, 'out.nc')
+            with open(path, 'wb') as fh:
+                fh.write(result)
+            ds = xr.open_dataset(path)
+            self.assertEqual({'lat': 400, 'lon': 400, 'bnds': 2}, ds.dims)
+            self.assertEqual(['conc_chl', 'kd489'], list(ds.data_vars))
+            self.assertEqual(
+                [
+                    'lat',
+                    'lat_bnds',
+                    'lon',
+                    'lon_bnds',
+                    'time',
+                    'time_bnds',
+                    'conc_chl',
+                    'kd489',
+                ],
+                list(ds.variables),
+            )
+
+    def test_get_crs_from_dataset(self):
+        ds = xr.Dataset({'crs': ([], None, {'spatial_ref': '3035'})})
+        self.assertEqual('EPSG:3035', get_crs_from_dataset(ds))
+
+    def test_dtype_to_opengis_datatype(self):
+        expected = [
+            (
+                np.uint16,
+                'http://www.opengis.net/def/dataType/OGC/0/unsignedShort',
+            ),
+            (np.int32, 'http://www.opengis.net/def/dataType/OGC/0/signedInt'),
+            (np.datetime64, 'http://www.opengis.net/def/bipm/UTC'),
+            (np.object_, ''),
+        ]
+        for dtype, opengis in expected:
+            self.assertEqual(opengis, dtype_to_opengis_datatype(dtype))
+
+    def test_get_dataarray_description(self):
+        name = 'foo'
+        da = xr.DataArray(data=[], coords=[('x', [])], dims=['x'], name=name)
+        self.assertEqual(name, get_dataarray_description(da))
diff --git a/test/webapi/ows/coverages/test_routes.py b/test/webapi/ows/coverages/test_routes.py
new file mode 100644
index 00000000..bb573a1d
--- /dev/null
+++ b/test/webapi/ows/coverages/test_routes.py
@@ -0,0 +1,63 @@
+# 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.
+
+from ...helpers import RoutesTestCase
+
+
+class CoveragesRoutesTest(RoutesTestCase):
+    def test_fetch_coverage(self):
+        response = self.fetch(
+            '/catalog/collections/demo/coverage?f=application/json'
+        )
+        self.assertResponseOK(response)
+
+    def test_fetch_coverage_html(self):
+        response = self.fetch(
+            '/catalog/collections/demo/coverage',
+            headers=dict(
+                Accept='text/nonexistent,application/json;q=0.9,text/html;q=1.0'
+            ),
+        )
+        self.assertResponseOK(response)
+        self.assertEqual('text/html', response.headers['Content-Type'])
+
+    def test_fetch_coverage_wrong_media_type(self):
+        response = self.fetch(
+            '/catalog/collections/demo/coverage',
+            headers=dict(Accept='text/nonexistent'),
+        )
+        self.assertEqual(response.status, 415)
+
+    def test_fetch_domainset(self):
+        response = self.fetch('/catalog/collections/demo/coverage/domainset')
+        self.assertResponseOK(response)
+
+    def test_fetch_rangetype(self):
+        response = self.fetch('/catalog/collections/demo/coverage/rangetype')
+        self.assertResponseOK(response)
+
+    def test_fetch_metadata(self):
+        response = self.fetch('/catalog/collections/demo/coverage/metadata')
+        self.assertResponseOK(response)
+
+    def test_fetch_rangeset(self):
+        response = self.fetch('/catalog/collections/demo/coverage/rangeset')
+        self.assertEqual(response.status, 501)