Warper: improve XSCALE/YSALE computation logic when the target shape is quite different from the source shape

The current logic to compute them is quite naive, and only reasonably
works when the shape of the output chunk is similar (in proportion and
orientation) to the one of the input chunk. Failing that, inappropriate
values are computed leading typically to blurring.

Rework the logic to be closer to a specific case that had already been
implemented to deal with discontinuities at the antimeridian, but with
modifications. So we compute the transformation of a maximum of 10x10
unit squares in the target pixel space to the source pixel space. If a
source polygon retains an approximative rectangular shape with modest
rotation, use the width/height of that rectangle to compute
XSCALE/YSCALE. Otherwise compute the area of the rectangle and make
XSCALE/YSCALE be the inverse of the square root of the area.
Average the values for each output polygon, except values that are less
than 10% of the maximum scale, to avoid issues with polar
discontinuities

Also revise the logic to decide whether we use a 4 input sample per output
pixel vs generalized resampling kernel for bilinear and bicubic.
Previously we only used the 4-sample formula if not downsampling (with a
margin error of 5%). Increase that tolerance to allow up to a 50%
downsampling. And emit a debug trace when the generalized formula is
used instead of the 4 sample one.

Due to those changes, a few existing expected checksums had to be
adjusted. And we also had to adjust a heuristics in the
average/mode/Q1/Q3/max/etc. resampling code path that dealt with
polar discontinuities to make a related test still pass.

Author: rouault

alg/gdalwarpkernel.cpp

@@ -2302,3 +2302,19 @@ def test_warp_zero_sized_target_extent():
     )
     assert out_ds.RasterXSize == 1
     assert out_ds.RasterYSize == 1
+
+
+@gdaltest.enable_exceptions()
+def test_warp_geolocation_array_with_rotation(tmp_path):
+    """Test that computed XSCALE/YSCALE are reasonable when a geolocation array
+    causes a ~ 90 degree rotation. Our input image is a checker of black and
+    white, and the expected output is also a non-blurred checker.
+    """
+    gdal.Warp(
+        tmp_path / "out.tif",
+        "data/input_for_geoloc_array_with_rotation.tif",
+        options="-r cubic -t_srs EPSG:32640 -to METHOD=GEOLOC_ARRAY -to GEOLOC_ARRAY=data/geoloc_array_with_rotation.tif -tr 30 30",
+    )
+    ds = gdal.Open(tmp_path / "out.tif")
+    ref_ds = gdal.Open("data/expected_output_for_geoloc_array_with_rotation.tif")
+    assert ds.GetRasterBand(1).Checksum() == ref_ds.GetRasterBand(1).Checksum()

autotest/alg/data/expected_output_for_geoloc_array_with_rotation.tif

@@ -485,11 +485,7 @@ def my_cbk(pct, _, arg):
         if gt[i] != pytest.approx(expected_gt[i], abs=1e-10 * abs(expected_gt[i])):
             assert False, gt
     got_cs = [ds.GetRasterBand(i + 1).Checksum() for i in range(3)]
-    assert got_cs in (
-        [26293, 23439, 14955],
-        [26228, 22085, 12992],
-        [25088, 23140, 13265],  # libjpeg 9e
-    )
+    assert got_cs == pytest.approx([20968, 23493, 14665], abs=100)
     assert ds.GetRasterBand(1).GetMaskBand().Checksum() == 17849
     assert ds.GetRasterBand(1).GetOverviewCount() == 0
     ds = None

autotest/alg/data/geoloc_array_with_rotation.tif

@@ -2270,16 +2270,14 @@ def test_gpkg_22(tile_drv_name):
 
 
 @pytest.mark.require_driver("PNG")
-def test_gpkg_26():
-
-    gdal.Unlink("/vsimem/tmp.gpkg")
-
-    tests = [
+@pytest.mark.parametrize(
+    "scheme,expected_cs,other_options",
+    [
         ("CUSTOM", [4672, 4672, 4672, 4873], None),
-        ("GoogleCRS84Quad", [3562, 3562, 3562, 3691], None),
-        ("GoogleCRS84Quad", [3562, 3562, 3562, 3691], ["RESAMPLING=BILINEAR"]),
-        ("GoogleCRS84Quad", [3417, 3417, 3417, 3691], ["RESAMPLING=CUBIC"]),
-        ("GoogleCRS84Quad", [3562, 3562, 3562, 3691], ["ZOOM_LEVEL_STRATEGY=AUTO"]),
+        ("GoogleCRS84Quad", [3439, 3439, 3439, 3691], None),
+        ("GoogleCRS84Quad", [3439, 3439, 3439, 3691], ["RESAMPLING=BILINEAR"]),
+        ("GoogleCRS84Quad", [3549, 3549, 3549, 3691], ["RESAMPLING=CUBIC"]),
+        ("GoogleCRS84Quad", [3439, 3439, 3439, 3691], ["ZOOM_LEVEL_STRATEGY=AUTO"]),
         (
             "GoogleCRS84Quad",
             [14445, 14445, 14445, 14448],
@@ -2295,77 +2293,65 @@ def test_gpkg_26():
             None,
             ["ZOOM_LEVEL=31"],
         ),
-        ("GoogleCRS84Quad", [3562, 3562, 3562, 3691], ["ZOOM_LEVEL_STRATEGY=LOWER"]),
+        ("GoogleCRS84Quad", [3439, 3439, 3439, 3691], ["ZOOM_LEVEL_STRATEGY=LOWER"]),
         ("GoogleMapsCompatible", [4118, 4118, 4118, 4406], None),
-        ("PseudoTMS_GlobalGeodetic", [3562, 3562, 3562, 3691], None),
+        ("PseudoTMS_GlobalGeodetic", [3439, 3439, 3439, 3691], None),
         ("PseudoTMS_GlobalMercator", [4118, 4118, 4118, 4406], None),
-    ]
-
-    for (scheme, expected_cs, other_options) in tests:
-
-        src_ds = gdal.Open("data/byte.tif")
-        options = ["TILE_FORMAT=PNG", "TILING_SCHEME=" + scheme]
-        if other_options:
-            options = options + other_options
-        if expected_cs is None:
-            with gdal.quiet_errors():
-                ds = gdaltest.gpkg_dr.CreateCopy(
-                    "/vsimem/tmp.gpkg", src_ds, options=options
-                )
-                assert ds is None
-                continue
-
-        ds = gdaltest.gpkg_dr.CreateCopy("/vsimem/tmp.gpkg", src_ds, options=options)
-        ds = None
+    ],
+)
+def test_gpkg_26(tmp_vsimem, scheme, expected_cs, other_options):
 
-        ds = gdal.OpenEx("/vsimem/tmp.gpkg", open_options=["BAND_COUNT=4"])
-        assert ds.GetMetadataItem("AREA_OR_POINT") == "Area"
-        got_cs = [ds.GetRasterBand(i + 1).Checksum() for i in range(4)]
-        # VC12 returns [3561, 3561, 3561, 3691] for GoogleCRS84Quad
-        # and For GoogleCRS84Quad RESAMPLING=CUBIC, got [3415, 3415, 3415, 3691]
-        if max([abs(got_cs[i] - expected_cs[i]) for i in range(4)]) > 2:
-            print(
-                "For %s, got %s, expected %s" % (scheme, str(got_cs), str(expected_cs))
+    src_ds = gdal.Open("data/byte.tif")
+    options = ["TILE_FORMAT=PNG", "TILING_SCHEME=" + scheme]
+    if other_options:
+        options = options + other_options
+    if expected_cs is None:
+        with gdal.quiet_errors():
+            ds = gdaltest.gpkg_dr.CreateCopy(
+                tmp_vsimem / "tmp.gpkg", src_ds, options=options
             )
-            assert gdal.GetConfigOption("APPVEYOR") is not None
-        ds = None
+            assert ds is None
+            return
 
-        gdal.Unlink("/vsimem/tmp.gpkg")
+    ds = gdaltest.gpkg_dr.CreateCopy(tmp_vsimem / "tmp.gpkg", src_ds, options=options)
+    ds = None
 
-    tests = [
+    ds = gdal.OpenEx(tmp_vsimem / "tmp.gpkg", open_options=["BAND_COUNT=4"])
+    assert ds.GetMetadataItem("AREA_OR_POINT") == "Area"
+    assert [ds.GetRasterBand(i + 1).Checksum() for i in range(4)] == pytest.approx(
+        expected_cs, abs=10
+    )
+    ds = None
+
+
+@pytest.mark.require_driver("PNG")
+@pytest.mark.parametrize(
+    "scheme,expected_cs,other_options",
+    [
         (
             "GoogleCRS84Quad",
-            [
-                [42255, 47336, 24963, 35707],
-                [42255, 47336, 24965, 35707],
-                [42253, 47333, 24961, 35707],
-                [42253, 47334, 24963, 35707],  # s390x
-            ],
+            [42255, 47336, 24963, 35707],
             None,
         ),
-        ("GoogleMapsCompatible", [[35429, 36787, 20035, 17849]], None),
-    ]
-
-    for (scheme, expected_cs, other_options) in tests:
+        ("GoogleMapsCompatible", [31249, 35596, 19001, 17849], None),
+    ],
+)
+def test_gpkg_26_bis(tmp_vsimem, scheme, expected_cs, other_options):
 
-        src_ds = gdal.Open("data/small_world.tif")
-        options = ["TILE_FORMAT=PNG", "TILING_SCHEME=" + scheme]
-        if other_options:
-            options = options + other_options
-        ds = gdaltest.gpkg_dr.CreateCopy("/vsimem/tmp.gpkg", src_ds, options=options)
-        ds = None
+    src_ds = gdal.Open("data/small_world.tif")
+    options = ["TILE_FORMAT=PNG", "TILING_SCHEME=" + scheme]
+    if other_options:
+        options = options + other_options
+    ds = gdaltest.gpkg_dr.CreateCopy(tmp_vsimem / "tmp.gpkg", src_ds, options=options)
+    ds = None
 
-        ds = gdal.Open("/vsimem/tmp.gpkg")
-        got_cs = [ds.GetRasterBand(i + 1).Checksum() for i in range(4)]
-        if got_cs not in expected_cs:
-            print(
-                "For %s, got %s, expected %s" % (scheme, str(got_cs), str(expected_cs))
-            )
-            assert gdal.GetConfigOption("APPVEYOR") is not None
-        ds = None
+    ds = gdal.Open(tmp_vsimem / "tmp.gpkg")
+    assert [ds.GetRasterBand(i + 1).Checksum() for i in range(4)] == pytest.approx(
+        expected_cs, abs=10
+    )
 
-        gdal.Unlink("/vsimem/tmp.gpkg")
 
+def test_gpkg_26_errors():
     # Test a few error cases
     with gdal.quiet_errors():
         ds = gdaltest.gpkg_dr.Create(

autotest/alg/data/input_for_geoloc_array_with_rotation.tif

@@ -1323,11 +1323,9 @@ def test_gdalalg_raster_tile_rgb(tmp_vsimem):
         assert ds.RasterCount == 3
         assert ds.RasterXSize == 256
         assert ds.RasterYSize == 256
-        assert [ds.GetRasterBand(i + 1).Checksum() for i in range(3)] == [
-            24650,
-            23280,
-            16559,
-        ]
+        assert [ds.GetRasterBand(i + 1).Checksum() for i in range(3)] == pytest.approx(
+            [22597, 22783, 16561], abs=30
+        )
 
 
 def test_gdalalg_raster_tile_rgba_all_opaque(tmp_vsimem):
@@ -1363,11 +1361,7 @@ def test_gdalalg_raster_tile_rgba_all_opaque(tmp_vsimem):
         assert ds.RasterXSize == 256
         assert ds.RasterYSize == 256
         assert [ds.GetRasterBand(i + 1).Checksum() for i in range(3)] == pytest.approx(
-            [
-                25111,
-                24737,
-                16108,
-            ],
+            [25111, 24737, 16107],
             abs=10,
         )
 
@@ -1504,11 +1498,9 @@ def test_gdalalg_raster_tile_rgba_no_alpha(tmp_vsimem):
         assert ds.RasterCount == 3
         assert ds.RasterXSize == 256
         assert ds.RasterYSize == 256
-        assert [ds.GetRasterBand(i + 1).Checksum() for i in range(3)] == [
-            24650,
-            23280,
-            16559,
-        ]
+        assert [ds.GetRasterBand(i + 1).Checksum() for i in range(3)] == pytest.approx(
+            [22597, 22783, 16561], abs=30
+        )
 
 
 def test_gdalalg_raster_tile_max_zoom(tmp_vsimem):

autotest/alg/warp.py

@@ -2497,7 +2497,7 @@ def test_gdalwarp_lib_ct():
         options='-r cubic -f MEM -t_srs EPSG:4326 -ct "proj=pipeline step inv proj=utm zone=11 ellps=clrk66 step proj=unitconvert xy_in=rad xy_out=deg step proj=axisswap order=2,1"',
     )
 
-    assert dstDS.GetRasterBand(1).Checksum() == 4705, "Bad checksum"
+    assert dstDS.GetRasterBand(1).Checksum() == 4772
 
 
 def test_gdalwarp_lib_ct_wkt():
@@ -2646,7 +2646,7 @@ def test_gdalwarp_lib_ct_wkt():
         coordinateOperation=wkt,
     )
 
-    assert dstDS.GetRasterBand(1).Checksum() == 4705, "Bad checksum"
+    assert dstDS.GetRasterBand(1).Checksum() == 4772
 
 
 ###############################################################################