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.

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

gdalwarp cutline performance #6905

Closed
blacha opened this issue Dec 12, 2022 · 4 comments · Fixed by #6907
Closed

gdalwarp cutline performance #6905

blacha opened this issue Dec 12, 2022 · 4 comments · Fixed by #6907
Assignees

Comments

@blacha
Copy link
Contributor

blacha commented Dec 12, 2022

Expected behavior and actual behavior.

I have a huge number of TIFFs I am processing into cloud optimised geotiffs, I am applying a cutline to remove bad information from the edge of the TIFFs (black/white spots/warped features).

If I apply a cutline to a tiff the processing time goes up significantly 2 seconds -> 72 seconds even if the cutline fully covers the TIFF.

Is there a easy way to determine if the cutline is needed before? Could GDAL automatically ignore cutlines if the the cutline fully covers the TIFF?

I am currently looking at adding this logic before applying the warp/COG

from osgeo import gdal
from osgeo import ogr


inds = gdal.Open("_SNC8650_BK32_5000_0810.tif", gdal.GA_ReadOnly)
incut = ogr.Open("./SNC8650-combined.fgb", gdal.GA_ReadOnly)

# Get raster geometry
transform = inds.GetGeoTransform()
pixelWidth = transform[1]
pixelHeight = transform[5]
cols = inds.RasterXSize
rows = inds.RasterYSize

xLeft = transform[0]
yTop = transform[3]
xRight = xLeft+cols*pixelWidth
yBottom = yTop+rows*pixelHeight

ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(xLeft, yTop)
ring.AddPoint(xLeft, yBottom)
ring.AddPoint(xRight, yBottom)
ring.AddPoint(xRight, yTop)
ring.AddPoint(xLeft, yTop)
rasterGeometry = ogr.Geometry(ogr.wkbPolygon)
rasterGeometry.AddGeometry(ring)

layer = incut.GetLayer()
feature = layer.GetFeature(0)
vectorGeometry = feature.GetGeometryRef()

inter = rasterGeometry.Intersection(vectorGeometry)
intersectionArea = inter.GetArea()
rasterArea = rasterGeometry.GetArea();

intersectionArea = intersectionArea / rasterArea 

if intersectionArea < 0.999:
  return apply_cutline(..)

Steps to reproduce the problem.

Data: https://public.lo.chard.com/gdal-6905/

I have a dataset of tiffs the black shaded boxes , and a cutline that is the red line
image

these commands operate on the red shaded tile

Create a VRT with a cutline and one without

gdalwarp -of VRT   images/_SNC8650_BK32_5000_0810.tif tmp/_SNC8650_BK32_5000_0810.nocut.vrt

gdalwarp -of VRT  -cutline images/SNC8650-combined.fgb images/_SNC8650_BK32_5000_0810.tif tmp/_SNC8650_BK32_5000_0810.cutline.vrt

Run the cog creation process on the cutline and non cutline variant

gdal_translate tmp/_SNC8650_BK32_5000_0810.cutline.vrt tmp/_SNC8650_BK32_5000_0810_cutline.tif -of COG -co num_threads=all_cpus -co COMPRESS=webp

gdal_translate tmp/_SNC8650_BK32_5000_0810.nocutline.vrt tmp/_SNC8650_BK32_5000_0810_nocutline.tif -of COG -co num_threads=all_cpus -co COMPRESS=webp

Operating system

Docker/Ubuntu

GDAL version and provenance

docker gdal
:ubuntu-small-latest
:ubuntu-small-3.6.0

rouault added a commit to rouault/gdal that referenced this issue Dec 12, 2022
…ith GeoTIFF keys override and northing, easting axis order (found when looking at OSGeo#6905)
@rouault rouault self-assigned this Dec 12, 2022
@blacha
Copy link
Contributor Author

blacha commented Dec 12, 2022

Another thought, We also found that if we optimize the cutline geometry down to just the tiff itself (or buffered slightly) it also significantly improved the performance of the warp/translate.

rouault added a commit to rouault/gdal that referenced this issue Dec 12, 2022
…ocessing chunks are fully contained within the cutline (fixes OSGeo#6905)
rouault added a commit to rouault/gdal that referenced this issue Dec 12, 2022
…ocessing chunks are fully contained within the cutline (fixes OSGeo#6905)
@rouault
Copy link
Member

rouault commented Dec 12, 2022

#6907 should fix your issue. With it:

$ time gdalwarp  _SNC8650_BK32_5000_0810.tif out.tif -overwrite 
2.228s
$ time gdalwarp -cutline SNC8650-combined.fgb _SNC8650_BK32_5000_0810.tif out.tif -overwrite 
2.366s

rouault added a commit to rouault/gdal that referenced this issue Dec 12, 2022
…ocessing chunks are fully contained within the cutline (fixes OSGeo#6905)
rouault added a commit to rouault/gdal that referenced this issue Dec 12, 2022
…ocessing chunks are fully contained within the cutline (fixes OSGeo#6905)
@blacha
Copy link
Contributor Author

blacha commented Dec 12, 2022

oh wow that was super fast thanks @rouault I am building a copy of that pr now!

@blacha
Copy link
Contributor Author

blacha commented Dec 12, 2022

Awesome 🎉 @rouault this does fix my problem!

$ time gdal_translate tmp/_SNC8650_BK32_5000_0810.cutline.vrt tmp/_SNC8650_BK32_5000_0810_cutline.tif -of COG -co num_threads=all_cpus -co COMPRESS=webp
Input file size is 6400, 9600
0...10...20...30...40...50...60...70...80...90...100 - done.

real    0m2.357s
user    0m19.790s
sys     0m1.690s

rouault added a commit that referenced this issue Dec 13, 2022
…ith GeoTIFF keys override and northing, easting axis order (found when looking at #6905)
rouault added a commit that referenced this issue Dec 13, 2022
…ocessing chunks are fully contained within the cutline (fixes #6905)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants