From 3204918e9a17cfe8f1749299ea8a5e6737932de2 Mon Sep 17 00:00:00 2001
From: Pontus Lurcock <pontus.lurcock@brockmann-consult.de>
Date: Thu, 7 Dec 2023 15:50:14 +0100
Subject: [PATCH 01/14] Coverages: don't parse absent scale-factor as 1

The OGC standard allows an implementation to downscale by default
if no scale-factor is given, so handling of absent scale-factor
parameter has been moved from the request parser to the controller
(though for the moment we are still defaulting to 1).
---
 test/webapi/ows/coverages/test_request.py |  2 +-
 xcube/webapi/ows/coverages/controllers.py | 10 +++++++---
 xcube/webapi/ows/coverages/request.py     |  4 +++-
 3 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/test/webapi/ows/coverages/test_request.py b/test/webapi/ows/coverages/test_request.py
index 1c342058e..9abe97803 100644
--- a/test/webapi/ows/coverages/test_request.py
+++ b/test/webapi/ows/coverages/test_request.py
@@ -95,7 +95,7 @@ def test_parse_properties(self):
         )
 
     def test_parse_scale_factor(self):
-        self.assertEqual(1, CoveragesRequest({}).scale_factor)
+        self.assertEqual(None, CoveragesRequest({}).scale_factor)
         self.assertEqual(
             1.5, CoveragesRequest({'scale-factor': ['1.5']}).scale_factor
         )
diff --git a/xcube/webapi/ows/coverages/controllers.py b/xcube/webapi/ows/coverages/controllers.py
index 1588f2d6a..04f43b052 100644
--- a/xcube/webapi/ows/coverages/controllers.py
+++ b/xcube/webapi/ows/coverages/controllers.py
@@ -134,17 +134,21 @@ 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)
+    # NB: a default scale factor of 1 is not compulsory. We may in future
+    # choose to downscale by default if a large coverage is requested.
+    scale_factor = 1 if request.scale_factor is None else request.scale_factor
+    _assert_coverage_size_ok(ds, scale_factor)
 
     source_gm = GridMapping.from_dataset(ds, crs=native_crs)
     target_gm = None
 
     if native_crs != final_crs:
         target_gm = source_gm.transform(final_crs).to_regular()
-    if request.scale_factor != 1:
+
+    if scale_factor != 1:
         if target_gm is None:
             target_gm = source_gm
-        target_gm = target_gm.scale(1. / request.scale_factor)
+        target_gm = target_gm.scale(1. / scale_factor)
     if request.scale_axes is not None:
         # TODO implement scale-axes
         raise ApiError.NotImplemented(
diff --git a/xcube/webapi/ows/coverages/request.py b/xcube/webapi/ows/coverages/request.py
index 7ef06c7d1..9a2241090 100644
--- a/xcube/webapi/ows/coverages/request.py
+++ b/xcube/webapi/ows/coverages/request.py
@@ -127,7 +127,9 @@ def _parse_scale_factor(self):
             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 = (

From d47f59e48902c99734a95810ec2cb2555fd8af5d Mon Sep 17 00:00:00 2001
From: Pontus Lurcock <pontus.lurcock@brockmann-consult.de>
Date: Fri, 8 Dec 2023 18:02:44 +0100
Subject: [PATCH 02/14] OGC Coverages: add CoverageScaling class

This is a helper class for implementing coverage scaling operations.
Not yet fully implemented and not yet in use.
---
 test/webapi/ows/coverages/test_scaling.py |  56 ++++++++++
 xcube/webapi/ows/coverages/controllers.py |  23 ++--
 xcube/webapi/ows/coverages/scaling.py     | 122 ++++++++++++++++++++++
 3 files changed, 192 insertions(+), 9 deletions(-)
 create mode 100644 test/webapi/ows/coverages/test_scaling.py
 create mode 100644 xcube/webapi/ows/coverages/scaling.py

diff --git a/test/webapi/ows/coverages/test_scaling.py b/test/webapi/ows/coverages/test_scaling.py
new file mode 100644
index 000000000..1cc6344dd
--- /dev/null
+++ b/test/webapi/ows/coverages/test_scaling.py
@@ -0,0 +1,56 @@
+import unittest
+
+import pyproj
+
+import xcube.core.new
+from xcube.server.api import ApiError
+from xcube.webapi.ows.coverages.request import CoveragesRequest
+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(CoveragesRequest({}), self.epsg4326, self.ds)
+        self.assertEqual((1, 1), scaling.scale)
+
+    def test_no_data(self):
+        with self.assertRaises(ApiError.NotFound):
+            CoverageScaling(
+                CoveragesRequest({}), self.epsg4326,
+                self.ds.isel(lat=slice(0, 0))
+            )
+
+    def test_crs_fallbacks(self):
+        # TODO: implement me
+        pass
+
+    def test_scale_factor(self):
+        scaling = CoverageScaling(
+            CoveragesRequest({'scale-factor': ['2']}),
+            self.epsg4326,
+            self.ds
+        )
+        self.assertEqual((2, 2), scaling.scale)
+
+    def test_scale_axes(self):
+        scaling = CoverageScaling(
+            CoveragesRequest({'scale-axes': ['Lat(3),Lon(1.2)']}),
+            self.epsg4326,
+            self.ds
+        )
+        self.assertEqual((1.2, 3), scaling.scale)
+        self.assertEqual((300, 60), scaling.size)
+
+    def test_scale_size(self):
+        scaling = CoverageScaling(
+            CoveragesRequest({'scale-size': ['Lat(90),Lon(240)']}),
+            self.epsg4326,
+            self.ds
+        )
+        self.assertEqual((240, 90), scaling.size)
+        self.assertEqual((1.5, 2), scaling.scale)
diff --git a/xcube/webapi/ows/coverages/controllers.py b/xcube/webapi/ows/coverages/controllers.py
index 04f43b052..34c98e5fd 100644
--- a/xcube/webapi/ows/coverages/controllers.py
+++ b/xcube/webapi/ows/coverages/controllers.py
@@ -206,8 +206,8 @@ def get_coverage_data(
 
 def _assert_coverage_size_ok(ds: xr.Dataset, scale_factor: float):
     size_limit = 4000 * 4000  # TODO make this configurable
-    h_dim = _get_h_dim(ds)
-    v_dim = _get_v_dim(ds)
+    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:
@@ -338,8 +338,8 @@ def _apply_geographic_subsetting(
 
     # 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]),
@@ -354,12 +354,17 @@ 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
 
 
+def apply_scaling(gm: GridMapping, ds: xr.Dataset, request: CoveragesRequest):
+    # TODO: implement me
+    pass
+
+
 def _find_geographic_parameters(
     names: list[str],
 ) -> tuple[Optional[str], Optional[str]]:
@@ -396,8 +401,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))
@@ -416,13 +421,13 @@ 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):
+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):
+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]
diff --git a/xcube/webapi/ows/coverages/scaling.py b/xcube/webapi/ows/coverages/scaling.py
new file mode 100644
index 000000000..8ac244c6a
--- /dev/null
+++ b/xcube/webapi/ows/coverages/scaling.py
@@ -0,0 +1,122 @@
+# 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.server.api import ApiError
+from xcube.webapi.ows.coverages.controllers import get_h_dim, get_v_dim
+from xcube.webapi.ows.coverages.request import CoveragesRequest
+
+
+class CoverageScaling:
+
+    _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: CoveragesRequest, crs: pyproj.CRS,
+                 ds: xr.Dataset):
+        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 scale(self) -> tuple[float, float]:
+        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]:
+        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]):
+        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

From 33073cfd8282e02d2d821946ca765cdeefd21878 Mon Sep 17 00:00:00 2001
From: Pontus Lurcock <pontus.lurcock@brockmann-consult.de>
Date: Mon, 11 Dec 2023 14:47:42 +0100
Subject: [PATCH 03/14] OGC Coverages: improve unit test coverage

---
 test/webapi/ows/coverages/test_scaling.py | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/test/webapi/ows/coverages/test_scaling.py b/test/webapi/ows/coverages/test_scaling.py
index 1cc6344dd..05de31528 100644
--- a/test/webapi/ows/coverages/test_scaling.py
+++ b/test/webapi/ows/coverages/test_scaling.py
@@ -1,4 +1,6 @@
+import collections
 import unittest
+from dataclasses import dataclass
 
 import pyproj
 
@@ -25,9 +27,15 @@ def test_no_data(self):
                 self.ds.isel(lat=slice(0, 0))
             )
 
-    def test_crs_fallbacks(self):
-        # TODO: implement me
-        pass
+    def test_crs_no_valid_axis(self):
+        @dataclass
+        class CrsMock:
+            axis_info = [object()]
+        # noinspection PyTypeChecker
+        self.assertIsNone(
+            CoverageScaling(CoveragesRequest({}), CrsMock(), self.ds)
+            .get_axis_from_crs(set())
+        )
 
     def test_scale_factor(self):
         scaling = CoverageScaling(

From 7c3a6094767d16f54538fed84ef1829bddf6b869 Mon Sep 17 00:00:00 2001
From: Pontus Lurcock <pontus.lurcock@brockmann-consult.de>
Date: Mon, 11 Dec 2023 16:11:57 +0100
Subject: [PATCH 04/14] OGC Coverages: add a unit test

---
 test/webapi/ows/coverages/test_controllers.py | 8 +++++++-
 test/webapi/ows/coverages/test_scaling.py     | 1 -
 xcube/webapi/ows/coverages/controllers.py     | 6 +++---
 3 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/test/webapi/ows/coverages/test_controllers.py b/test/webapi/ows/coverages/test_controllers.py
index d1307d5fa..c8770d43f 100644
--- a/test/webapi/ows/coverages/test_controllers.py
+++ b/test/webapi/ows/coverages/test_controllers.py
@@ -38,7 +38,7 @@
     dtype_to_opengis_datatype,
     get_dataarray_description,
     get_units,
-    is_xy_order,
+    is_xy_order, transform_bbox,
 )
 
 
@@ -316,3 +316,9 @@ def test_is_xy(self):
                 )
             )
         )
+
+    def test_transform_bbox_same_crs(self):
+        self.assertEqual(
+            bbox := [1, 2, 3, 4],
+            transform_bbox(bbox, crs := pyproj.CRS('EPSG:4326'), crs)
+        )
\ No newline at end of file
diff --git a/test/webapi/ows/coverages/test_scaling.py b/test/webapi/ows/coverages/test_scaling.py
index 05de31528..4d17d972b 100644
--- a/test/webapi/ows/coverages/test_scaling.py
+++ b/test/webapi/ows/coverages/test_scaling.py
@@ -1,4 +1,3 @@
-import collections
 import unittest
 from dataclasses import dataclass
 
diff --git a/xcube/webapi/ows/coverages/controllers.py b/xcube/webapi/ows/coverages/controllers.py
index 34c98e5fd..11f822ee4 100644
--- a/xcube/webapi/ows/coverages/controllers.py
+++ b/xcube/webapi/ows/coverages/controllers.py
@@ -312,7 +312,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
     )
 
@@ -334,7 +334,7 @@ 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.
@@ -377,7 +377,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:

From 843a6fb92262973a3c6a9034cdbadbfd771d7f60 Mon Sep 17 00:00:00 2001
From: Pontus Lurcock <pontus.lurcock@brockmann-consult.de>
Date: Mon, 11 Dec 2023 16:19:33 +0100
Subject: [PATCH 05/14] Rename OGC CoveragesRequest to CoverageRequest

---
 test/webapi/ows/coverages/test_request.py | 78 +++++++++++------------
 test/webapi/ows/coverages/test_scaling.py | 14 ++--
 xcube/webapi/ows/coverages/controllers.py |  6 +-
 xcube/webapi/ows/coverages/request.py     |  2 +-
 xcube/webapi/ows/coverages/scaling.py     |  4 +-
 5 files changed, 52 insertions(+), 52 deletions(-)

diff --git a/test/webapi/ows/coverages/test_request.py b/test/webapi/ows/coverages/test_request.py
index 9abe97803..bd391ab98 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 @@ def test_parse_subset(self):
             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 @@ def test_parse_subset(self):
             ).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(None, 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
index 4d17d972b..c01dff775 100644
--- a/test/webapi/ows/coverages/test_scaling.py
+++ b/test/webapi/ows/coverages/test_scaling.py
@@ -5,7 +5,7 @@
 
 import xcube.core.new
 from xcube.server.api import ApiError
-from xcube.webapi.ows.coverages.request import CoveragesRequest
+from xcube.webapi.ows.coverages.request import CoverageRequest
 from xcube.webapi.ows.coverages.scaling import CoverageScaling
 
 
@@ -16,13 +16,13 @@ def setUp(self):
         self.ds = xcube.core.new.new_cube()
 
     def test_default_scaling(self):
-        scaling = CoverageScaling(CoveragesRequest({}), self.epsg4326, self.ds)
+        scaling = CoverageScaling(CoverageRequest({}), self.epsg4326, self.ds)
         self.assertEqual((1, 1), scaling.scale)
 
     def test_no_data(self):
         with self.assertRaises(ApiError.NotFound):
             CoverageScaling(
-                CoveragesRequest({}), self.epsg4326,
+                CoverageRequest({}), self.epsg4326,
                 self.ds.isel(lat=slice(0, 0))
             )
 
@@ -32,13 +32,13 @@ class CrsMock:
             axis_info = [object()]
         # noinspection PyTypeChecker
         self.assertIsNone(
-            CoverageScaling(CoveragesRequest({}), CrsMock(), self.ds)
+            CoverageScaling(CoverageRequest({}), CrsMock(), self.ds)
             .get_axis_from_crs(set())
         )
 
     def test_scale_factor(self):
         scaling = CoverageScaling(
-            CoveragesRequest({'scale-factor': ['2']}),
+            CoverageRequest({'scale-factor': ['2']}),
             self.epsg4326,
             self.ds
         )
@@ -46,7 +46,7 @@ def test_scale_factor(self):
 
     def test_scale_axes(self):
         scaling = CoverageScaling(
-            CoveragesRequest({'scale-axes': ['Lat(3),Lon(1.2)']}),
+            CoverageRequest({'scale-axes': ['Lat(3),Lon(1.2)']}),
             self.epsg4326,
             self.ds
         )
@@ -55,7 +55,7 @@ def test_scale_axes(self):
 
     def test_scale_size(self):
         scaling = CoverageScaling(
-            CoveragesRequest({'scale-size': ['Lat(90),Lon(240)']}),
+            CoverageRequest({'scale-size': ['Lat(90),Lon(240)']}),
             self.epsg4326,
             self.ds
         )
diff --git a/xcube/webapi/ows/coverages/controllers.py b/xcube/webapi/ows/coverages/controllers.py
index 11f822ee4..2c37fc5d6 100644
--- a/xcube/webapi/ows/coverages/controllers.py
+++ b/xcube/webapi/ows/coverages/controllers.py
@@ -32,7 +32,7 @@
 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
 
 
 def get_coverage_as_json(ctx: DatasetsContext, collection_id: str):
@@ -89,7 +89,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))
 
@@ -360,7 +360,7 @@ def get_bbox_from_ds(ds: xr.Dataset):
     return bbox
 
 
-def apply_scaling(gm: GridMapping, ds: xr.Dataset, request: CoveragesRequest):
+def apply_scaling(gm: GridMapping, ds: xr.Dataset, request: CoverageRequest):
     # TODO: implement me
     pass
 
diff --git a/xcube/webapi/ows/coverages/request.py b/xcube/webapi/ows/coverages/request.py
index 9a2241090..0a921ddf0 100644
--- a/xcube/webapi/ows/coverages/request.py
+++ b/xcube/webapi/ows/coverages/request.py
@@ -24,7 +24,7 @@
 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
diff --git a/xcube/webapi/ows/coverages/scaling.py b/xcube/webapi/ows/coverages/scaling.py
index 8ac244c6a..9ee4c9043 100644
--- a/xcube/webapi/ows/coverages/scaling.py
+++ b/xcube/webapi/ows/coverages/scaling.py
@@ -25,7 +25,7 @@
 
 from xcube.server.api import ApiError
 from xcube.webapi.ows.coverages.controllers import get_h_dim, get_v_dim
-from xcube.webapi.ows.coverages.request import CoveragesRequest
+from xcube.webapi.ows.coverages.request import CoverageRequest
 
 
 class CoverageScaling:
@@ -37,7 +37,7 @@ class CoverageScaling:
     _x_name: str
     _y_name: str
 
-    def __init__(self, request: CoveragesRequest, crs: pyproj.CRS,
+    def __init__(self, request: CoverageRequest, crs: pyproj.CRS,
                  ds: xr.Dataset):
         h_dim = get_h_dim(ds)
         v_dim = get_v_dim(ds)

From 597f8ae844e4159246810505ccce34fff39415dc Mon Sep 17 00:00:00 2001
From: Pontus Lurcock <pontus.lurcock@brockmann-consult.de>
Date: Mon, 11 Dec 2023 17:20:28 +0100
Subject: [PATCH 06/14] OGC API - Coverages: add full support for scaling

scale-size and scale-axes are now supported in addition to the
previously supported scale-factor. The recently added CoverageScaling
class is used for the implementation.
---
 test/webapi/ows/coverages/test_controllers.py | 22 +----------
 xcube/webapi/ows/coverages/controllers.py     | 38 ++++---------------
 xcube/webapi/ows/coverages/scaling.py         |  9 ++++-
 xcube/webapi/ows/coverages/util.py            | 35 +++++++++++++++++
 4 files changed, 53 insertions(+), 51 deletions(-)
 create mode 100644 xcube/webapi/ows/coverages/util.py

diff --git a/test/webapi/ows/coverages/test_controllers.py b/test/webapi/ows/coverages/test_controllers.py
index c8770d43f..d93d36f95 100644
--- a/test/webapi/ows/coverages/test_controllers.py
+++ b/test/webapi/ows/coverages/test_controllers.py
@@ -76,7 +76,7 @@ def test_get_coverage_data_netcdf(self):
             '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(
@@ -170,7 +170,7 @@ def test_get_coverage_data_png(self):
         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 +262,6 @@ def test_get_crs_from_dataset(self):
         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 = [
             (
diff --git a/xcube/webapi/ows/coverages/controllers.py b/xcube/webapi/ows/coverages/controllers.py
index 2c37fc5d6..996b84d87 100644
--- a/xcube/webapi/ows/coverages/controllers.py
+++ b/xcube/webapi/ows/coverages/controllers.py
@@ -33,6 +33,8 @@
 from xcube.util.timeindex import ensure_time_index_compatible
 from xcube.webapi.datasets.context import DatasetsContext
 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):
@@ -140,27 +142,15 @@ def get_coverage_data(
     _assert_coverage_size_ok(ds, scale_factor)
 
     source_gm = GridMapping.from_dataset(ds, crs=native_crs)
-    target_gm = None
+    target_gm = source_gm
 
     if native_crs != final_crs:
-        target_gm = source_gm.transform(final_crs).to_regular()
-
-    if scale_factor != 1:
-        if target_gm is None:
-            target_gm = source_gm
-        target_gm = target_gm.scale(1. / 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.'
-        )
+        target_gm = target_gm.transform(final_crs).to_regular()
 
-    if target_gm is not None:
+    scaling = CoverageScaling(request, final_crs, ds)
+    target_gm = scaling.apply(target_gm)
+
+    if target_gm is not source_gm:
         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
@@ -421,18 +411,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]:
diff --git a/xcube/webapi/ows/coverages/scaling.py b/xcube/webapi/ows/coverages/scaling.py
index 9ee4c9043..45c568572 100644
--- a/xcube/webapi/ows/coverages/scaling.py
+++ b/xcube/webapi/ows/coverages/scaling.py
@@ -23,8 +23,9 @@
 import pyproj
 import xarray as xr
 
+from xcube.core.gridmapping import GridMapping
 from xcube.server.api import ApiError
-from xcube.webapi.ows.coverages.controllers import get_h_dim, get_v_dim
+from xcube.webapi.ows.coverages.util import get_h_dim, get_v_dim
 from xcube.webapi.ows.coverages.request import CoverageRequest
 
 
@@ -120,3 +121,9 @@ def get_axis_from_crs(self, valid_identifiers: set[str]):
             if identifiers & valid_identifiers:
                 return axis.abbrev
         return None
+
+    def apply(self, gm: GridMapping):
+        if self.scale == (1, 1):
+            return gm
+        else:
+            return gm.scale((1 / self.scale[0], 1 / self.scale[1]))
diff --git a/xcube/webapi/ows/coverages/util.py b/xcube/webapi/ows/coverages/util.py
new file mode 100644
index 000000000..0c23ee3a4
--- /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):
+    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]

From f2edb8bbab370776f1990189a56eede722f2ef34 Mon Sep 17 00:00:00 2001
From: Pontus Lurcock <pontus.lurcock@brockmann-consult.de>
Date: Tue, 12 Dec 2023 12:20:54 +0100
Subject: [PATCH 07/14] OGC Coverages: minor refactoring and reformatting

---
 test/webapi/ows/coverages/test_controllers.py | 13 ++--
 test/webapi/ows/coverages/test_scaling.py     | 20 +++---
 xcube/webapi/ows/coverages/controllers.py     | 69 +++++++------------
 xcube/webapi/ows/coverages/scaling.py         | 18 ++---
 4 files changed, 51 insertions(+), 69 deletions(-)

diff --git a/test/webapi/ows/coverages/test_controllers.py b/test/webapi/ows/coverages/test_controllers.py
index d93d36f95..00e0ec45e 100644
--- a/test/webapi/ows/coverages/test_controllers.py
+++ b/test/webapi/ows/coverages/test_controllers.py
@@ -38,7 +38,8 @@
     dtype_to_opengis_datatype,
     get_dataarray_description,
     get_units,
-    is_xy_order, transform_bbox,
+    is_xy_order,
+    transform_bbox,
 )
 
 
@@ -58,7 +59,7 @@ def test_get_coverage_data_tiff(self):
             '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'
@@ -116,7 +117,7 @@ def test_get_coverage_data_netcdf(self):
                     'time_bnds',
                     'conc_chl',
                     'kd489',
-                    'crs'
+                    'crs',
                 },
                 set(ds.variables),
             )
@@ -160,7 +161,7 @@ def test_get_coverage_data_time_slice_subset(self):
                     'time',
                     'time_bnds',
                     'conc_chl',
-                    'spatial_ref'
+                    'spatial_ref',
                 },
                 set(ds.variables),
             )
@@ -302,5 +303,5 @@ def test_is_xy(self):
     def test_transform_bbox_same_crs(self):
         self.assertEqual(
             bbox := [1, 2, 3, 4],
-            transform_bbox(bbox, crs := pyproj.CRS('EPSG:4326'), crs)
-        )
\ No newline at end of file
+            transform_bbox(bbox, crs := pyproj.CRS('EPSG:4326'), crs),
+        )
diff --git a/test/webapi/ows/coverages/test_scaling.py b/test/webapi/ows/coverages/test_scaling.py
index c01dff775..16ad7b889 100644
--- a/test/webapi/ows/coverages/test_scaling.py
+++ b/test/webapi/ows/coverages/test_scaling.py
@@ -10,7 +10,6 @@
 
 
 class ScalingTest(unittest.TestCase):
-
     def setUp(self):
         self.epsg4326 = pyproj.CRS('EPSG:4326')
         self.ds = xcube.core.new.new_cube()
@@ -22,25 +21,26 @@ def test_default_scaling(self):
     def test_no_data(self):
         with self.assertRaises(ApiError.NotFound):
             CoverageScaling(
-                CoverageRequest({}), self.epsg4326,
-                self.ds.isel(lat=slice(0, 0))
+                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())
+            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
+            CoverageRequest({'scale-factor': ['2']}), self.epsg4326, self.ds
         )
         self.assertEqual((2, 2), scaling.scale)
 
@@ -48,7 +48,7 @@ def test_scale_axes(self):
         scaling = CoverageScaling(
             CoverageRequest({'scale-axes': ['Lat(3),Lon(1.2)']}),
             self.epsg4326,
-            self.ds
+            self.ds,
         )
         self.assertEqual((1.2, 3), scaling.scale)
         self.assertEqual((300, 60), scaling.size)
@@ -57,7 +57,7 @@ def test_scale_size(self):
         scaling = CoverageScaling(
             CoverageRequest({'scale-size': ['Lat(90),Lon(240)']}),
             self.epsg4326,
-            self.ds
+            self.ds,
         )
         self.assertEqual((240, 90), scaling.size)
         self.assertEqual((1.5, 2), scaling.scale)
diff --git a/xcube/webapi/ows/coverages/controllers.py b/xcube/webapi/ows/coverages/controllers.py
index 996b84d87..1848c28fb 100644
--- a/xcube/webapi/ows/coverages/controllers.py
+++ b/xcube/webapi/ows/coverages/controllers.py
@@ -102,18 +102,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
@@ -136,20 +125,12 @@ def get_coverage_data(
     if request.bbox is not None:
         ds = _apply_bbox(ds, request.bbox, bbox_crs, always_xy=False)
 
-    # NB: a default scale factor of 1 is not compulsory. We may in future
-    # choose to downscale by default if a large coverage is requested.
-    scale_factor = 1 if request.scale_factor is None else request.scale_factor
-    _assert_coverage_size_ok(ds, scale_factor)
-
-    source_gm = GridMapping.from_dataset(ds, crs=native_crs)
-    target_gm = source_gm
-
+    scaling = CoverageScaling(request, final_crs, ds)
+    _assert_coverage_size_ok(scaling)
+    target_gm = source_gm = GridMapping.from_dataset(ds, crs=native_crs)
     if native_crs != final_crs:
         target_gm = target_gm.transform(final_crs).to_regular()
-
-    scaling = CoverageScaling(request, final_crs, ds)
     target_gm = scaling.apply(target_gm)
-
     if target_gm is not source_gm:
         ds = resample_in_space(ds, source_gm=source_gm, target_gm=target_gm)
 
@@ -194,25 +175,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}.'
         )
 
 
@@ -350,11 +334,6 @@ def get_bbox_from_ds(ds: xr.Dataset):
     return bbox
 
 
-def apply_scaling(gm: GridMapping, ds: xr.Dataset, request: CoverageRequest):
-    # TODO: implement me
-    pass
-
-
 def _find_geographic_parameters(
     names: list[str],
 ) -> tuple[Optional[str], Optional[str]]:
diff --git a/xcube/webapi/ows/coverages/scaling.py b/xcube/webapi/ows/coverages/scaling.py
index 45c568572..64d53bd93 100644
--- a/xcube/webapi/ows/coverages/scaling.py
+++ b/xcube/webapi/ows/coverages/scaling.py
@@ -30,7 +30,6 @@
 
 
 class CoverageScaling:
-
     _scale: Optional[tuple[float, float]] = None
     _final_size: Optional[tuple[int, int]] = None
     _initial_size: tuple[int, int]
@@ -38,8 +37,9 @@ class CoverageScaling:
     _x_name: str
     _y_name: str
 
-    def __init__(self, request: CoverageRequest, crs: pyproj.CRS,
-                 ds: xr.Dataset):
+    def __init__(
+        self, request: CoverageRequest, crs: pyproj.CRS, ds: xr.Dataset
+    ):
         h_dim = get_h_dim(ds)
         v_dim = get_v_dim(ds)
         for d in h_dim, v_dim:
@@ -100,7 +100,7 @@ def size(self) -> tuple[float, float]:
             return x_initial / x_scale, y_initial / y_scale
 
     def _get_xy_values(
-            self, axis_to_value: dict[str, float]
+        self, axis_to_value: dict[str, float]
     ) -> tuple[float, float]:
         x, y = None, None
         for axis in axis_to_value:
@@ -114,10 +114,12 @@ def get_axis_from_crs(self, valid_identifiers: set[str]):
         for axis in self._crs.axis_info:
             if not hasattr(axis, 'abbrev'):
                 continue
-            identifiers = set(map(
-                lambda attr: getattr(axis, attr, '').lower(),
-                ['name', 'abbrev']
-            ))
+            identifiers = set(
+                map(
+                    lambda attr: getattr(axis, attr, '').lower(),
+                    ['name', 'abbrev'],
+                )
+            )
             if identifiers & valid_identifiers:
                 return axis.abbrev
         return None

From de3ddba1bef12c7c61b6fab2f142cf2092f97a7c Mon Sep 17 00:00:00 2001
From: Pontus Lurcock <pontus.lurcock@brockmann-consult.de>
Date: Tue, 12 Dec 2023 12:26:13 +0100
Subject: [PATCH 08/14] Update changelog

---
 CHANGES.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CHANGES.md b/CHANGES.md
index 5863576eb..2b5482fa3 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

From cd3364fbe8d441f8e2ad4c8751499b43f2db46d0 Mon Sep 17 00:00:00 2001
From: Pontus Lurcock <pontus.lurcock@brockmann-consult.de>
Date: Wed, 13 Dec 2023 10:25:52 +0100
Subject: [PATCH 09/14] Coverages: regularize grid mapping before scaling

---
 xcube/webapi/ows/coverages/scaling.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/xcube/webapi/ows/coverages/scaling.py b/xcube/webapi/ows/coverages/scaling.py
index 64d53bd93..21440f757 100644
--- a/xcube/webapi/ows/coverages/scaling.py
+++ b/xcube/webapi/ows/coverages/scaling.py
@@ -128,4 +128,5 @@ def apply(self, gm: GridMapping):
         if self.scale == (1, 1):
             return gm
         else:
-            return gm.scale((1 / self.scale[0], 1 / self.scale[1]))
+            return \
+                gm.to_regular().scale((1 / self.scale[0], 1 / self.scale[1]))

From c14e9a9fcb6c82a96bb21f3f3df74b18069623a6 Mon Sep 17 00:00:00 2001
From: Pontus Lurcock <pontus.lurcock@brockmann-consult.de>
Date: Wed, 13 Dec 2023 15:16:44 +0100
Subject: [PATCH 10/14] OGC Coverages: improve scaling

Scaling is now applied after final application of bounding boxes, to
ensure that the final size is correct.
---
 xcube/webapi/ows/coverages/controllers.py | 36 +++++++++++++++--------
 xcube/webapi/ows/coverages/scaling.py     |  4 ++-
 2 files changed, 27 insertions(+), 13 deletions(-)

diff --git a/xcube/webapi/ows/coverages/controllers.py b/xcube/webapi/ows/coverages/controllers.py
index 1848c28fb..a7ea82ff6 100644
--- a/xcube/webapi/ows/coverages/controllers.py
+++ b/xcube/webapi/ows/coverages/controllers.py
@@ -127,19 +127,31 @@ def get_coverage_data(
 
     scaling = CoverageScaling(request, final_crs, ds)
     _assert_coverage_size_ok(scaling)
-    target_gm = source_gm = GridMapping.from_dataset(ds, crs=native_crs)
+    transformed_gm = source_gm = GridMapping.from_dataset(ds, crs=native_crs)
     if native_crs != final_crs:
-        target_gm = target_gm.transform(final_crs).to_regular()
-    target_gm = scaling.apply(target_gm)
-    if target_gm is not source_gm:
-        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)
+        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 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 scaling after bbox, subsetting, and CRS transformation, to make
+    # sure that the final size is correct.
+    if scaling.scale != (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():
diff --git a/xcube/webapi/ows/coverages/scaling.py b/xcube/webapi/ows/coverages/scaling.py
index 21440f757..d2387a316 100644
--- a/xcube/webapi/ows/coverages/scaling.py
+++ b/xcube/webapi/ows/coverages/scaling.py
@@ -128,5 +128,7 @@ def apply(self, gm: GridMapping):
         if self.scale == (1, 1):
             return gm
         else:
+            regular = gm.to_regular()
+            source = regular.size
             return \
-                gm.to_regular().scale((1 / self.scale[0], 1 / self.scale[1]))
+                regular.scale((self.size[0] / source[0], self.size[1] / source[1]))

From 2573d09a26e502be895c6b6ddea4b06d0a807f3c Mon Sep 17 00:00:00 2001
From: Pontus Lurcock <pontus.lurcock@brockmann-consult.de>
Date: Wed, 13 Dec 2023 16:10:15 +0100
Subject: [PATCH 11/14] OGC Coverages: more reliable scaling

- scaling is recalculated after reprojection, bounding box application,
  and subsetting, to make sure that scaling factors are calculated with
  reference to the correct initial area.

- The checks for an empty or excessively large coverage are now carried
  out twice -- once before initial determination of the grid mapping
  (to avoid giving an obscure error if the grid mapping can't be
  determined) and once before final application of the scaling (to
  catch problems caused by bbox application etc.)
---
 xcube/webapi/ows/coverages/controllers.py | 13 +++++++++----
 xcube/webapi/ows/coverages/scaling.py     |  5 +++--
 2 files changed, 12 insertions(+), 6 deletions(-)

diff --git a/xcube/webapi/ows/coverages/controllers.py b/xcube/webapi/ows/coverages/controllers.py
index a7ea82ff6..3c92e6d11 100644
--- a/xcube/webapi/ows/coverages/controllers.py
+++ b/xcube/webapi/ows/coverages/controllers.py
@@ -125,8 +125,11 @@ def get_coverage_data(
     if request.bbox is not None:
         ds = _apply_bbox(ds, request.bbox, bbox_crs, always_xy=False)
 
-    scaling = CoverageScaling(request, final_crs, ds)
-    _assert_coverage_size_ok(scaling)
+    # 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:
         transformed_gm = transformed_gm.transform(final_crs).to_regular()
@@ -146,8 +149,10 @@ def get_coverage_data(
         if subset_bbox is not None:
             ds = _apply_bbox(ds, subset_bbox, subset_crs, always_xy=True)
 
-    # Apply scaling after bbox, subsetting, and CRS transformation, to make
-    # sure that the final size is correct.
+    # 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.scale != (1, 1):
         cropped_gm = GridMapping.from_dataset(ds, crs=final_crs)
         scaled_gm = scaling.apply(cropped_gm)
diff --git a/xcube/webapi/ows/coverages/scaling.py b/xcube/webapi/ows/coverages/scaling.py
index d2387a316..f679c6c13 100644
--- a/xcube/webapi/ows/coverages/scaling.py
+++ b/xcube/webapi/ows/coverages/scaling.py
@@ -130,5 +130,6 @@ def apply(self, gm: GridMapping):
         else:
             regular = gm.to_regular()
             source = regular.size
-            return \
-                regular.scale((self.size[0] / source[0], self.size[1] / source[1]))
+            return regular.scale(
+                (self.size[0] / source[0], self.size[1] / source[1])
+            )

From 37b04130f626b9af4e8e89c2040663710f2466b2 Mon Sep 17 00:00:00 2001
From: Pontus Lurcock <pontus.lurcock@brockmann-consult.de>
Date: Thu, 14 Dec 2023 14:27:03 +0100
Subject: [PATCH 12/14] OGC Coverages: make sure GeoTIFFs have CRS info

---
 xcube/webapi/ows/coverages/controllers.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/xcube/webapi/ows/coverages/controllers.py b/xcube/webapi/ows/coverages/controllers.py
index 3c92e6d11..44f7e4928 100644
--- a/xcube/webapi/ows/coverages/controllers.py
+++ b/xcube/webapi/ows/coverages/controllers.py
@@ -25,6 +25,7 @@
 
 import numpy as np
 import pyproj
+import rasterio
 import xarray as xr
 
 from xcube.core.gridmapping import GridMapping
@@ -170,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:
@@ -423,7 +424,9 @@ 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
@@ -452,6 +455,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

From 11ff37f8b3abc3b9d0639ae2e3af5a2e33743dbb Mon Sep 17 00:00:00 2001
From: Pontus Lurcock <pontus.lurcock@brockmann-consult.de>
Date: Thu, 14 Dec 2023 17:56:34 +0100
Subject: [PATCH 13/14] OGC Coverages: improve test coverage

---
 test/webapi/ows/coverages/test_controllers.py | 18 ++++++++++++++++++
 test/webapi/ows/coverages/test_scaling.py     | 11 +++++++++++
 2 files changed, 29 insertions(+)

diff --git a/test/webapi/ows/coverages/test_controllers.py b/test/webapi/ows/coverages/test_controllers.py
index 00e0ec45e..cc8509f9b 100644
--- a/test/webapi/ows/coverages/test_controllers.py
+++ b/test/webapi/ows/coverages/test_controllers.py
@@ -71,6 +71,24 @@ def test_get_coverage_data_tiff(self):
             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'
+        )
+        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_netcdf(self):
         crs = 'OGC:CRS84'
         query = {
diff --git a/test/webapi/ows/coverages/test_scaling.py b/test/webapi/ows/coverages/test_scaling.py
index 16ad7b889..e8dbd4edf 100644
--- a/test/webapi/ows/coverages/test_scaling.py
+++ b/test/webapi/ows/coverages/test_scaling.py
@@ -61,3 +61,14 @@ def test_scale_size(self):
         )
         self.assertEqual((240, 90), scaling.size)
         self.assertEqual((1.5, 2), scaling.scale)
+
+    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),
+        )

From 37929816c7fb020dfec0c55128743d8da6aae5b5 Mon Sep 17 00:00:00 2001
From: Pontus Lurcock <pontus.lurcock@brockmann-consult.de>
Date: Wed, 3 Jan 2024 11:58:02 +0100
Subject: [PATCH 14/14] Improve OGC Coverages scaling implementation

Mainly in response to PR review

- Rename CoverageScaling.scale property to "factor"

- Add docstrings and comments to CoverageScaling class

- Add some missing type hints

- Add unit tests for the coverages.util module

- Update docstring for controllers.dataset_to_image

- Make ScalingTest.test_scale_factor more thorough
---
 test/webapi/ows/coverages/test_controllers.py |  1 +
 test/webapi/ows/coverages/test_scaling.py     |  9 +--
 test/webapi/ows/coverages/test_util.py        | 18 +++++
 xcube/webapi/ows/coverages/controllers.py     |  3 +-
 xcube/webapi/ows/coverages/scaling.py         | 70 +++++++++++++++++--
 xcube/webapi/ows/coverages/util.py            |  4 +-
 6 files changed, 94 insertions(+), 11 deletions(-)
 create mode 100644 test/webapi/ows/coverages/test_util.py

diff --git a/test/webapi/ows/coverages/test_controllers.py b/test/webapi/ows/coverages/test_controllers.py
index cc8509f9b..35898056b 100644
--- a/test/webapi/ows/coverages/test_controllers.py
+++ b/test/webapi/ows/coverages/test_controllers.py
@@ -91,6 +91,7 @@ def test_get_coverage_data_geo_subset(self):
 
     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'],
diff --git a/test/webapi/ows/coverages/test_scaling.py b/test/webapi/ows/coverages/test_scaling.py
index e8dbd4edf..115ea59d6 100644
--- a/test/webapi/ows/coverages/test_scaling.py
+++ b/test/webapi/ows/coverages/test_scaling.py
@@ -16,7 +16,7 @@ def setUp(self):
 
     def test_default_scaling(self):
         scaling = CoverageScaling(CoverageRequest({}), self.epsg4326, self.ds)
-        self.assertEqual((1, 1), scaling.scale)
+        self.assertEqual((1, 1), scaling.factor)
 
     def test_no_data(self):
         with self.assertRaises(ApiError.NotFound):
@@ -42,7 +42,8 @@ def test_scale_factor(self):
         scaling = CoverageScaling(
             CoverageRequest({'scale-factor': ['2']}), self.epsg4326, self.ds
         )
-        self.assertEqual((2, 2), scaling.scale)
+        self.assertEqual((2, 2), scaling.factor)
+        self.assertEqual((180, 90), scaling.size)
 
     def test_scale_axes(self):
         scaling = CoverageScaling(
@@ -50,7 +51,7 @@ def test_scale_axes(self):
             self.epsg4326,
             self.ds,
         )
-        self.assertEqual((1.2, 3), scaling.scale)
+        self.assertEqual((1.2, 3), scaling.factor)
         self.assertEqual((300, 60), scaling.size)
 
     def test_scale_size(self):
@@ -60,7 +61,7 @@ def test_scale_size(self):
             self.ds,
         )
         self.assertEqual((240, 90), scaling.size)
-        self.assertEqual((1.5, 2), scaling.scale)
+        self.assertEqual((1.5, 2), scaling.factor)
 
     def test_apply_identity_scaling(self):
         # noinspection PyTypeChecker
diff --git a/test/webapi/ows/coverages/test_util.py b/test/webapi/ows/coverages/test_util.py
new file mode 100644
index 000000000..260431171
--- /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))
diff --git a/xcube/webapi/ows/coverages/controllers.py b/xcube/webapi/ows/coverages/controllers.py
index 44f7e4928..940ba9bbd 100644
--- a/xcube/webapi/ows/coverages/controllers.py
+++ b/xcube/webapi/ows/coverages/controllers.py
@@ -154,7 +154,7 @@ def get_coverage_data(
     # 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.scale != (1, 1):
+    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)
@@ -433,6 +433,7 @@ def dataset_to_image(
 
     :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
     """
 
diff --git a/xcube/webapi/ows/coverages/scaling.py b/xcube/webapi/ows/coverages/scaling.py
index f679c6c13..d041fa3c9 100644
--- a/xcube/webapi/ows/coverages/scaling.py
+++ b/xcube/webapi/ows/coverages/scaling.py
@@ -30,6 +30,13 @@
 
 
 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]
@@ -40,6 +47,15 @@ class CoverageScaling:
     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:
@@ -82,7 +98,21 @@ def __init__(
             self._scale = (1, 1)
 
     @property
-    def scale(self) -> tuple[float, float]:
+    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:
@@ -92,6 +122,11 @@ def scale(self) -> tuple[float, float]:
 
     @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:
@@ -110,7 +145,19 @@ def _get_xy_values(
                 y = axis_to_value[axis]
         return x, y
 
-    def get_axis_from_crs(self, valid_identifiers: set[str]):
+    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
@@ -124,12 +171,27 @@ def get_axis_from_crs(self, valid_identifiers: set[str]):
                 return axis.abbrev
         return None
 
-    def apply(self, gm: GridMapping):
-        if self.scale == (1, 1):
+    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
index 0c23ee3a4..72e41316c 100644
--- a/xcube/webapi/ows/coverages/util.py
+++ b/xcube/webapi/ows/coverages/util.py
@@ -23,13 +23,13 @@
 import xarray as xr
 
 
-def get_h_dim(ds: xr.Dataset):
+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):
+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]