Merge pull request OSGeo#9062 from rouault/fix_9059

gdalwarp: do not enable blank line detection when -tap and -te are specified (fixes OSGeo#9059)
rouault · Jan 15, 2024 · b6ab258 · b6ab258
2 parents b44c448 + 579845b
commit b6ab258
Show file tree

Hide file tree

Showing 2 changed files with 49 additions and 8 deletions.
diff --git a/apps/gdalwarp_lib.cpp b/apps/gdalwarp_lib.cpp
@@ -4237,16 +4237,27 @@ static GDALDatasetH GDALWarpCreateOutput(
             (psOptions->bCropToCutline &&
              psOptions->aosWarpOptions.FetchBool("CUTLINE_ALL_TOUCHED", false)))
         {
-            bDetectBlankBorders = true;
-
-            psOptions->dfMinX = floor(psOptions->dfMinX / psOptions->dfXRes) *
-                                psOptions->dfXRes;
+            if ((psOptions->bTargetAlignedPixels &&
+                 bNeedsSuggestedWarpOutput) ||
+                (psOptions->bCropToCutline &&
+                 psOptions->aosWarpOptions.FetchBool("CUTLINE_ALL_TOUCHED",
+                                                     false)))
+            {
+                bDetectBlankBorders = true;
+            }
+            constexpr double EPS = 1e-8;
+            psOptions->dfMinX =
+                floor(psOptions->dfMinX / psOptions->dfXRes + EPS) *
+                psOptions->dfXRes;
             psOptions->dfMaxX =
-                ceil(psOptions->dfMaxX / psOptions->dfXRes) * psOptions->dfXRes;
-            psOptions->dfMinY = floor(psOptions->dfMinY / psOptions->dfYRes) *
-                                psOptions->dfYRes;
+                ceil(psOptions->dfMaxX / psOptions->dfXRes - EPS) *
+                psOptions->dfXRes;
+            psOptions->dfMinY =
+                floor(psOptions->dfMinY / psOptions->dfYRes + EPS) *
+                psOptions->dfYRes;
             psOptions->dfMaxY =
-                ceil(psOptions->dfMaxY / psOptions->dfYRes) * psOptions->dfYRes;
+                ceil(psOptions->dfMaxY / psOptions->dfYRes - EPS) *
+                psOptions->dfYRes;
         }
 
         const auto UpdateGeoTransformandAndPixelLines = [&]()

diff --git a/autotest/utilities/test_gdalwarp_lib.py b/autotest/utilities/test_gdalwarp_lib.py
@@ -665,6 +665,36 @@ def test_gdalwarp_lib_32():
     ds = None
 
 
+###############################################################################
+# Test -tap, -tr and -te
+
+
+def test_gdalwarp_lib_tap_tr_te(tmp_vsimem):
+
+    src_filename = str(tmp_vsimem / "src.tif")
+    src_ds = gdal.GetDriverByName("GTiff").Create(
+        src_filename, 10980, 10980, 1, options=["SPARSE_OK=YES"]
+    )
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(32631)
+    src_ds.SetSpatialRef(srs)
+    src_ds.SetGeoTransform([600000, 10, 0, 5700000, 0, -10])
+
+    ds = gdal.Warp(
+        "",
+        src_ds,
+        format="VRT",
+        targetAlignedPixels=True,
+        xRes=10,
+        yRes=10,
+        outputBounds=[599800.0, 5590200.0, 709800.0, 5700000.0],
+    )
+    assert ds is not None
+    assert ds.GetGeoTransform() == pytest.approx(
+        (599800.0, 10.0, 0.0, 5700000.0, 0.0, -10.0)
+    )
+
+
 ###############################################################################
 # Test warping multiple sources