Merge pull request #6914 from OSGeo/backport-6907-to-release/3.6

[Backport release/3.6] gdalwarp: speed-up warping with cutline when the source dataset or …
OSGeo · Dec 13, 2022 · 25ecdaa · 25ecdaa
2 parents f25aee2 + f2183a0
commit 25ecdaa
Show file tree

Hide file tree

Showing 3 changed files with 109 additions and 0 deletions.
diff --git a/alg/gdalcutline.cpp b/alg/gdalcutline.cpp
@@ -333,6 +333,38 @@ GDALWarpCutlineMasker( void *pMaskFuncArg,
         return CE_None;
     }
 
+    // And now check if the chunk to warp is fully contained within the cutline
+    // to save rasterization.
+    if( OGRGeometryFactory::haveGEOS()
+#ifdef DEBUG
+        // Env var just for debugging purposes
+         && !CPLTestBool(CPLGetConfigOption("GDALCUTLINE_SKIP_CONTAINMENT_TEST", "NO"))
+#endif
+    )
+    {
+        OGRLinearRing* poRing = new OGRLinearRing();
+        poRing->addPoint(-psWO->dfCutlineBlendDist + nXOff,
+                         -psWO->dfCutlineBlendDist + nYOff);
+        poRing->addPoint(-psWO->dfCutlineBlendDist + nXOff,
+                         psWO->dfCutlineBlendDist + nYOff + nYSize);
+        poRing->addPoint(psWO->dfCutlineBlendDist + nXOff + nXSize,
+                         psWO->dfCutlineBlendDist + nYOff + nYSize);
+        poRing->addPoint(psWO->dfCutlineBlendDist + nXOff + nXSize,
+                         -psWO->dfCutlineBlendDist + nYOff);
+        poRing->addPoint(-psWO->dfCutlineBlendDist + nXOff,
+                         -psWO->dfCutlineBlendDist + nYOff);
+        OGRPolygon oChunkFootprint;
+        oChunkFootprint.addRingDirectly(poRing);
+        OGREnvelope sChunkEnvelope;
+        oChunkFootprint.getEnvelope(&sChunkEnvelope );
+        if( sEnvelope.Contains(sChunkEnvelope) &&
+            OGRGeometry::FromHandle(hPolygon)->Contains(&oChunkFootprint) )
+        {
+            CPLDebug("WARP", "Source chunk fully contained within cutline.");
+            return CE_None;
+        }
+    }
+
 /* -------------------------------------------------------------------- */
 /*      Create a byte buffer into which we can burn the                 */
 /*      mask polygon and wrap it up as a memory dataset.                */

diff --git a/apps/gdalwarp_lib.cpp b/apps/gdalwarp_lib.cpp
@@ -4422,6 +4422,43 @@ TransformCutlineToSource( GDALDatasetH hSrcDS, OGRGeometryH hCutline,
         return CE_Failure;
     }
 
+    // Optimization: if the cutline contains the footprint of the source
+    // dataset, no need to use the cutline.
+    if( OGRGeometryFactory::haveGEOS()
+#ifdef DEBUG
+        // Env var just for debugging purposes
+        && !CPLTestBool(CPLGetConfigOption("GDALWARP_SKIP_CUTLINE_CONTAINMENT_TEST", "NO"))
+#endif
+        )
+    {
+        const double dfCutlineBlendDist = CPLAtof(
+            CSLFetchNameValueDef( *ppapszWarpOptions, "CUTLINE_BLEND_DIST", "0") );
+        OGRLinearRing* poRing = new OGRLinearRing();
+        poRing->addPoint(-dfCutlineBlendDist,
+                         -dfCutlineBlendDist);
+        poRing->addPoint(-dfCutlineBlendDist,
+                         dfCutlineBlendDist + GDALGetRasterYSize(hSrcDS));
+        poRing->addPoint(dfCutlineBlendDist + GDALGetRasterXSize(hSrcDS),
+                         dfCutlineBlendDist + GDALGetRasterYSize(hSrcDS));
+        poRing->addPoint(dfCutlineBlendDist + GDALGetRasterXSize(hSrcDS),
+                         -dfCutlineBlendDist);
+        poRing->addPoint(-dfCutlineBlendDist,
+                         -dfCutlineBlendDist);
+        OGRPolygon oSrcDSFootprint;
+        oSrcDSFootprint.addRingDirectly(poRing);
+        OGREnvelope sSrcDSEnvelope;
+        oSrcDSFootprint.getEnvelope(&sSrcDSEnvelope );
+        OGREnvelope sCutlineEnvelope;
+        OGRGeometry::FromHandle(hMultiPolygon)->getEnvelope(&sCutlineEnvelope );
+        if( sCutlineEnvelope.Contains(sSrcDSEnvelope) &&
+            OGRGeometry::FromHandle(hMultiPolygon)->Contains(&oSrcDSFootprint) )
+        {
+            CPLDebug("WARP", "Source dataset fully contained within cutline.");
+            OGR_G_DestroyGeometry( hMultiPolygon );
+            return CE_None;
+        }
+    }
+
 /* -------------------------------------------------------------------- */
 /*      Convert aggregate geometry into WKT.                            */
 /* -------------------------------------------------------------------- */

diff --git a/autotest/utilities/test_gdalwarp_lib.py b/autotest/utilities/test_gdalwarp_lib.py
@@ -462,6 +462,46 @@ def test_gdalwarp_lib_21():
     ds = None
 
 
+###############################################################################
+# Test cutline whose extent is larger than the source data
+
+
+@pytest.mark.parametrize(
+    "options", [{}, {"GDALWARP_SKIP_CUTLINE_CONTAINMENT_TEST": "YES"}]
+)
+def test_gdalwarp_lib_cutline_larger_source_dataset(options):
+
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(26711)
+    cutline_ds = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(
+        "/vsimem/cutline.shp"
+    )
+    cutline_lyr = cutline_ds.CreateLayer("cutline", srs=srs)
+    f = ogr.Feature(cutline_lyr.GetLayerDefn())
+    f.SetGeometry(
+        ogr.CreateGeometryFromWkt(
+            "POLYGON((400000 3000000,400000 4000000,500000 4000000,500000 3000000,400000 3000000))"
+        )
+    )
+    cutline_lyr.CreateFeature(f)
+    cutline_ds = None
+
+    with gdaltest.config_options(options):
+        ds = gdal.Warp(
+            "",
+            "../gcore/data/byte.tif",
+            format="MEM",
+            cutlineDSName="/vsimem/cutline.shp",
+            cutlineLayer="cutline",
+        )
+    assert ds is not None
+
+    assert ds.GetRasterBand(1).Checksum() == 4672
+
+    ds = None
+    ogr.GetDriverByName("ESRI Shapefile").DeleteDataSource("/vsimem/cutline.shp")
+
+
 ###############################################################################
 # Test cutline with ALL_TOUCHED enabled.