diff --git a/alg/gdalwarpkernel.cpp b/alg/gdalwarpkernel.cpp index 2c046b67f055..47c5bf0f58dd 100644 --- a/alg/gdalwarpkernel.cpp +++ b/alg/gdalwarpkernel.cpp @@ -6685,9 +6685,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 index 2d55e07c212a..8cbb25f2909e 100755 --- 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)