Always use XSCALE=1,YSCALE=1 in warp. (#1457)

* Use SCALEX=1,SCALEY=1 in both warp code-paths.
opendatacube · Jun 19, 2023 · 7c7f460 · 7c7f460
1 parent c8795e3
commit 7c7f460
Show file tree

Hide file tree

Showing 3 changed files with 24 additions and 4 deletions.
diff --git 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 norm_read_args(roi, shape, extra_dim_index):
 
         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
@@ -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`)
 
 

diff --git 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