We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
BUG: Scale in rectilinear reproject produces different results (without dask)
Related:
- #1448
- #1450

### Expected behaviour

reproject behaves similar to GDAL defaults

### Actual behaviour

The behavior has slight differences with GDAL defaults

### Steps to reproduce the behaviour
- Original CRS = EPSG:32615
- Original Transform = Affine(5.0, 0.0, 451560.0, 0.0, -5.0, 4640690.0)

```python
import datacube
import rioxarray
import xarray
from rasterio.enum import Resampling

dc = datacube.Datacube()
old = xarray.open_dataset("old_demsm.nc", decode_coords="all")  # loaded with datacube==1.8.12 
new = dc.load( # datacube==1.8.13 
    product=...,
    measurements=["demsm"],
    geom=geom,
    output_crs="EPSG:6933",
    resolution=[-10, 10],
    resampling="cubic_spline",
).squeeze()
```

![image](https://github.com/opendatacube/datacube-core/assets/8699967/28737a10-f813-4b6d-a374-22801875ef36)

```python
orignal_rasters = dc.find_datasets(
    product=...,
    measurements=["demsm"],
    geom=geom,
)
original_demsm = rioxarray.open_rasterio(
    orignal_rasters[0].raster_path,
    masked=True,
).rio.reproject_match(
    match_data_array=new["demsm"],
    resampling=Resampling.cubic_spline,
).rio.clip(
    [geom],
    crs="EPSG:4326",
).squeeze()
```
![image](https://github.com/opendatacube/datacube-core/assets/8699967/c9430328-60d4-414e-a17e-2a2c17035bc0)

### Environment information

* Which ``datacube --version`` are you using? 1.8.13

cc @alfredoahds 

--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `datacube/storage/_read.py`
Content:
```
1 # This file is part of the Open Data Cube, see https://opendatacube.org for more information
2 #
3 # Copyright (c) 2015-2023 ODC Contributors
4 # SPDX-License-Identifier: Apache-2.0
5 """ Dataset -> Raster
6 """
7 from affine import Affine
8 import numpy as np
9 from typing import Optional, Tuple
10 
11 from ..utils.math import is_almost_int, valid_mask
12 
13 from ..utils.geometry import (
14     roi_shape,
15     roi_is_empty,
16     roi_is_full,
17     roi_pad,
18     GeoBox,
19     w_,
20     warp_affine,
21     rio_reproject,
22     compute_reproject_roi)
23 
24 from ..utils.geometry._warp import is_resampling_nn, Resampling, Nodata
25 from ..utils.geometry import gbox as gbx
26 
27 
28 def rdr_geobox(rdr) -> GeoBox:
29     """ Construct GeoBox from opened dataset reader.
30     """
31     h, w = rdr.shape
32     return GeoBox(w, h, rdr.transform, rdr.crs)
33 
34 
35 def can_paste(rr, stol=1e-3, ttol=1e-2):
36     """
37     Take result of compute_reproject_roi and check if can read(possibly with scale) and paste,
38     or do we need to read then reproject.
39 
40     :returns: (True, None) if one can just read and paste
41     :returns: (False, Reason) if pasting is not possible, so need to reproject after reading
42     """
43     if not rr.is_st:  # not linear or not Scale + Translation
44         return False, "not ST"
45 
46     scale = rr.scale
47     if not is_almost_int(scale, stol):  # non-integer scaling
48         return False, "non-integer scale"
49 
50     scale = np.round(scale)
51     A = rr.transform.linear           # src -> dst
52     A = A*Affine.scale(scale, scale)  # src.overview[scale] -> dst
53 
54     (sx, _, tx,  # tx, ty are in dst pixel space
55      _, sy, ty,
56      *_) = A
57 
58     if any(abs(abs(s) - 1) > stol
59            for s in (sx, sy)):  # not equal scaling across axis?
60         return False, "sx!=sy, probably"
61 
62     ny, nx = (n/scale
63               for n in roi_shape(rr.roi_src))
64 
65     # src_roi doesn't divide by scale properly:
66     #  example 3x7 scaled down by factor of 2
67     if not all(is_almost_int(n, stol) for n in (nx, ny)):
68         return False, "src_roi doesn't align for scale"
69 
70     # TODO: probably need to deal with sub-pixel translation here, if we want
71     # to ignore sub-pixel translation and dst roi is 1 pixel bigger than src it
72     # should still be ok to paste after cropping dst roi by one pixel on the
73     # appropriate side. As it stands sub-pixel translation will be ignored only
74     # in some cases.
75 
76     # scaled down shape doesn't match dst shape
77     s_shape = (int(ny), int(nx))
78     if s_shape != roi_shape(rr.roi_dst):
79         return False, "src_roi/scale != dst_roi"
80 
81     # final check: sub-pixel translation
82     if not all(is_almost_int(t, ttol) for t in (tx, ty)):
83         return False, "sub-pixel translation"
84 
85     return True, None
86 
87 
88 def pick_read_scale(scale: float, rdr=None, tol=1e-3):
89     assert scale > 0
90     # First find nearest integer scale
91     #    Scale down to nearest integer, unless we can scale up by less than tol
92     #
93     # 2.999999 -> 3
94     # 2.8 -> 2
95     # 0.3 -> 1
96 
97     if scale < 1:
98         return 1
99 
100     if is_almost_int(scale, tol):
101         scale = np.round(scale)
102 
103     scale = int(scale)
104 
105     if rdr is not None:
106         # TODO: check available overviews in rdr
107         pass
108 
109     return scale
110 
111 
112 def read_time_slice(rdr,
113                     dst: np.ndarray,
114                     dst_gbox: GeoBox,
115                     resampling: Resampling,
116                     dst_nodata: Nodata,
117                     extra_dim_index: Optional[int] = None) -> Tuple[slice, slice]:
118     """ From opened reader object read into `dst`
119 
120     :returns: affected destination region
121     """
122     assert dst.shape == dst_gbox.shape
123     src_gbox = rdr_geobox(rdr)
124 
125     rr = compute_reproject_roi(src_gbox, dst_gbox)
126 
127     if roi_is_empty(rr.roi_dst):
128         return rr.roi_dst
129 
130     is_nn = is_resampling_nn(resampling)
131     scale = pick_read_scale(rr.scale, rdr)
132     scale_x, scale_y = [pick_read_scale(s) for s in rr.scale2]
133 
134     paste_ok, _ = can_paste(rr, ttol=0.9 if is_nn else 0.01)
135 
136     def norm_read_args(roi, shape, extra_dim_index):
137         if roi_is_full(roi, rdr.shape):
138             roi = None
139 
140         if roi is None and shape == rdr.shape:
141             shape = None
142 
143         w = w_[roi]
144 
145         # Build 3D read window
146         # Note: Might be a good idea to natively support nD read windows.
147         if extra_dim_index is not None:
148             if w is None:
149                 w = ()
150             return (extra_dim_index,) + w, shape
151         else:
152             # 2D read window
153             return w, shape
154 
155     if paste_ok:
156         A = rr.transform.linear
157         sx, sy = A.a, A.e
158 
159         dst = dst[rr.roi_dst]
160         pix = rdr.read(*norm_read_args(rr.roi_src, dst.shape, extra_dim_index))
161 
162         if sx < 0:
163             pix = pix[:, ::-1]
164         if sy < 0:
165             pix = pix[::-1, :]
166 
167         if rdr.nodata is None:
168             np.copyto(dst, pix)
169         else:
170             np.copyto(dst, pix, where=valid_mask(pix, rdr.nodata))
171     else:
172         if rr.is_st:
173             # add padding on src/dst ROIs, it was set to tight bounds
174             # TODO: this should probably happen inside compute_reproject_roi
175             rr.roi_dst = roi_pad(rr.roi_dst, 1, dst_gbox.shape)
176             rr.roi_src = roi_pad(rr.roi_src, 1, src_gbox.shape)
177 
178         dst = dst[rr.roi_dst]
179         dst_gbox = dst_gbox[rr.roi_dst]
180         src_gbox = src_gbox[rr.roi_src]
181         if scale > 1:
182             src_gbox = gbx.zoom_out(src_gbox, scale)
183 
184         pix = rdr.read(*norm_read_args(rr.roi_src, src_gbox.shape, extra_dim_index))
185 
186         if rr.transform.linear is not None:
187             A = (~src_gbox.transform)*dst_gbox.transform
188             warp_affine(pix, dst, A, resampling,
189                         src_nodata=rdr.nodata, dst_nodata=dst_nodata)
190         else:
191             rio_reproject(pix, dst, src_gbox, dst_gbox, resampling,
192                           src_nodata=rdr.nodata, dst_nodata=dst_nodata,
193                           XSCALE=scale_x, YSCALE=scale_y)
194 
195     return rr.roi_dst
196 
197 
198 def read_time_slice_v2(rdr,
199                        dst_gbox: GeoBox,
200                        resampling: Resampling,
201                        dst_nodata: Nodata) -> Tuple[Optional[np.ndarray],
202                                                     Tuple[slice, slice]]:
203     """ From opened reader object read into `dst`
204 
205     :returns: pixels read and ROI of dst_gbox that was affected
206     """
207     # pylint: disable=too-many-locals
208     src_gbox = rdr_geobox(rdr)
209 
210     rr = compute_reproject_roi(src_gbox, dst_gbox)
211 
212     if roi_is_empty(rr.roi_dst):
213         return None, rr.roi_dst
214 
215     is_nn = is_resampling_nn(resampling)
216     scale = pick_read_scale(rr.scale, rdr)
217 
218     paste_ok, _ = can_paste(rr, ttol=0.9 if is_nn else 0.01)
219 
220     def norm_read_args(roi, shape):
221         if roi_is_full(roi, rdr.shape):
222             roi = None
223 
224         if roi is None and shape == rdr.shape:
225             shape = None
226 
227         return roi, shape
228 
229     if paste_ok:
230         read_shape = roi_shape(rr.roi_dst)
231         A = rr.transform.linear
232         sx, sy = A.a, A.e
233 
234         pix = rdr.read(*norm_read_args(rr.roi_src, read_shape)).result()
235 
236         if sx < 0:
237             pix = pix[:, ::-1]
238         if sy < 0:
239             pix = pix[::-1, :]
240 
241         # normalise nodata to be equal to `dst_nodata`
242         if rdr.nodata is not None and rdr.nodata != dst_nodata:
243             pix[pix == rdr.nodata] = dst_nodata
244 
245         dst = pix
246     else:
247         if rr.is_st:
248             # add padding on src/dst ROIs, it was set to tight bounds
249             # TODO: this should probably happen inside compute_reproject_roi
250             rr.roi_dst = roi_pad(rr.roi_dst, 1, dst_gbox.shape)
251             rr.roi_src = roi_pad(rr.roi_src, 1, src_gbox.shape)
252 
253         dst_gbox = dst_gbox[rr.roi_dst]
254         src_gbox = src_gbox[rr.roi_src]
255         if scale > 1:
256             src_gbox = gbx.zoom_out(src_gbox, scale)
257 
258         dst = np.full(dst_gbox.shape, dst_nodata, dtype=rdr.dtype)
259         pix = rdr.read(*norm_read_args(rr.roi_src, src_gbox.shape)).result()
260 
261         if rr.transform.linear is not None:
262             A = (~src_gbox.transform)*dst_gbox.transform
263             warp_affine(pix, dst, A, resampling,
264                         src_nodata=rdr.nodata, dst_nodata=dst_nodata)
265         else:
266             rio_reproject(pix, dst, src_gbox, dst_gbox, resampling,
267                           src_nodata=rdr.nodata, dst_nodata=dst_nodata)
268 
269     return dst, rr.roi_dst
270
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/datacube/storage/_read.py b/datacube/storage/_read.py
--- a/datacube/storage/_read.py
+++ b/datacube/storage/_read.py
@@ -129,7 +129,6 @@
 
     is_nn = is_resampling_nn(resampling)
     scale = pick_read_scale(rr.scale, rdr)
-    scale_x, scale_y = [pick_read_scale(s) for s in rr.scale2]
 
     paste_ok, _ = can_paste(rr, ttol=0.9 if is_nn else 0.01)
 
@@ -183,14 +182,32 @@
 
         pix = rdr.read(*norm_read_args(rr.roi_src, src_gbox.shape, extra_dim_index))
 
+        # XSCALE and YSCALE are (currently) undocumented arguments that rasterio passed through to
+        # GDAL.  Not using them results in very inaccurate warping in images with highly
+        # non-square (i.e. long and thin) aspect ratios.
+        #
+        # See https://github.com/OSGeo/gdal/issues/7750 as well as
+        # https://github.com/opendatacube/datacube-core/pull/1450 and
+        # https://github.com/opendatacube/datacube-core/issues/1456
+        #
+        # In theory we might be able to get better results for queries with significantly
+        # different vertical and horizontal scales, but explicitly using XSCALE=1, YSCALE=1
+        # appears to be most appropriate for most requests, and is demonstrably better
+        # than not setting them at all.
+        gdal_scale_params = {
+            "XSCALE": 1,
+            "YSCALE": 1,
+        }
         if rr.transform.linear is not None:
             A = (~src_gbox.transform)*dst_gbox.transform
             warp_affine(pix, dst, A, resampling,
-                        src_nodata=rdr.nodata, dst_nodata=dst_nodata)
+                        src_nodata=rdr.nodata, dst_nodata=dst_nodata,
+                        **gdal_scale_params)
         else:
+
             rio_reproject(pix, dst, src_gbox, dst_gbox, resampling,
                           src_nodata=rdr.nodata, dst_nodata=dst_nodata,
-                          XSCALE=scale_x, YSCALE=scale_y)
+                          **gdal_scale_params)
 
     return rr.roi_dst

{"golden_diff": "diff --git a/datacube/storage/_read.py b/datacube/storage/_read.py\n--- a/datacube/storage/_read.py\n+++ b/datacube/storage/_read.py\n@@ -129,7 +129,6 @@\n \n     is_nn = is_resampling_nn(resampling)\n     scale = pick_read_scale(rr.scale, rdr)\n-    scale_x, scale_y = [pick_read_scale(s) for s in rr.scale2]\n \n     paste_ok, _ = can_paste(rr, ttol=0.9 if is_nn else 0.01)\n \n@@ -183,14 +182,32 @@\n \n         pix = rdr.read(*norm_read_args(rr.roi_src, src_gbox.shape, extra_dim_index))\n \n+        # XSCALE and YSCALE are (currently) undocumented arguments that rasterio passed through to\n+        # GDAL.  Not using them results in very inaccurate warping in images with highly\n+        # non-square (i.e. long and thin) aspect ratios.\n+        #\n+        # See https://github.com/OSGeo/gdal/issues/7750 as well as\n+        # https://github.com/opendatacube/datacube-core/pull/1450 and\n+        # https://github.com/opendatacube/datacube-core/issues/1456\n+        #\n+        # In theory we might be able to get better results for queries with significantly\n+        # different vertical and horizontal scales, but explicitly using XSCALE=1, YSCALE=1\n+        # appears to be most appropriate for most requests, and is demonstrably better\n+        # than not setting them at all.\n+        gdal_scale_params = {\n+            \"XSCALE\": 1,\n+            \"YSCALE\": 1,\n+        }\n         if rr.transform.linear is not None:\n             A = (~src_gbox.transform)*dst_gbox.transform\n             warp_affine(pix, dst, A, resampling,\n-                        src_nodata=rdr.nodata, dst_nodata=dst_nodata)\n+                        src_nodata=rdr.nodata, dst_nodata=dst_nodata,\n+                        **gdal_scale_params)\n         else:\n+\n             rio_reproject(pix, dst, src_gbox, dst_gbox, resampling,\n                           src_nodata=rdr.nodata, dst_nodata=dst_nodata,\n-                          XSCALE=scale_x, YSCALE=scale_y)\n+                          **gdal_scale_params)\n \n     return rr.roi_dst\n", "issue": "BUG: Scale in rectilinear reproject produces different results (without dask)\nRelated:\r\n- #1448\r\n- #1450\r\n\r\n### Expected behaviour\r\n\r\nreproject behaves similar to GDAL defaults\r\n\r\n### Actual behaviour\r\n\r\nThe behavior has slight differences with GDAL defaults\r\n\r\n### Steps to reproduce the behaviour\r\n- Original CRS = EPSG:32615\r\n- Original Transform = Affine(5.0, 0.0, 451560.0, 0.0, -5.0, 4640690.0)\r\n\r\n```python\r\nimport datacube\r\nimport rioxarray\r\nimport xarray\r\nfrom rasterio.enum import Resampling\r\n\r\ndc = datacube.Datacube()\r\nold = xarray.open_dataset(\"old_demsm.nc\", decode_coords=\"all\")  # loaded with datacube==1.8.12 \r\nnew = dc.load( # datacube==1.8.13 \r\n    product=...,\r\n    measurements=[\"demsm\"],\r\n    geom=geom,\r\n    output_crs=\"EPSG:6933\",\r\n    resolution=[-10, 10],\r\n    resampling=\"cubic_spline\",\r\n).squeeze()\r\n```\r\n\r\n![image](https://github.com/opendatacube/datacube-core/assets/8699967/28737a10-f813-4b6d-a374-22801875ef36)\r\n\r\n```python\r\norignal_rasters = dc.find_datasets(\r\n    product=...,\r\n    measurements=[\"demsm\"],\r\n    geom=geom,\r\n)\r\noriginal_demsm = rioxarray.open_rasterio(\r\n    orignal_rasters[0].raster_path,\r\n    masked=True,\r\n).rio.reproject_match(\r\n    match_data_array=new[\"demsm\"],\r\n    resampling=Resampling.cubic_spline,\r\n).rio.clip(\r\n    [geom],\r\n    crs=\"EPSG:4326\",\r\n).squeeze()\r\n```\r\n![image](https://github.com/opendatacube/datacube-core/assets/8699967/c9430328-60d4-414e-a17e-2a2c17035bc0)\r\n\r\n### Environment information\r\n\r\n* Which ``datacube --version`` are you using? 1.8.13\r\n\r\ncc @alfredoahds \n", "before_files": [{"content": "# This file is part of the Open Data Cube, see https://opendatacube.org for more information\n#\n# Copyright (c) 2015-2023 ODC Contributors\n# SPDX-License-Identifier: Apache-2.0\n\"\"\" Dataset -> Raster\n\"\"\"\nfrom affine import Affine\nimport numpy as np\nfrom typing import Optional, Tuple\n\nfrom ..utils.math import is_almost_int, valid_mask\n\nfrom ..utils.geometry import (\n    roi_shape,\n    roi_is_empty,\n    roi_is_full,\n    roi_pad,\n    GeoBox,\n    w_,\n    warp_affine,\n    rio_reproject,\n    compute_reproject_roi)\n\nfrom ..utils.geometry._warp import is_resampling_nn, Resampling, Nodata\nfrom ..utils.geometry import gbox as gbx\n\n\ndef rdr_geobox(rdr) -> GeoBox:\n    \"\"\" Construct GeoBox from opened dataset reader.\n    \"\"\"\n    h, w = rdr.shape\n    return GeoBox(w, h, rdr.transform, rdr.crs)\n\n\ndef can_paste(rr, stol=1e-3, ttol=1e-2):\n    \"\"\"\n    Take result of compute_reproject_roi and check if can read(possibly with scale) and paste,\n    or do we need to read then reproject.\n\n    :returns: (True, None) if one can just read and paste\n    :returns: (False, Reason) if pasting is not possible, so need to reproject after reading\n    \"\"\"\n    if not rr.is_st:  # not linear or not Scale + Translation\n        return False, \"not ST\"\n\n    scale = rr.scale\n    if not is_almost_int(scale, stol):  # non-integer scaling\n        return False, \"non-integer scale\"\n\n    scale = np.round(scale)\n    A = rr.transform.linear           # src -> dst\n    A = A*Affine.scale(scale, scale)  # src.overview[scale] -> dst\n\n    (sx, _, tx,  # tx, ty are in dst pixel space\n     _, sy, ty,\n     *_) = A\n\n    if any(abs(abs(s) - 1) > stol\n           for s in (sx, sy)):  # not equal scaling across axis?\n        return False, \"sx!=sy, probably\"\n\n    ny, nx = (n/scale\n              for n in roi_shape(rr.roi_src))\n\n    # src_roi doesn't divide by scale properly:\n    #  example 3x7 scaled down by factor of 2\n    if not all(is_almost_int(n, stol) for n in (nx, ny)):\n        return False, \"src_roi doesn't align for scale\"\n\n    # TODO: probably need to deal with sub-pixel translation here, if we want\n    # to ignore sub-pixel translation and dst roi is 1 pixel bigger than src it\n    # should still be ok to paste after cropping dst roi by one pixel on the\n    # appropriate side. As it stands sub-pixel translation will be ignored only\n    # in some cases.\n\n    # scaled down shape doesn't match dst shape\n    s_shape = (int(ny), int(nx))\n    if s_shape != roi_shape(rr.roi_dst):\n        return False, \"src_roi/scale != dst_roi\"\n\n    # final check: sub-pixel translation\n    if not all(is_almost_int(t, ttol) for t in (tx, ty)):\n        return False, \"sub-pixel translation\"\n\n    return True, None\n\n\ndef pick_read_scale(scale: float, rdr=None, tol=1e-3):\n    assert scale > 0\n    # First find nearest integer scale\n    #    Scale down to nearest integer, unless we can scale up by less than tol\n    #\n    # 2.999999 -> 3\n    # 2.8 -> 2\n    # 0.3 -> 1\n\n    if scale < 1:\n        return 1\n\n    if is_almost_int(scale, tol):\n        scale = np.round(scale)\n\n    scale = int(scale)\n\n    if rdr is not None:\n        # TODO: check available overviews in rdr\n        pass\n\n    return scale\n\n\ndef read_time_slice(rdr,\n                    dst: np.ndarray,\n                    dst_gbox: GeoBox,\n                    resampling: Resampling,\n                    dst_nodata: Nodata,\n                    extra_dim_index: Optional[int] = None) -> Tuple[slice, slice]:\n    \"\"\" From opened reader object read into `dst`\n\n    :returns: affected destination region\n    \"\"\"\n    assert dst.shape == dst_gbox.shape\n    src_gbox = rdr_geobox(rdr)\n\n    rr = compute_reproject_roi(src_gbox, dst_gbox)\n\n    if roi_is_empty(rr.roi_dst):\n        return rr.roi_dst\n\n    is_nn = is_resampling_nn(resampling)\n    scale = pick_read_scale(rr.scale, rdr)\n    scale_x, scale_y = [pick_read_scale(s) for s in rr.scale2]\n\n    paste_ok, _ = can_paste(rr, ttol=0.9 if is_nn else 0.01)\n\n    def norm_read_args(roi, shape, extra_dim_index):\n        if roi_is_full(roi, rdr.shape):\n            roi = None\n\n        if roi is None and shape == rdr.shape:\n            shape = None\n\n        w = w_[roi]\n\n        # Build 3D read window\n        # Note: Might be a good idea to natively support nD read windows.\n        if extra_dim_index is not None:\n            if w is None:\n                w = ()\n            return (extra_dim_index,) + w, shape\n        else:\n            # 2D read window\n            return w, shape\n\n    if paste_ok:\n        A = rr.transform.linear\n        sx, sy = A.a, A.e\n\n        dst = dst[rr.roi_dst]\n        pix = rdr.read(*norm_read_args(rr.roi_src, dst.shape, extra_dim_index))\n\n        if sx < 0:\n            pix = pix[:, ::-1]\n        if sy < 0:\n            pix = pix[::-1, :]\n\n        if rdr.nodata is None:\n            np.copyto(dst, pix)\n        else:\n            np.copyto(dst, pix, where=valid_mask(pix, rdr.nodata))\n    else:\n        if rr.is_st:\n            # add padding on src/dst ROIs, it was set to tight bounds\n            # TODO: this should probably happen inside compute_reproject_roi\n            rr.roi_dst = roi_pad(rr.roi_dst, 1, dst_gbox.shape)\n            rr.roi_src = roi_pad(rr.roi_src, 1, src_gbox.shape)\n\n        dst = dst[rr.roi_dst]\n        dst_gbox = dst_gbox[rr.roi_dst]\n        src_gbox = src_gbox[rr.roi_src]\n        if scale > 1:\n            src_gbox = gbx.zoom_out(src_gbox, scale)\n\n        pix = rdr.read(*norm_read_args(rr.roi_src, src_gbox.shape, extra_dim_index))\n\n        if rr.transform.linear is not None:\n            A = (~src_gbox.transform)*dst_gbox.transform\n            warp_affine(pix, dst, A, resampling,\n                        src_nodata=rdr.nodata, dst_nodata=dst_nodata)\n        else:\n            rio_reproject(pix, dst, src_gbox, dst_gbox, resampling,\n                          src_nodata=rdr.nodata, dst_nodata=dst_nodata,\n                          XSCALE=scale_x, YSCALE=scale_y)\n\n    return rr.roi_dst\n\n\ndef read_time_slice_v2(rdr,\n                       dst_gbox: GeoBox,\n                       resampling: Resampling,\n                       dst_nodata: Nodata) -> Tuple[Optional[np.ndarray],\n                                                    Tuple[slice, slice]]:\n    \"\"\" From opened reader object read into `dst`\n\n    :returns: pixels read and ROI of dst_gbox that was affected\n    \"\"\"\n    # pylint: disable=too-many-locals\n    src_gbox = rdr_geobox(rdr)\n\n    rr = compute_reproject_roi(src_gbox, dst_gbox)\n\n    if roi_is_empty(rr.roi_dst):\n        return None, rr.roi_dst\n\n    is_nn = is_resampling_nn(resampling)\n    scale = pick_read_scale(rr.scale, rdr)\n\n    paste_ok, _ = can_paste(rr, ttol=0.9 if is_nn else 0.01)\n\n    def norm_read_args(roi, shape):\n        if roi_is_full(roi, rdr.shape):\n            roi = None\n\n        if roi is None and shape == rdr.shape:\n            shape = None\n\n        return roi, shape\n\n    if paste_ok:\n        read_shape = roi_shape(rr.roi_dst)\n        A = rr.transform.linear\n        sx, sy = A.a, A.e\n\n        pix = rdr.read(*norm_read_args(rr.roi_src, read_shape)).result()\n\n        if sx < 0:\n            pix = pix[:, ::-1]\n        if sy < 0:\n            pix = pix[::-1, :]\n\n        # normalise nodata to be equal to `dst_nodata`\n        if rdr.nodata is not None and rdr.nodata != dst_nodata:\n            pix[pix == rdr.nodata] = dst_nodata\n\n        dst = pix\n    else:\n        if rr.is_st:\n            # add padding on src/dst ROIs, it was set to tight bounds\n            # TODO: this should probably happen inside compute_reproject_roi\n            rr.roi_dst = roi_pad(rr.roi_dst, 1, dst_gbox.shape)\n            rr.roi_src = roi_pad(rr.roi_src, 1, src_gbox.shape)\n\n        dst_gbox = dst_gbox[rr.roi_dst]\n        src_gbox = src_gbox[rr.roi_src]\n        if scale > 1:\n            src_gbox = gbx.zoom_out(src_gbox, scale)\n\n        dst = np.full(dst_gbox.shape, dst_nodata, dtype=rdr.dtype)\n        pix = rdr.read(*norm_read_args(rr.roi_src, src_gbox.shape)).result()\n\n        if rr.transform.linear is not None:\n            A = (~src_gbox.transform)*dst_gbox.transform\n            warp_affine(pix, dst, A, resampling,\n                        src_nodata=rdr.nodata, dst_nodata=dst_nodata)\n        else:\n            rio_reproject(pix, dst, src_gbox, dst_gbox, resampling,\n                          src_nodata=rdr.nodata, dst_nodata=dst_nodata)\n\n    return dst, rr.roi_dst\n", "path": "datacube/storage/_read.py"}], "after_files": [{"content": "# This file is part of the Open Data Cube, see https://opendatacube.org for more information\n#\n# Copyright (c) 2015-2023 ODC Contributors\n# SPDX-License-Identifier: Apache-2.0\n\"\"\" Dataset -> Raster\n\"\"\"\nfrom affine import Affine\nimport numpy as np\nfrom typing import Optional, Tuple\n\nfrom ..utils.math import is_almost_int, valid_mask\n\nfrom ..utils.geometry import (\n    roi_shape,\n    roi_is_empty,\n    roi_is_full,\n    roi_pad,\n    GeoBox,\n    w_,\n    warp_affine,\n    rio_reproject,\n    compute_reproject_roi)\n\nfrom ..utils.geometry._warp import is_resampling_nn, Resampling, Nodata\nfrom ..utils.geometry import gbox as gbx\n\n\ndef rdr_geobox(rdr) -> GeoBox:\n    \"\"\" Construct GeoBox from opened dataset reader.\n    \"\"\"\n    h, w = rdr.shape\n    return GeoBox(w, h, rdr.transform, rdr.crs)\n\n\ndef can_paste(rr, stol=1e-3, ttol=1e-2):\n    \"\"\"\n    Take result of compute_reproject_roi and check if can read(possibly with scale) and paste,\n    or do we need to read then reproject.\n\n    :returns: (True, None) if one can just read and paste\n    :returns: (False, Reason) if pasting is not possible, so need to reproject after reading\n    \"\"\"\n    if not rr.is_st:  # not linear or not Scale + Translation\n        return False, \"not ST\"\n\n    scale = rr.scale\n    if not is_almost_int(scale, stol):  # non-integer scaling\n        return False, \"non-integer scale\"\n\n    scale = np.round(scale)\n    A = rr.transform.linear           # src -> dst\n    A = A*Affine.scale(scale, scale)  # src.overview[scale] -> dst\n\n    (sx, _, tx,  # tx, ty are in dst pixel space\n     _, sy, ty,\n     *_) = A\n\n    if any(abs(abs(s) - 1) > stol\n           for s in (sx, sy)):  # not equal scaling across axis?\n        return False, \"sx!=sy, probably\"\n\n    ny, nx = (n/scale\n              for n in roi_shape(rr.roi_src))\n\n    # src_roi doesn't divide by scale properly:\n    #  example 3x7 scaled down by factor of 2\n    if not all(is_almost_int(n, stol) for n in (nx, ny)):\n        return False, \"src_roi doesn't align for scale\"\n\n    # TODO: probably need to deal with sub-pixel translation here, if we want\n    # to ignore sub-pixel translation and dst roi is 1 pixel bigger than src it\n    # should still be ok to paste after cropping dst roi by one pixel on the\n    # appropriate side. As it stands sub-pixel translation will be ignored only\n    # in some cases.\n\n    # scaled down shape doesn't match dst shape\n    s_shape = (int(ny), int(nx))\n    if s_shape != roi_shape(rr.roi_dst):\n        return False, \"src_roi/scale != dst_roi\"\n\n    # final check: sub-pixel translation\n    if not all(is_almost_int(t, ttol) for t in (tx, ty)):\n        return False, \"sub-pixel translation\"\n\n    return True, None\n\n\ndef pick_read_scale(scale: float, rdr=None, tol=1e-3):\n    assert scale > 0\n    # First find nearest integer scale\n    #    Scale down to nearest integer, unless we can scale up by less than tol\n    #\n    # 2.999999 -> 3\n    # 2.8 -> 2\n    # 0.3 -> 1\n\n    if scale < 1:\n        return 1\n\n    if is_almost_int(scale, tol):\n        scale = np.round(scale)\n\n    scale = int(scale)\n\n    if rdr is not None:\n        # TODO: check available overviews in rdr\n        pass\n\n    return scale\n\n\ndef read_time_slice(rdr,\n                    dst: np.ndarray,\n                    dst_gbox: GeoBox,\n                    resampling: Resampling,\n                    dst_nodata: Nodata,\n                    extra_dim_index: Optional[int] = None) -> Tuple[slice, slice]:\n    \"\"\" From opened reader object read into `dst`\n\n    :returns: affected destination region\n    \"\"\"\n    assert dst.shape == dst_gbox.shape\n    src_gbox = rdr_geobox(rdr)\n\n    rr = compute_reproject_roi(src_gbox, dst_gbox)\n\n    if roi_is_empty(rr.roi_dst):\n        return rr.roi_dst\n\n    is_nn = is_resampling_nn(resampling)\n    scale = pick_read_scale(rr.scale, rdr)\n\n    paste_ok, _ = can_paste(rr, ttol=0.9 if is_nn else 0.01)\n\n    def norm_read_args(roi, shape, extra_dim_index):\n        if roi_is_full(roi, rdr.shape):\n            roi = None\n\n        if roi is None and shape == rdr.shape:\n            shape = None\n\n        w = w_[roi]\n\n        # Build 3D read window\n        # Note: Might be a good idea to natively support nD read windows.\n        if extra_dim_index is not None:\n            if w is None:\n                w = ()\n            return (extra_dim_index,) + w, shape\n        else:\n            # 2D read window\n            return w, shape\n\n    if paste_ok:\n        A = rr.transform.linear\n        sx, sy = A.a, A.e\n\n        dst = dst[rr.roi_dst]\n        pix = rdr.read(*norm_read_args(rr.roi_src, dst.shape, extra_dim_index))\n\n        if sx < 0:\n            pix = pix[:, ::-1]\n        if sy < 0:\n            pix = pix[::-1, :]\n\n        if rdr.nodata is None:\n            np.copyto(dst, pix)\n        else:\n            np.copyto(dst, pix, where=valid_mask(pix, rdr.nodata))\n    else:\n        if rr.is_st:\n            # add padding on src/dst ROIs, it was set to tight bounds\n            # TODO: this should probably happen inside compute_reproject_roi\n            rr.roi_dst = roi_pad(rr.roi_dst, 1, dst_gbox.shape)\n            rr.roi_src = roi_pad(rr.roi_src, 1, src_gbox.shape)\n\n        dst = dst[rr.roi_dst]\n        dst_gbox = dst_gbox[rr.roi_dst]\n        src_gbox = src_gbox[rr.roi_src]\n        if scale > 1:\n            src_gbox = gbx.zoom_out(src_gbox, scale)\n\n        pix = rdr.read(*norm_read_args(rr.roi_src, src_gbox.shape, extra_dim_index))\n\n        # XSCALE and YSCALE are (currently) undocumented arguments that rasterio passed through to\n        # GDAL.  Not using them results in very inaccurate warping in images with highly\n        # non-square (i.e. long and thin) aspect ratios.\n        #\n        # See https://github.com/OSGeo/gdal/issues/7750 as well as\n        # https://github.com/opendatacube/datacube-core/pull/1450 and\n        # https://github.com/opendatacube/datacube-core/issues/1456\n        #\n        # In theory we might be able to get better results for queries with significantly\n        # different vertical and horizontal scales, but explicitly using XSCALE=1, YSCALE=1\n        # appears to be most appropriate for most requests, and is demonstrably better\n        # than not setting them at all.\n        gdal_scale_params = {\n            \"XSCALE\": 1,\n            \"YSCALE\": 1,\n        }\n        if rr.transform.linear is not None:\n            A = (~src_gbox.transform)*dst_gbox.transform\n            warp_affine(pix, dst, A, resampling,\n                        src_nodata=rdr.nodata, dst_nodata=dst_nodata,\n                        **gdal_scale_params)\n        else:\n\n            rio_reproject(pix, dst, src_gbox, dst_gbox, resampling,\n                          src_nodata=rdr.nodata, dst_nodata=dst_nodata,\n                          **gdal_scale_params)\n\n    return rr.roi_dst\n\n\ndef read_time_slice_v2(rdr,\n                       dst_gbox: GeoBox,\n                       resampling: Resampling,\n                       dst_nodata: Nodata) -> Tuple[Optional[np.ndarray],\n                                                    Tuple[slice, slice]]:\n    \"\"\" From opened reader object read into `dst`\n\n    :returns: pixels read and ROI of dst_gbox that was affected\n    \"\"\"\n    # pylint: disable=too-many-locals\n    src_gbox = rdr_geobox(rdr)\n\n    rr = compute_reproject_roi(src_gbox, dst_gbox)\n\n    if roi_is_empty(rr.roi_dst):\n        return None, rr.roi_dst\n\n    is_nn = is_resampling_nn(resampling)\n    scale = pick_read_scale(rr.scale, rdr)\n\n    paste_ok, _ = can_paste(rr, ttol=0.9 if is_nn else 0.01)\n\n    def norm_read_args(roi, shape):\n        if roi_is_full(roi, rdr.shape):\n            roi = None\n\n        if roi is None and shape == rdr.shape:\n            shape = None\n\n        return roi, shape\n\n    if paste_ok:\n        read_shape = roi_shape(rr.roi_dst)\n        A = rr.transform.linear\n        sx, sy = A.a, A.e\n\n        pix = rdr.read(*norm_read_args(rr.roi_src, read_shape)).result()\n\n        if sx < 0:\n            pix = pix[:, ::-1]\n        if sy < 0:\n            pix = pix[::-1, :]\n\n        # normalise nodata to be equal to `dst_nodata`\n        if rdr.nodata is not None and rdr.nodata != dst_nodata:\n            pix[pix == rdr.nodata] = dst_nodata\n\n        dst = pix\n    else:\n        if rr.is_st:\n            # add padding on src/dst ROIs, it was set to tight bounds\n            # TODO: this should probably happen inside compute_reproject_roi\n            rr.roi_dst = roi_pad(rr.roi_dst, 1, dst_gbox.shape)\n            rr.roi_src = roi_pad(rr.roi_src, 1, src_gbox.shape)\n\n        dst_gbox = dst_gbox[rr.roi_dst]\n        src_gbox = src_gbox[rr.roi_src]\n        if scale > 1:\n            src_gbox = gbx.zoom_out(src_gbox, scale)\n\n        dst = np.full(dst_gbox.shape, dst_nodata, dtype=rdr.dtype)\n        pix = rdr.read(*norm_read_args(rr.roi_src, src_gbox.shape)).result()\n\n        if rr.transform.linear is not None:\n            A = (~src_gbox.transform)*dst_gbox.transform\n            warp_affine(pix, dst, A, resampling,\n                        src_nodata=rdr.nodata, dst_nodata=dst_nodata)\n        else:\n            rio_reproject(pix, dst, src_gbox, dst_gbox, resampling,\n                          src_nodata=rdr.nodata, dst_nodata=dst_nodata)\n\n    return dst, rr.roi_dst\n", "path": "datacube/storage/_read.py"}]}