Repository: dcs4cop/xcube

Issue:
Implement OGC Coverages scaling requirements class
The requirements class is defined [here](https://docs.ogc.org/DRAFTS/19-087.html#rc-scaling). Currently we have a partial implementation, supporting only the parameter `scale-factor`. 

Discussion:

Patch:
diff --git a/CHANGES.md b/CHANGES.md
index 5863576e..2b5482fa 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -6,7 +6,7 @@
   * Provide links for multiple coverages data formats
   * Add `crs` and `crs_storage` properties to STAC data
 * OGC Coverages:
-  * Support `scale-factor` parameter
+  * Support scaling parameters `scale-factor`, `scale-axes`, and `scale-size`
   * Improve handling of bbox parameters
   * Handle half-open datetime intervals
   * More robust and standard-compliant parameter parsing and checking
diff --git a/xcube/webapi/ows/coverages/controllers.py b/xcube/webapi/ows/coverages/controllers.py
index 1588f2d6..940ba9bb 100644
--- a/xcube/webapi/ows/coverages/controllers.py
+++ b/xcube/webapi/ows/coverages/controllers.py
@@ -25,6 +25,7 @@ from typing import Mapping, Sequence, Optional, Any, Literal, NamedTuple, Union
 
 import numpy as np
 import pyproj
+import rasterio
 import xarray as xr
 
 from xcube.core.gridmapping import GridMapping
@@ -32,7 +33,9 @@ from xcube.core.resampling import resample_in_space
 from xcube.server.api import ApiError
 from xcube.util.timeindex import ensure_time_index_compatible
 from xcube.webapi.datasets.context import DatasetsContext
-from xcube.webapi.ows.coverages.request import CoveragesRequest
+from xcube.webapi.ows.coverages.request import CoverageRequest
+from xcube.webapi.ows.coverages.scaling import CoverageScaling
+from xcube.webapi.ows.coverages.util import get_h_dim, get_v_dim
 
 
 def get_coverage_as_json(ctx: DatasetsContext, collection_id: str):
@@ -89,7 +92,7 @@ def get_coverage_data(
     ds = get_dataset(ctx, collection_id)
 
     try:
-        request = CoveragesRequest(query)
+        request = CoverageRequest(query)
     except ValueError as e:
         raise ApiError.BadRequest(str(e))
 
@@ -100,18 +103,7 @@ def get_coverage_data(
     subset_crs = request.subset_crs
 
     if request.properties is not None:
-        requested_vars = set(request.properties)
-        data_vars = set(map(str, ds.data_vars))
-        unrecognized_vars = requested_vars - data_vars
-        if unrecognized_vars == set():
-            ds = ds.drop_vars(
-                list(data_vars - requested_vars - {'crs', 'spatial_ref'})
-            )
-        else:
-            raise ApiError.BadRequest(
-                f'The following properties are not present in the coverage '
-                f'{collection_id}: {", ".join(unrecognized_vars)}'
-            )
+        ds = _apply_properties(collection_id, ds, request.properties)
 
     # https://docs.ogc.org/DRAFTS/19-087.html#datetime-parameter-subset-requirements
     # requirement 7D: "If a datetime parameter is specified requesting a
@@ -134,37 +126,38 @@ def get_coverage_data(
     if request.bbox is not None:
         ds = _apply_bbox(ds, request.bbox, bbox_crs, always_xy=False)
 
-    _assert_coverage_size_ok(ds, request.scale_factor)
-
-    source_gm = GridMapping.from_dataset(ds, crs=native_crs)
-    target_gm = None
+    # Do a provisional size check with an approximate scaling before attempting
+    # to determine a grid mapping, so the client gets a comprehensible error
+    # if the coverage is empty or too large.
+    _assert_coverage_size_ok(CoverageScaling(request, final_crs, ds))
 
+    transformed_gm = source_gm = GridMapping.from_dataset(ds, crs=native_crs)
     if native_crs != final_crs:
-        target_gm = source_gm.transform(final_crs).to_regular()
-    if request.scale_factor != 1:
-        if target_gm is None:
-            target_gm = source_gm
-        target_gm = target_gm.scale(1. / request.scale_factor)
-    if request.scale_axes is not None:
-        # TODO implement scale-axes
-        raise ApiError.NotImplemented(
-            'The scale-axes parameter is not yet supported.'
-        )
-    if request.scale_size:
-        # TODO implement scale-size
-        raise ApiError.NotImplemented(
-            'The scale-size parameter is not yet supported.'
+        transformed_gm = transformed_gm.transform(final_crs).to_regular()
+    if transformed_gm is not source_gm:
+        # We can't combine the scaling operation with this CRS transformation,
+        # since the size may end up wrong after the re-application of bounding
+        # boxes below.
+        ds = resample_in_space(
+            ds, source_gm=source_gm, target_gm=transformed_gm
         )
 
-    if target_gm is not None:
-        ds = resample_in_space(ds, source_gm=source_gm, target_gm=target_gm)
-
-    # In this case, the transformed native CRS bbox[es] may have been
-    # too big, so we re-crop in the final CRS.
-    if native_crs != final_crs and request.bbox is not None:
-        ds = _apply_bbox(ds, request.bbox, bbox_crs, always_xy=False)
-    if subset_bbox is not None:
-        ds = _apply_bbox(ds, subset_bbox, subset_crs, always_xy=True)
+    if native_crs != final_crs:
+        # If we've resampled into a new CRS, the transformed native-CRS
+        # bbox[es] may have been too big, so we re-crop in the final CRS.
+        if request.bbox is not None:
+            ds = _apply_bbox(ds, request.bbox, bbox_crs, always_xy=False)
+        if subset_bbox is not None:
+            ds = _apply_bbox(ds, subset_bbox, subset_crs, always_xy=True)
+
+    # Apply final size check and scaling operation after bbox, subsetting,
+    # and CRS transformation, to make sure that the final size is correct.
+    scaling = CoverageScaling(request, final_crs, ds)
+    _assert_coverage_size_ok(scaling)
+    if scaling.factor != (1, 1):
+        cropped_gm = GridMapping.from_dataset(ds, crs=final_crs)
+        scaled_gm = scaling.apply(cropped_gm)
+        ds = resample_in_space(ds, source_gm=cropped_gm, target_gm=scaled_gm)
 
     ds.rio.write_crs(final_crs, inplace=True)
     for var in ds.data_vars.values():
@@ -178,9 +171,9 @@ def get_coverage_data(
         netcdf={'netcdf', 'application/netcdf', 'application/x-netcdf'},
     )
     if content_type in media_types['tiff']:
-        content = dataset_to_image(ds, 'tiff')
+        content = dataset_to_image(ds, 'tiff', final_crs)
     elif content_type in media_types['png']:
-        content = dataset_to_image(ds, 'png')
+        content = dataset_to_image(ds, 'png', final_crs)
     elif content_type in media_types['netcdf']:
         content = dataset_to_netcdf(ds)
     else:
@@ -200,25 +193,28 @@ def get_coverage_data(
     return content, final_bbox, final_crs
 
 
-def _assert_coverage_size_ok(ds: xr.Dataset, scale_factor: float):
+def _apply_properties(collection_id, ds, properties):
+    requested_vars = set(properties)
+    data_vars = set(map(str, ds.data_vars))
+    unrecognized_vars = requested_vars - data_vars
+    if unrecognized_vars == set():
+        ds = ds.drop_vars(
+            list(data_vars - requested_vars - {'crs', 'spatial_ref'})
+        )
+    else:
+        raise ApiError.BadRequest(
+            f'The following properties are not present in the coverage '
+            f'{collection_id}: {", ".join(unrecognized_vars)}'
+        )
+    return ds
+
+
+def _assert_coverage_size_ok(scaling: CoverageScaling):
     size_limit = 4000 * 4000  # TODO make this configurable
-    h_dim = _get_h_dim(ds)
-    v_dim = _get_v_dim(ds)
-    for d in h_dim, v_dim:
-        size = ds.dims[d]
-        if size == 0:
-            # Requirement 8C currently specifies a 204 rather than 404 here,
-            # but spec will soon be updated to allow 404 as an alternative.
-            # (J. Jacovella-St-Louis, pers. comm., 2023-11-27).
-            raise ApiError.NotFound(
-                f'Requested coverage contains no data: {d} has zero size.'
-            )
-    if (h_size := ds.dims[h_dim] / scale_factor) * (
-        y_size := ds.dims[v_dim] / scale_factor
-    ) > size_limit:
+    x, y = scaling.size
+    if (x * y) > size_limit:
         raise ApiError.ContentTooLarge(
-            f'Requested coverage is too large:'
-            f'{h_size} × {y_size} > {size_limit}.'
+            f'Requested coverage is too large:' f'{x} × {y} > {size_limit}.'
         )
 
 
@@ -308,7 +304,7 @@ def _apply_geographic_subsetting(
     # subsetting is only specified in one dimension.
     full_bbox_native = get_bbox_from_ds(ds)
     native_crs = get_crs_from_dataset(ds)
-    full_bbox_subset_crs = _transform_bbox(
+    full_bbox_subset_crs = transform_bbox(
         full_bbox_native, native_crs, subset_crs
     )
 
@@ -330,12 +326,12 @@ def _apply_geographic_subsetting(
     bbox_subset_crs = [x0, y0, x1, y1]
 
     # 4. Transform requested bbox from subsetting CRS to dataset-native CRS.
-    bbox_native_crs = _transform_bbox(bbox_subset_crs, subset_crs, native_crs)
+    bbox_native_crs = transform_bbox(bbox_subset_crs, subset_crs, native_crs)
 
     # 6. Apply the dataset-native bbox using sel, making sure that y/latitude
     # slice has the same ordering as the corresponding co-ordinate.
-    h_dim = _get_h_dim(ds)
-    v_dim = _get_v_dim(ds)
+    h_dim = get_h_dim(ds)
+    v_dim = get_v_dim(ds)
     ds = ds.sel(
         indexers={
             h_dim: slice(bbox_native_crs[0], bbox_native_crs[2]),
@@ -350,7 +346,7 @@ def _apply_geographic_subsetting(
 
 
 def get_bbox_from_ds(ds: xr.Dataset):
-    h, v = ds[_get_h_dim(ds)], ds[_get_v_dim(ds)]
+    h, v = ds[get_h_dim(ds)], ds[get_v_dim(ds)]
     bbox = list(map(float, [h[0], v[0], h[-1], v[-1]]))
     _ensure_bbox_y_ascending(bbox)
     return bbox
@@ -368,7 +364,7 @@ def _find_geographic_parameters(
     return x, y
 
 
-def _transform_bbox(
+def transform_bbox(
     bbox: list[float], source_crs: pyproj.CRS, dest_crs: pyproj.CRS
 ) -> list[float]:
     if source_crs == dest_crs:
@@ -392,8 +388,8 @@ def _apply_bbox(
         )
         _ensure_bbox_y_ascending(bbox, always_xy or is_xy_order(bbox_crs))
         bbox = transformer.transform_bounds(*bbox)
-    h_dim = _get_h_dim(ds)
-    v_dim = _get_v_dim(ds)
+    h_dim = get_h_dim(ds)
+    v_dim = get_v_dim(ds)
     x0, y0, x1, y1 = (
         (0, 1, 2, 3)
         if (always_xy or is_xy_order(native_crs))
@@ -412,18 +408,6 @@ def _ensure_bbox_y_ascending(bbox: list, xy_order: bool = True):
         bbox[y0], bbox[y1] = bbox[y1], bbox[y0]
 
 
-def _get_h_dim(ds: xr.Dataset):
-    return [
-        d for d in list(map(str, ds.dims)) if d[:3].lower() in {'x', 'lon'}
-    ][0]
-
-
-def _get_v_dim(ds: xr.Dataset):
-    return [
-        d for d in list(map(str, ds.dims)) if d[:3].lower() in {'y', 'lat'}
-    ][0]
-
-
 def _correct_inverted_y_range_if_necessary(
     ds: xr.Dataset, axis: str, range_: tuple[float, float]
 ) -> tuple[float, float]:
@@ -440,13 +424,16 @@ def _correct_inverted_y_range_if_necessary(
 
 
 def dataset_to_image(
-    ds: xr.Dataset, image_format: Literal['png', 'tiff'] = 'png'
+    ds: xr.Dataset,
+    image_format: Literal['png', 'tiff'] = 'png',
+    crs: pyproj.CRS = None,
 ) -> bytes:
     """
     Return an in-memory bitmap (TIFF or PNG) representing a dataset
 
     :param ds: a dataset
     :param image_format: image format to generate ("png" or "tiff")
+    :param crs: CRS of the dataset
     :return: TIFF-formatted bytes representing the dataset
     """
 
@@ -469,6 +456,9 @@ def dataset_to_image(
             ds[list(ds.data_vars)[0]].rio.to_raster(path)
         else:
             ds.rio.to_raster(path)
+        if crs is not None:
+            with rasterio.open(path, mode='r+') as src:
+                src.crs = crs
         with open(path, 'rb') as fh:
             data = fh.read()
     return data
diff --git a/xcube/webapi/ows/coverages/request.py b/xcube/webapi/ows/coverages/request.py
index 7ef06c7d..0a921ddf 100644
--- a/xcube/webapi/ows/coverages/request.py
+++ b/xcube/webapi/ows/coverages/request.py
@@ -24,7 +24,7 @@ import pyproj
 import rfc3339_validator
 
 
-class CoveragesRequest:
+class CoverageRequest:
     """A representation of a parsed OGC API - Coverages request
 
     As defined in https://docs.ogc.org/DRAFTS/19-087.html
@@ -127,7 +127,9 @@ class CoveragesRequest:
             except ValueError:
                 raise ValueError(f'Invalid scale-factor {scale_factor_str}')
         else:
-            self.scale_factor = 1
+            # We don't default to 1, since (per the standard) an implementation
+            # may choose to downscale by default.
+            self.scale_factor = None
 
     def _parse_scale_axes(self):
         self.scale_axes = (
diff --git a/xcube/webapi/ows/coverages/scaling.py b/xcube/webapi/ows/coverages/scaling.py
new file mode 100644
index 00000000..d041fa3c
--- /dev/null
+++ b/xcube/webapi/ows/coverages/scaling.py
@@ -0,0 +1,197 @@
+# 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 typing import Optional
+
+import pyproj
+import xarray as xr
+
+from xcube.core.gridmapping import GridMapping
+from xcube.server.api import ApiError
+from xcube.webapi.ows.coverages.util import get_h_dim, get_v_dim
+from xcube.webapi.ows.coverages.request import CoverageRequest
+
+
+class CoverageScaling:
+    """Representation of a scaling applied to an OGC coverage
+
+    This class represents the scaling specified in an OGC coverage request.
+    It is instantiated using a `CoverageRequest` instance, and can apply
+    itself to a `GridMapping` instance.
+    """
+
+    _scale: Optional[tuple[float, float]] = None
+    _final_size: Optional[tuple[int, int]] = None
+    _initial_size: tuple[int, int]
+    _crs: pyproj.CRS
+    _x_name: str
+    _y_name: str
+
+    def __init__(
+        self, request: CoverageRequest, crs: pyproj.CRS, ds: xr.Dataset
+    ):
+        """Create a new scaling from a coverages request object
+
+        :param request: a request optionally including scaling parameters.
+          If any scaling parameters are present, the returned instance will
+          correspond to them. If no scaling parameters are present, a
+          default scaling will be used (currently 1:1)
+        :param crs: the CRS of the dataset to be scaled
+        :param ds: the dataset to be scaled
+        """
+        h_dim = get_h_dim(ds)
+        v_dim = get_v_dim(ds)
+        for d in h_dim, v_dim:
+            size = ds.dims[d]
+            if size == 0:
+                # Requirement 8C currently specifies a 204 rather than 404 here,
+                # but spec will soon be updated to allow 404 as an alternative.
+                # (J. Jacovella-St-Louis, pers. comm., 2023-11-27).
+                raise ApiError.NotFound(
+                    f'Requested coverage contains no data: {d} has zero size.'
+                )
+        self._initial_size = ds.dims[h_dim], ds.dims[v_dim]
+
+        self._crs = crs
+        self._y_name = self.get_axis_from_crs(
+            {'lat', 'latitude', 'geodetic latitude', 'n', 'north', 'y'}
+        )
+        self._x_name = self.get_axis_from_crs(
+            {'longitude', 'geodetic longitude', 'lon', 'e', 'east', 'x'}
+        )
+
+        # The standard doesn't define behaviour when multiple scale
+        # parameters are given. We choose to handle one and ignore the
+        # others in such cases.
+        if request.scale_factor is not None:
+            self._scale = request.scale_factor, request.scale_factor
+        elif request.scale_axes is not None:
+            self._scale = self._get_xy_values(request.scale_axes)
+        elif request.scale_size is not None:
+            # The standard allows floats for "scale-size" but mandates:
+            # "The returned coverage SHALL contain exactly the specified number
+            # of sample values". We can't return a fractional number of sample
+            # values, so truncate to int here.
+            x, y = self._get_xy_values(request.scale_size)
+            self._final_size = int(x), int(y)
+        else:
+            # The standard doesn't mandate this as a default; in future, we
+            # might choose to downscale automatically if a large coverage
+            # is requested without an explicit scaling parameter.
+            self._scale = (1, 1)
+
+    @property
+    def factor(self) -> tuple[float, float]:
+        """Return the two-dimensional scale factor of this scaling
+
+        The components of the scale tuple are expressed as downscaling
+        factors: values greater than 1 imply that the rescaled size
+        of the coverage in the corresponding dimension will be smaller than
+        the original size, and vice versa.
+
+        If the scaling was initially specified as a final size rather than
+        a factor, the factor property is an estimate based on the dataset
+        dimensions; the effective factor may be different when the scaling
+        is applied to a GridMapping.
+
+        :return: a 2-tuple of the x and y scale factors, in that order
+        """
+        if self._scale is not None:
+            return self._scale
+        else:
+            x_initial, y_initial = self._initial_size
+            x_final, y_final = self._final_size
+            return x_initial / x_final, y_initial / y_final
+
+    @property
+    def size(self) -> tuple[float, float]:
+        """Return the final coverage size produced by this scaling
+
+        :return: a 2-tuple of the scaled x and y sizes, in that order,
+          in pixels
+        """
+        if self._final_size is not None:
+            return self._final_size
+        else:
+            x_initial, y_initial = self._initial_size
+            x_scale, y_scale = self._scale
+            return x_initial / x_scale, y_initial / y_scale
+
+    def _get_xy_values(
+        self, axis_to_value: dict[str, float]
+    ) -> tuple[float, float]:
+        x, y = None, None
+        for axis in axis_to_value:
+            if axis.lower()[:3] in ['x', 'e', 'eas', 'lon', self._x_name]:
+                x = axis_to_value[axis]
+            if axis.lower()[:3] in ['y', 'n', 'nor', 'lat', self._y_name]:
+                y = axis_to_value[axis]
+        return x, y
+
+    def get_axis_from_crs(self, valid_identifiers: set[str]) -> Optional[str]:
+        """Find an axis abbreviation via a set of possible axis identifiers
+
+        This method operates on the CRS with which this scaling was
+        instantiated. It returns the abbreviation of the first axis in
+        the CRS which has either a name or abbreviation matching any string
+        in the supplied set.
+
+        :param valid_identifiers: a set of axis identifiers
+        :return: the abbreviation of the first axis in this scaling's CRS
+          whose name or abbreviation is in the supplied set, or `None` if
+          no such axis exists
+        """
+        for axis in self._crs.axis_info:
+            if not hasattr(axis, 'abbrev'):
+                continue
+            identifiers = set(
+                map(
+                    lambda attr: getattr(axis, attr, '').lower(),
+                    ['name', 'abbrev'],
+                )
+            )
+            if identifiers & valid_identifiers:
+                return axis.abbrev
+        return None
+
+    def apply(self, gm: GridMapping) -> GridMapping:
+        """Apply this scaling to a grid mapping
+
+        The supplied grid mapping is regularized before being scaled.
+
+        :param gm: a grid mapping to be scaled
+        :return: the supplied grid mapping, scaled according to this scaling.
+          If this scaling is 1:1, the returned grid mapping may be the
+          original object.
+        """
+        if self.factor == (1, 1):
+            return gm
+        else:
+            regular = gm.to_regular()
+            source = regular.size
+            # Even if the scaling was specified as a factor, we calculate
+            # from the (inferred) final size. If a factor was given,
+            # self.size is the final size as calculated from the originally
+            # specified dataset, which is what the client would expect. The
+            # regularized GridMapping might have a different size,
+            # so we don't want to apply a specified factor to it directly.
+            return regular.scale(
+                (self.size[0] / source[0], self.size[1] / source[1])
+            )
diff --git a/xcube/webapi/ows/coverages/util.py b/xcube/webapi/ows/coverages/util.py
new file mode 100644
index 00000000..72e41316
--- /dev/null
+++ b/xcube/webapi/ows/coverages/util.py
@@ -0,0 +1,35 @@
+# 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 xarray as xr
+
+
+def get_h_dim(ds: xr.Dataset) -> str:
+    return [
+        d for d in list(map(str, ds.dims)) if d[:3].lower() in {'x', 'lon'}
+    ][0]
+
+
+def get_v_dim(ds: xr.Dataset) -> str:
+    return [
+        d for d in list(map(str, ds.dims)) if d[:3].lower() in {'y', 'lat'}
+    ][0]


Test patch:
diff --git a/test/webapi/ows/coverages/test_controllers.py b/test/webapi/ows/coverages/test_controllers.py
index d1307d5f..35898056 100644
--- a/test/webapi/ows/coverages/test_controllers.py
+++ b/test/webapi/ows/coverages/test_controllers.py
@@ -39,6 +39,7 @@ from xcube.webapi.ows.coverages.controllers import (
     get_dataarray_description,
     get_units,
     is_xy_order,
+    transform_bbox,
 )
 
 
@@ -58,7 +59,25 @@ class CoveragesControllersTest(unittest.TestCase):
             'bbox': ['51,1,52,2'],
             'bbox-crs': ['[EPSG:4326]'],
             'datetime': ['2017-01-25T00:00:00Z'],
-            'properties': ['conc_chl']
+            'properties': ['conc_chl'],
+        }
+        content, content_bbox, content_crs = get_coverage_data(
+            get_coverages_ctx().datasets_ctx, 'demo', query, 'image/tiff'
+        )
+        with BytesIO(content) 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_geo_subset(self):
+        query = {
+            'subset': ['Lat(51:52),Lon(1:2)'],
+            'subset-crs': ['[EPSG:4326]'],
+            'datetime': ['2017-01-25T00:00:00Z'],
+            'properties': ['conc_chl'],
+            'crs': ['[OGC:CRS84]']
         }
         content, content_bbox, content_crs = get_coverage_data(
             get_coverages_ctx().datasets_ctx, 'demo', query, 'image/tiff'
@@ -72,11 +91,12 @@ class CoveragesControllersTest(unittest.TestCase):
 
     def test_get_coverage_data_netcdf(self):
         crs = 'OGC:CRS84'
+        # Unscaled size is 400, 400
         query = {
             'bbox': ['1,51,2,52'],
             'datetime': ['2017-01-24T00:00:00Z/2017-01-27T00:00:00Z'],
             'properties': ['conc_chl,kd489'],
-            'scale-factor': [2],
+            'scale-axes': ['lat(2),lon(2)'],
             'crs': [crs],
         }
         content, content_bbox, content_crs = get_coverage_data(
@@ -116,7 +136,7 @@ class CoveragesControllersTest(unittest.TestCase):
                     'time_bnds',
                     'conc_chl',
                     'kd489',
-                    'crs'
+                    'crs',
                 },
                 set(ds.variables),
             )
@@ -160,7 +180,7 @@ class CoveragesControllersTest(unittest.TestCase):
                     'time',
                     'time_bnds',
                     'conc_chl',
-                    'spatial_ref'
+                    'spatial_ref',
                 },
                 set(ds.variables),
             )
@@ -170,7 +190,7 @@ class CoveragesControllersTest(unittest.TestCase):
         query = {
             'subset': ['lat(52:51),lon(1:2),time(2017-01-25)'],
             'properties': ['conc_chl'],
-            'scale-factor': ['4'],
+            'scale-size': ['lat(100),lon(100)'],
         }
         content, content_bbox, content_crs = get_coverage_data(
             get_coverages_ctx().datasets_ctx, 'demo', query, 'png'
@@ -262,24 +282,6 @@ class CoveragesControllersTest(unittest.TestCase):
         ds = xr.Dataset({'crs': ([], None, {'spatial_ref': '3035'})})
         self.assertEqual('EPSG:3035', get_crs_from_dataset(ds).to_string())
 
-    def test_get_coverage_scale_axes(self):
-        with self.assertRaises(ApiError.NotImplemented):
-            get_coverage_data(
-                get_coverages_ctx().datasets_ctx,
-                'demo',
-                {'scale-axes': ['x(2)']},
-                'application/netcdf',
-            )
-
-    def test_get_coverage_scale_size(self):
-        with self.assertRaises(ApiError.NotImplemented):
-            get_coverage_data(
-                get_coverages_ctx().datasets_ctx,
-                'demo',
-                {'scale-size': ['x(2)']},
-                'application/netcdf',
-            )
-
     def test_dtype_to_opengis_datatype(self):
         expected = [
             (
@@ -316,3 +318,9 @@ class CoveragesControllersTest(unittest.TestCase):
                 )
             )
         )
+
+    def test_transform_bbox_same_crs(self):
+        self.assertEqual(
+            bbox := [1, 2, 3, 4],
+            transform_bbox(bbox, crs := pyproj.CRS('EPSG:4326'), crs),
+        )
diff --git a/test/webapi/ows/coverages/test_request.py b/test/webapi/ows/coverages/test_request.py
index 1c342058..bd391ab9 100644
--- a/test/webapi/ows/coverages/test_request.py
+++ b/test/webapi/ows/coverages/test_request.py
@@ -1,62 +1,62 @@
 import unittest
 import pyproj
 
-from xcube.webapi.ows.coverages.request import CoveragesRequest
+from xcube.webapi.ows.coverages.request import CoverageRequest
 
 
 class CoveragesRequestTest(unittest.TestCase):
     def test_parse_bbox(self):
-        self.assertIsNone(CoveragesRequest({}).bbox)
+        self.assertIsNone(CoverageRequest({}).bbox)
         self.assertEqual(
             [1.1, 2.2, 3.3, 4.4],
-            CoveragesRequest(dict(bbox=['1.1,2.2,3.3,4.4'])).bbox,
+            CoverageRequest(dict(bbox=['1.1,2.2,3.3,4.4'])).bbox,
         )
         with self.assertRaises(ValueError):
-            CoveragesRequest(dict(bbox=['foo,bar,baz']))
+            CoverageRequest(dict(bbox=['foo,bar,baz']))
         with self.assertRaises(ValueError):
-            CoveragesRequest(dict(bbox=['1.1,2.2,3.3']))
+            CoverageRequest(dict(bbox=['1.1,2.2,3.3']))
 
     def test_parse_bbox_crs(self):
         self.assertEqual(
             pyproj.CRS('OGC:CRS84'),
-            CoveragesRequest({}).bbox_crs,
+            CoverageRequest({}).bbox_crs,
         )
         self.assertEqual(
             pyproj.CRS(crs_spec := 'EPSG:4326'),
-            CoveragesRequest({'bbox-crs': [crs_spec]}).bbox_crs,
+            CoverageRequest({'bbox-crs': [crs_spec]}).bbox_crs,
         )
         self.assertEqual(
             pyproj.CRS(crs_spec := 'OGC:CRS84'),
-            CoveragesRequest({'bbox-crs': [f'[{crs_spec}]']}).bbox_crs,
+            CoverageRequest({'bbox-crs': [f'[{crs_spec}]']}).bbox_crs,
         )
         with self.assertRaises(ValueError):
-            CoveragesRequest({'bbox-crs': ['not a CRS specifier']})
+            CoverageRequest({'bbox-crs': ['not a CRS specifier']})
 
     def test_parse_datetime(self):
         dt0 = '2018-02-12T23:20:52Z'
         dt1 = '2019-02-12T23:20:52Z'
-        self.assertIsNone(CoveragesRequest({}).datetime)
-        self.assertEqual(dt0, CoveragesRequest({'datetime': [dt0]}).datetime)
+        self.assertIsNone(CoverageRequest({}).datetime)
+        self.assertEqual(dt0, CoverageRequest({'datetime': [dt0]}).datetime)
         self.assertEqual(
-            (dt0, None), CoveragesRequest({'datetime': [f'{dt0}/..']}).datetime
+            (dt0, None), CoverageRequest({'datetime': [f'{dt0}/..']}).datetime
         )
         self.assertEqual(
-            (None, dt1), CoveragesRequest({'datetime': [f'../{dt1}']}).datetime
+            (None, dt1), CoverageRequest({'datetime': [f'../{dt1}']}).datetime
         )
         self.assertEqual(
             (dt0, dt1),
-            CoveragesRequest({'datetime': [f'{dt0}/{dt1}']}).datetime,
+            CoverageRequest({'datetime': [f'{dt0}/{dt1}']}).datetime,
         )
         with self.assertRaises(ValueError):
-            CoveragesRequest({'datetime': [f'{dt0}/{dt0}/{dt1}']})
+            CoverageRequest({'datetime': [f'{dt0}/{dt0}/{dt1}']})
         with self.assertRaises(ValueError):
-            CoveragesRequest({'datetime': ['not a valid time string']})
+            CoverageRequest({'datetime': ['not a valid time string']})
 
     def test_parse_subset(self):
-        self.assertIsNone(CoveragesRequest({}).subset)
+        self.assertIsNone(CoverageRequest({}).subset)
         self.assertEqual(
             dict(Lat=('10', '20'), Lon=('30', None), time='2019-03-27'),
-            CoveragesRequest(
+            CoverageRequest(
                 dict(subset=['Lat(10:20),Lon(30:*),time("2019-03-27")'])
             ).subset,
         )
@@ -64,7 +64,7 @@ class CoveragesRequestTest(unittest.TestCase):
             dict(
                 Lat=(None, '20'), Lon='30', time=('2019-03-27', '2020-03-27')
             ),
-            CoveragesRequest(
+            CoverageRequest(
                 dict(
                     subset=[
                         'Lat(*:20),Lon(30),time("2019-03-27":"2020-03-27")'
@@ -73,62 +73,62 @@ class CoveragesRequestTest(unittest.TestCase):
             ).subset,
         )
         with self.assertRaises(ValueError):
-            CoveragesRequest({'subset': ['not a valid specifier']})
+            CoverageRequest({'subset': ['not a valid specifier']})
 
     def test_parse_subset_crs(self):
         self.assertEqual(
             pyproj.CRS('OGC:CRS84'),
-            CoveragesRequest({}).subset_crs,
+            CoverageRequest({}).subset_crs,
         )
         self.assertEqual(
             pyproj.CRS(crs_spec := 'EPSG:4326'),
-            CoveragesRequest({'subset-crs': [crs_spec]}).subset_crs,
+            CoverageRequest({'subset-crs': [crs_spec]}).subset_crs,
         )
         with self.assertRaises(ValueError):
-            CoveragesRequest({'subset-crs': ['not a CRS specifier']})
+            CoverageRequest({'subset-crs': ['not a CRS specifier']})
 
     def test_parse_properties(self):
-        self.assertIsNone(CoveragesRequest({}).properties)
+        self.assertIsNone(CoverageRequest({}).properties)
         self.assertEqual(
             ['foo', 'bar', 'baz'],
-            CoveragesRequest(dict(properties=['foo,bar,baz'])).properties,
+            CoverageRequest(dict(properties=['foo,bar,baz'])).properties,
         )
 
     def test_parse_scale_factor(self):
-        self.assertEqual(1, CoveragesRequest({}).scale_factor)
+        self.assertEqual(None, CoverageRequest({}).scale_factor)
         self.assertEqual(
-            1.5, CoveragesRequest({'scale-factor': ['1.5']}).scale_factor
+            1.5, CoverageRequest({'scale-factor': ['1.5']}).scale_factor
         )
         with self.assertRaises(ValueError):
-            CoveragesRequest({'scale-factor': ['this is not a number']})
+            CoverageRequest({'scale-factor': ['this is not a number']})
 
     def test_parse_scale_axes(self):
-        self.assertIsNone(CoveragesRequest({}).scale_axes)
+        self.assertIsNone(CoverageRequest({}).scale_axes)
         self.assertEqual(
             dict(Lat=1.5, Lon=2.5),
-            CoveragesRequest({'scale-axes': ['Lat(1.5),Lon(2.5)']}).scale_axes
+            CoverageRequest({'scale-axes': ['Lat(1.5),Lon(2.5)']}).scale_axes
         )
         with self.assertRaises(ValueError):
-            CoveragesRequest({'scale-axes': ['Lat(1.5']})
+            CoverageRequest({'scale-axes': ['Lat(1.5']})
         with self.assertRaises(ValueError):
-            CoveragesRequest({'scale-axes': ['Lat(not a number)']})
+            CoverageRequest({'scale-axes': ['Lat(not a number)']})
 
     def test_parse_scale_size(self):
-        self.assertIsNone(CoveragesRequest({}).scale_size)
+        self.assertIsNone(CoverageRequest({}).scale_size)
         self.assertEqual(
             dict(Lat=12.3, Lon=45.6),
-            CoveragesRequest({'scale-size': ['Lat(12.3),Lon(45.6)']}).scale_size
+            CoverageRequest({'scale-size': ['Lat(12.3),Lon(45.6)']}).scale_size
         )
         with self.assertRaises(ValueError):
-            CoveragesRequest({'scale-size': ['Lat(1.5']})
+            CoverageRequest({'scale-size': ['Lat(1.5']})
         with self.assertRaises(ValueError):
-            CoveragesRequest({'scale-size': ['Lat(not a number)']})
+            CoverageRequest({'scale-size': ['Lat(not a number)']})
 
     def test_parse_crs(self):
-        self.assertIsNone(CoveragesRequest({}).crs)
+        self.assertIsNone(CoverageRequest({}).crs)
         self.assertEqual(
             pyproj.CRS(crs := 'EPSG:4326'),
-            CoveragesRequest({'crs': [crs]}).crs
+            CoverageRequest({'crs': [crs]}).crs
         )
         with self.assertRaises(ValueError):
-            CoveragesRequest({'crs': ['an invalid CRS specifier']})
+            CoverageRequest({'crs': ['an invalid CRS specifier']})
diff --git a/test/webapi/ows/coverages/test_scaling.py b/test/webapi/ows/coverages/test_scaling.py
new file mode 100644
index 00000000..115ea59d
--- /dev/null
+++ b/test/webapi/ows/coverages/test_scaling.py
@@ -0,0 +1,75 @@
+import unittest
+from dataclasses import dataclass
+
+import pyproj
+
+import xcube.core.new
+from xcube.server.api import ApiError
+from xcube.webapi.ows.coverages.request import CoverageRequest
+from xcube.webapi.ows.coverages.scaling import CoverageScaling
+
+
+class ScalingTest(unittest.TestCase):
+    def setUp(self):
+        self.epsg4326 = pyproj.CRS('EPSG:4326')
+        self.ds = xcube.core.new.new_cube()
+
+    def test_default_scaling(self):
+        scaling = CoverageScaling(CoverageRequest({}), self.epsg4326, self.ds)
+        self.assertEqual((1, 1), scaling.factor)
+
+    def test_no_data(self):
+        with self.assertRaises(ApiError.NotFound):
+            CoverageScaling(
+                CoverageRequest({}),
+                self.epsg4326,
+                self.ds.isel(lat=slice(0, 0)),
+            )
+
+    def test_crs_no_valid_axis(self):
+        @dataclass
+        class CrsMock:
+            axis_info = [object()]
+
+        # noinspection PyTypeChecker
+        self.assertIsNone(
+            CoverageScaling(
+                CoverageRequest({}), CrsMock(), self.ds
+            ).get_axis_from_crs(set())
+        )
+
+    def test_scale_factor(self):
+        scaling = CoverageScaling(
+            CoverageRequest({'scale-factor': ['2']}), self.epsg4326, self.ds
+        )
+        self.assertEqual((2, 2), scaling.factor)
+        self.assertEqual((180, 90), scaling.size)
+
+    def test_scale_axes(self):
+        scaling = CoverageScaling(
+            CoverageRequest({'scale-axes': ['Lat(3),Lon(1.2)']}),
+            self.epsg4326,
+            self.ds,
+        )
+        self.assertEqual((1.2, 3), scaling.factor)
+        self.assertEqual((300, 60), scaling.size)
+
+    def test_scale_size(self):
+        scaling = CoverageScaling(
+            CoverageRequest({'scale-size': ['Lat(90),Lon(240)']}),
+            self.epsg4326,
+            self.ds,
+        )
+        self.assertEqual((240, 90), scaling.size)
+        self.assertEqual((1.5, 2), scaling.factor)
+
+    def test_apply_identity_scaling(self):
+        # noinspection PyTypeChecker
+        self.assertEqual(
+            gm_mock := object(),
+            CoverageScaling(
+                CoverageRequest({'scale-factor': ['1']}),
+                self.epsg4326,
+                self.ds,
+            ).apply(gm_mock),
+        )
diff --git a/test/webapi/ows/coverages/test_util.py b/test/webapi/ows/coverages/test_util.py
new file mode 100644
index 00000000..26043117
--- /dev/null
+++ b/test/webapi/ows/coverages/test_util.py
@@ -0,0 +1,18 @@
+import unittest
+import xcube.core.new
+import xcube.webapi.ows.coverages.util as util
+
+
+class UtilTest(unittest.TestCase):
+
+    def setUp(self):
+        self.ds_latlon = xcube.core.new.new_cube()
+        self.ds_xy = xcube.core.new.new_cube(x_name='x', y_name='y')
+
+    def test_get_h_dim(self):
+        self.assertEqual('lon', util.get_h_dim(self.ds_latlon))
+        self.assertEqual('x', util.get_h_dim(self.ds_xy))
+
+    def test_get_v_dim(self):
+        self.assertEqual('lat', util.get_v_dim(self.ds_latlon))
+        self.assertEqual('y', util.get_v_dim(self.ds_xy))