Repository: opendatacube/datacube-core

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 

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/opendatacube/datacube-core/pull/1457?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube) Report
Patch coverage has no change and project coverage change: **`-0.01`** :warning:
> Comparison is base [(`c8795e3`)](https://app.codecov.io/gh/opendatacube/datacube-core/commit/c8795e3144fcf6cdabed23eecc9e8cb04626bcef?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube) 91.82% compared to head [(`19461c3`)](https://app.codecov.io/gh/opendatacube/datacube-core/pull/1457?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube) 91.82%.

<details><summary>Additional details and impacted files</summary>


```diff
@@             Coverage Diff             @@
##           develop    #1457      +/-   ##
===========================================
- Coverage    91.82%   91.82%   -0.01%     
===========================================
  Files          132      132              
  Lines        14516    14515       -1     
===========================================
- Hits         13330    13329       -1     
  Misses        1186     1186              
```


| [Impacted Files](https://app.codecov.io/gh/opendatacube/datacube-core/pull/1457?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube) | Coverage Δ | |
|---|---|---|
| [datacube/storage/\_read.py](https://app.codecov.io/gh/opendatacube/datacube-core/pull/1457?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube#diff-ZGF0YWN1YmUvc3RvcmFnZS9fcmVhZC5weQ==) | `100.00% <ø> (ø)` | |


</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/opendatacube/datacube-core/pull/1457?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube).   
:loudspeaker: Do you have feedback about the report comment? [Let us know in this issue](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube).


Patch:
diff --git a/datacube/storage/_read.py b/datacube/storage/_read.py
index 2c1dd05f..fb5b03b8 100644
--- a/datacube/storage/_read.py
+++ b/datacube/storage/_read.py
@@ -129,7 +129,6 @@ def read_time_slice(rdr,
 
     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 @@ def read_time_slice(rdr,
 
         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
 
diff --git a/docs/about/whats_new.rst b/docs/about/whats_new.rst
index 0b5d561a..3747a54c 100644
--- a/docs/about/whats_new.rst
+++ b/docs/about/whats_new.rst
@@ -8,6 +8,8 @@ What's New
 v1.8.next
 =========
 
+- Second attempt to address unexpected handling of image aspect ratios in rasterio and
+  GDAL. (:pull:`1457`)
 - Fix broken pypi publishing Github action (:pull:`1454`)
 
 


Test patch:
diff --git a/tests/storage/test_storage.py b/tests/storage/test_storage.py
index 37c981b7..b784f876 100644
--- a/tests/storage/test_storage.py
+++ b/tests/storage/test_storage.py
@@ -233,7 +233,8 @@ def assert_same_read_results(source, dst_shape, dst_dtype, dst_transform, dst_no
                                 dst_transform=dst_transform,
                                 dst_crs=str(dst_projection),
                                 dst_nodata=dst_nodata,
-                                resampling=resampling)
+                                resampling=resampling,
+                                XSCALE=1, YSCALE=1)
 
     result = np.full(dst_shape, dst_nodata, dtype=dst_dtype)
     H, W = dst_shape