Merge pull request #10822 from rouault/fix_gdalwarp_one_one_blank_edge

gdalwarp: fix crash/infinite loop when using -tr one a 1x1 blank raster (3.8.0 regression)

apps/gdalwarp_lib.cpp

@@ -4428,10 +4428,14 @@ static GDALDatasetH GDALWarpCreateOutput(
             nPixels = static_cast<int>((psOptions->dfMaxX - psOptions->dfMinX +
                                         (psOptions->dfXRes / 2.0)) /
                                        psOptions->dfXRes);
+            if (nPixels == 0)
+                nPixels = 1;
             nLines = static_cast<int>(
                 (std::fabs(psOptions->dfMaxY - psOptions->dfMinY) +
                  (psOptions->dfYRes / 2.0)) /
                 psOptions->dfYRes);
+            if (nLines == 0)
+                nLines = 1;
             adfDstGeoTransform[0] = psOptions->dfMinX;
             adfDstGeoTransform[3] = psOptions->dfMaxY;
             adfDstGeoTransform[1] = psOptions->dfXRes;
@@ -4445,7 +4449,7 @@ static GDALDatasetH GDALWarpCreateOutput(
         {
             // Try to detect if the edge of the raster would be blank
             // Cf https://github.com/OSGeo/gdal/issues/7905
-            while (true)
+            while (nPixels > 1 || nLines > 1)
             {
                 UpdateGeoTransformandAndPixelLines();
 
@@ -4547,18 +4551,30 @@ static GDALDatasetH GDALWarpCreateOutput(
 
                 if (bTopBlankLine)
                 {
+                    if (psOptions->dfMaxY - psOptions->dfMinY <=
+                        2 * psOptions->dfYRes)
+                        break;
                     psOptions->dfMaxY -= psOptions->dfYRes;
                 }
                 if (bBottomBlankLine)
                 {
+                    if (psOptions->dfMaxY - psOptions->dfMinY <=
+                        2 * psOptions->dfYRes)
+                        break;
                     psOptions->dfMinY += psOptions->dfYRes;
                 }
                 if (bLeftBlankCol)
                 {
+                    if (psOptions->dfMaxX - psOptions->dfMinX <=
+                        2 * psOptions->dfXRes)
+                        break;
                     psOptions->dfMinX += psOptions->dfXRes;
                 }
                 if (bRightBlankCol)
                 {
+                    if (psOptions->dfMaxX - psOptions->dfMinX <=
+                        2 * psOptions->dfXRes)
+                        break;
                     psOptions->dfMaxX -= psOptions->dfXRes;
                 }
             }

autotest/utilities/test_gdalwarp_lib.py

@@ -4259,3 +4259,31 @@ def test_gdalwarp_lib_minus_180_plus_180_to_span_over_180(tmp_vsimem, extra_colu
     ) == src_ds.GetRasterBand(1).ReadRaster(
         0, 0, src_ds.RasterXSize // 2, src_ds.RasterYSize
     )
+
+
+###############################################################################
+# Test bugfix for https://lists.osgeo.org/pipermail/gdal-dev/2024-September/059512.html
+
+
+@pytest.mark.parametrize("with_tap", [True, False])
+def test_gdalwarp_lib_blank_edge_one_by_one(with_tap):
+
+    src_ds = gdal.GetDriverByName("MEM").Create("", 1, 1)
+    src_ds.SetGeoTransform([6.8688, 0.0009, 0, 51.3747, 0, -0.0009])
+    srs = osr.SpatialReference()
+    srs.SetFromUserInput("WGS84")
+    srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    src_ds.SetSpatialRef(srs)
+    options = "-f MEM -tr 1000 1000 -t_srs EPSG:32631"
+    if with_tap:
+        options += " -tap"
+    out_ds = gdal.Warp("", src_ds, options=options)
+    assert out_ds.RasterXSize == 1
+    assert out_ds.RasterYSize == 1
+    gt = out_ds.GetGeoTransform()
+    if with_tap:
+        assert gt == pytest.approx((769000.0, 1000.0, 0.0, 5699000.0, 0.0, -1000.0))
+    else:
+        assert gt == pytest.approx(
+            (769234.6506516202, 1000.0, 0.0, 5698603.782217737, 0.0, -1000.0)
+        )