Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[Backport release/3.6] gdalwarp: speed-up warping with cutline when the source dataset or … #6914

Merged
merged 1 commit into from Dec 13, 2022
Merged

[Backport release/3.6] gdalwarp: speed-up warping with cutline when the source dataset or … #6914

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
32 changes: 32 additions & 0 deletions alg/gdalcutline.cpp
View file
Open in desktop
Expand Up @@ -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. */
Expand Down
37 changes: 37 additions & 0 deletions apps/gdalwarp_lib.cpp
View file
Open in desktop
Expand Up @@ -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. */
/* -------------------------------------------------------------------- */
Expand Down
40 changes: 40 additions & 0 deletions autotest/utilities/test_gdalwarp_lib.py
View file
Open in desktop
Expand Up @@ -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.

Expand Down