Merge pull request #10899 from OSGeo/backport-10894-to-release/3.9

[Backport release/3.9] Warper: fix too lax heuristics about antimeridian warping for Avg/Sum/Q1/Q3/Mode algorithms
OSGeo · Sep 30, 2024 · a5a8514 · a5a8514
2 parents 88cedd8 + 4f96dee
commit a5a8514
Show file tree

Hide file tree

Showing 2 changed files with 39 additions and 3 deletions.
diff --git a/alg/gdalwarpkernel.cpp b/alg/gdalwarpkernel.cpp
@@ -6697,9 +6697,26 @@ static void GWKAverageOrModeThread(void *pData)
                 (nSrcXSize - padfX2[iDstX]) * poWK->dfXScale <
                     nThresholdWrapOverX)
             {
-                bWrapOverX = true;
-                std::swap(padfX[iDstX], padfX2[iDstX]);
-                padfX2[iDstX] += nSrcXSize;
+                // Check there is a discontinuity by checking at mid-pixel.
+                // NOTE: all this remains fragile. To confidently
+                // detect antimeridian warping we should probably try to access
+                // georeferenced coordinates, and not rely only on tests on
+                // image space coordinates. But accessing georeferenced
+                // coordinates from here is not trivial, and we would for example
+                // have to handle both geographic, Mercator, etc.
+                // Let's hope this heuristics is good enough for now.
+                double x = iDstX + 0.5 + poWK->nDstXOff;
+                double y = iDstY + poWK->nDstYOff;
+                double z = 0;
+                int bSuccess = FALSE;
+                poWK->pfnTransformer(psJob->pTransformerArg, TRUE, 1, &x, &y,
+                                     &z, &bSuccess);
+                if (bSuccess && x < padfX[iDstX])
+                {
+                    bWrapOverX = true;
+                    std::swap(padfX[iDstX], padfX2[iDstX]);
+                    padfX2[iDstX] += nSrcXSize;
+                }
             }
 
             const double dfXMin = padfX[iDstX] - poWK->nSrcXOff;

diff --git a/autotest/utilities/test_gdalwarp_lib.py b/autotest/utilities/test_gdalwarp_lib.py
@@ -4287,3 +4287,22 @@ def test_gdalwarp_lib_blank_edge_one_by_one(with_tap):
         assert gt == pytest.approx(
             (769234.6506516202, 1000.0, 0.0, 5698603.782217737, 0.0, -1000.0)
         )
+
+
+###############################################################################
+# Test bugfix for https://github.com/OSGeo/gdal/issues/10892
+
+
+def test_gdalwarp_lib_average_ten_ten_to_one_one():
+
+    src_ds = gdal.GetDriverByName("MEM").Create("", 10, 10)
+    src_ds.SetGeoTransform([0, 1, 0, 0, 0, -1])
+    srs = osr.SpatialReference()
+    srs.SetFromUserInput("WGS84")
+    srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    src_ds.SetSpatialRef(srs)
+    src_ds.GetRasterBand(1).Fill(1)
+    out_ds = gdal.Warp(
+        "", src_ds, width=1, height=1, resampleAlg=gdal.GRIORA_Average, format="MEM"
+    )
+    assert out_ds.GetRasterBand(1).ComputeRasterMinMax() == (1, 1)