diff --git a/alg/gdalcutline.cpp b/alg/gdalcutline.cpp index fd039fff8dc3..66389128cffe 100644 --- 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 index 85a72c9f61a1..19b2f1626f9a 100644 --- 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 index 4288593f83df..d2974bff3c92 100755 --- 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.