We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Distorted data along fractional Dask chunks when loading data with `dc.load` and Dask
### Expected behaviour
When loading data with Open Data Cube, we expect to be able to use Dask to efficiently parallelise workflows and reduce peak memory usage. Part of using Dask is an implicit assumption that data loaded by Dask will be identical to data loaded without Dask.

### Actual behaviour
There appears to be a bug in `datacube.load` that is causing data loaded with Dask to be distorted along edges of Dask chunks. We first noticed data being blurred/stretched along the right side of our data loads when loaded with Dask - see animation below which shows the effect along the right side of the image:
 
![dask_vs_nodask](https://github.com/opendatacube/datacube-core/assets/17680388/61d9b798-62af-42ee-b600-d4dc5a283be7)

Looking into this deeper, it appears that if a non-"nearest" resampling method is used (e.g. "cubic", "bilinear") and a Dask chunk extends off the right side of a data load, the data from that chunk will be strongly distorted. For example, if `dask_chunks={"time": 1, "x": 2048, "y": 2048}` and our data load is 2100 pixels wide, the right-most 52 pixels will be affected. If Dask is not used, `dc.load` works perfectly. 

### Steps to reproduce the behaviour
We've put together a Jupyter Notebook with a reproducible example here - it includes several examples of the issue, including some visualisations of how much this affects loaded data: https://gist.github.com/robbibt/4d6905d5d13137bb55362c4a6f1c18b8

For example, here is the difference between data loaded with and without Dask - note the band of strong differences along the right-hand side of the image:
![image](https://github.com/opendatacube/datacube-core/assets/17680388/a9609c86-ff5e-451c-a78c-771c65bb346b)


### Environment information

* Which ``datacube --version`` are you using? 1.8.12
* What datacube deployment/enviornment are you running against?
Prod DEA Sandbox


--- 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 
133     paste_ok, _ = can_paste(rr, ttol=0.9 if is_nn else 0.01)
134 
135     def norm_read_args(roi, shape, extra_dim_index):
136         if roi_is_full(roi, rdr.shape):
137             roi = None
138 
139         if roi is None and shape == rdr.shape:
140             shape = None
141 
142         w = w_[roi]
143 
144         # Build 3D read window
145         # Note: Might be a good idea to natively support nD read windows.
146         if extra_dim_index is not None:
147             if w is None:
148                 w = ()
149             return (extra_dim_index,) + w, shape
150         else:
151             # 2D read window
152             return w, shape
153 
154     if paste_ok:
155         A = rr.transform.linear
156         sx, sy = A.a, A.e
157 
158         dst = dst[rr.roi_dst]
159         pix = rdr.read(*norm_read_args(rr.roi_src, dst.shape, extra_dim_index))
160 
161         if sx < 0:
162             pix = pix[:, ::-1]
163         if sy < 0:
164             pix = pix[::-1, :]
165 
166         if rdr.nodata is None:
167             np.copyto(dst, pix)
168         else:
169             np.copyto(dst, pix, where=valid_mask(pix, rdr.nodata))
170     else:
171         if rr.is_st:
172             # add padding on src/dst ROIs, it was set to tight bounds
173             # TODO: this should probably happen inside compute_reproject_roi
174             rr.roi_dst = roi_pad(rr.roi_dst, 1, dst_gbox.shape)
175             rr.roi_src = roi_pad(rr.roi_src, 1, src_gbox.shape)
176 
177         dst = dst[rr.roi_dst]
178         dst_gbox = dst_gbox[rr.roi_dst]
179         src_gbox = src_gbox[rr.roi_src]
180         if scale > 1:
181             src_gbox = gbx.zoom_out(src_gbox, scale)
182 
183         pix = rdr.read(*norm_read_args(rr.roi_src, src_gbox.shape, extra_dim_index))
184 
185         if rr.transform.linear is not None:
186             A = (~src_gbox.transform)*dst_gbox.transform
187             warp_affine(pix, dst, A, resampling,
188                         src_nodata=rdr.nodata, dst_nodata=dst_nodata)
189         else:
190             rio_reproject(pix, dst, src_gbox, dst_gbox, resampling,
191                           src_nodata=rdr.nodata, dst_nodata=dst_nodata)
192 
193     return rr.roi_dst
194 
195 
196 def read_time_slice_v2(rdr,
197                        dst_gbox: GeoBox,
198                        resampling: Resampling,
199                        dst_nodata: Nodata) -> Tuple[Optional[np.ndarray],
200                                                     Tuple[slice, slice]]:
201     """ From opened reader object read into `dst`
202 
203     :returns: pixels read and ROI of dst_gbox that was affected
204     """
205     # pylint: disable=too-many-locals
206     src_gbox = rdr_geobox(rdr)
207 
208     rr = compute_reproject_roi(src_gbox, dst_gbox)
209 
210     if roi_is_empty(rr.roi_dst):
211         return None, rr.roi_dst
212 
213     is_nn = is_resampling_nn(resampling)
214     scale = pick_read_scale(rr.scale, rdr)
215 
216     paste_ok, _ = can_paste(rr, ttol=0.9 if is_nn else 0.01)
217 
218     def norm_read_args(roi, shape):
219         if roi_is_full(roi, rdr.shape):
220             roi = None
221 
222         if roi is None and shape == rdr.shape:
223             shape = None
224 
225         return roi, shape
226 
227     if paste_ok:
228         read_shape = roi_shape(rr.roi_dst)
229         A = rr.transform.linear
230         sx, sy = A.a, A.e
231 
232         pix = rdr.read(*norm_read_args(rr.roi_src, read_shape)).result()
233 
234         if sx < 0:
235             pix = pix[:, ::-1]
236         if sy < 0:
237             pix = pix[::-1, :]
238 
239         # normalise nodata to be equal to `dst_nodata`
240         if rdr.nodata is not None and rdr.nodata != dst_nodata:
241             pix[pix == rdr.nodata] = dst_nodata
242 
243         dst = pix
244     else:
245         if rr.is_st:
246             # add padding on src/dst ROIs, it was set to tight bounds
247             # TODO: this should probably happen inside compute_reproject_roi
248             rr.roi_dst = roi_pad(rr.roi_dst, 1, dst_gbox.shape)
249             rr.roi_src = roi_pad(rr.roi_src, 1, src_gbox.shape)
250 
251         dst_gbox = dst_gbox[rr.roi_dst]
252         src_gbox = src_gbox[rr.roi_src]
253         if scale > 1:
254             src_gbox = gbx.zoom_out(src_gbox, scale)
255 
256         dst = np.full(dst_gbox.shape, dst_nodata, dtype=rdr.dtype)
257         pix = rdr.read(*norm_read_args(rr.roi_src, src_gbox.shape)).result()
258 
259         if rr.transform.linear is not None:
260             A = (~src_gbox.transform)*dst_gbox.transform
261             warp_affine(pix, dst, A, resampling,
262                         src_nodata=rdr.nodata, dst_nodata=dst_nodata)
263         else:
264             rio_reproject(pix, dst, src_gbox, dst_gbox, resampling,
265                           src_nodata=rdr.nodata, dst_nodata=dst_nodata)
266 
267     return dst, rr.roi_dst
268
```
--- 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,6 +129,7 @@
 
     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)
 
@@ -188,7 +189,8 @@
                         src_nodata=rdr.nodata, dst_nodata=dst_nodata)
         else:
             rio_reproject(pix, dst, src_gbox, dst_gbox, resampling,
-                          src_nodata=rdr.nodata, dst_nodata=dst_nodata)
+                          src_nodata=rdr.nodata, dst_nodata=dst_nodata,
+                          XSCALE=scale_x, YSCALE=scale_y)
 
     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,6 +129,7 @@\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@@ -188,7 +189,8 @@\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+                          src_nodata=rdr.nodata, dst_nodata=dst_nodata,\n+                          XSCALE=scale_x, YSCALE=scale_y)\n \n     return rr.roi_dst\n", "issue": "Distorted data along fractional Dask chunks when loading data with `dc.load` and Dask\n### Expected behaviour\r\nWhen loading data with Open Data Cube, we expect to be able to use Dask to efficiently parallelise workflows and reduce peak memory usage. Part of using Dask is an implicit assumption that data loaded by Dask will be identical to data loaded without Dask.\r\n\r\n### Actual behaviour\r\nThere appears to be a bug in `datacube.load` that is causing data loaded with Dask to be distorted along edges of Dask chunks. We first noticed data being blurred/stretched along the right side of our data loads when loaded with Dask - see animation below which shows the effect along the right side of the image:\r\n \r\n![dask_vs_nodask](https://github.com/opendatacube/datacube-core/assets/17680388/61d9b798-62af-42ee-b600-d4dc5a283be7)\r\n\r\nLooking into this deeper, it appears that if a non-\"nearest\" resampling method is used (e.g. \"cubic\", \"bilinear\") and a Dask chunk extends off the right side of a data load, the data from that chunk will be strongly distorted. For example, if `dask_chunks={\"time\": 1, \"x\": 2048, \"y\": 2048}` and our data load is 2100 pixels wide, the right-most 52 pixels will be affected. If Dask is not used, `dc.load` works perfectly. \r\n\r\n### Steps to reproduce the behaviour\r\nWe've put together a Jupyter Notebook with a reproducible example here - it includes several examples of the issue, including some visualisations of how much this affects loaded data: https://gist.github.com/robbibt/4d6905d5d13137bb55362c4a6f1c18b8\r\n\r\nFor example, here is the difference between data loaded with and without Dask - note the band of strong differences along the right-hand side of the image:\r\n![image](https://github.com/opendatacube/datacube-core/assets/17680388/a9609c86-ff5e-451c-a78c-771c65bb346b)\r\n\r\n\r\n### Environment information\r\n\r\n* Which ``datacube --version`` are you using? 1.8.12\r\n* What datacube deployment/enviornment are you running against?\r\nProd DEA Sandbox\r\n\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\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\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    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"}]}