Warping with cutline and blending for large and complex cut polygons

[asgerpetersen]

When warping with a cutline and using the blend option the processing time is heavily dependent on the complexity of the cut polygon. If a large part of the cut polygon falls outside the raster a lot of time is wasted in the GEOSDistance calls.

I have for instance tried cutting together two rasters of 2km x 2km using a coastline polygon which is very complex and has a bounding box of approx 300km x 300km. In this case it takes ~120 minutes to do the warp. If I reduce the coastline polygon to the relevant area the running time is ~14 minutes.

This approach could be incorporated in GDAL cutline algorithm by using something like (in pseudo code):

env = envelope(raster)
expand(env, blending distance)
reduced_cut_poly = intersection(env, cut_poly)

...and the using the reduced_cut_poly instead of the original polygon.

This would dramatically decrease the running time of the case with large and complex polygon and it would hardly affect the running time of other cases (assuming that GEOS has reasonable shortcut tests built in).

[asgerpetersen]

Looked at the coastline again, and the geometry that is actually used for cutting has a bounding box more like 30km x 30km. So the cut polygon doesn't have to be extremely larger than the raster to make the above suggested change relevant.

[warmerdam]

I have implemented this in trunk (r15946) and it seems to work though I haven't evaluated the performance implications.

[asgerpetersen]

Hi Frank,

I have another suggestion for a performance improvement which is related to this ticket:

When testing for each pixel if it falls outside the cutline bbox it would be smart to use the bbox of the intersected poLines instead of hPolygon.

[warmerdam]

I have reviewed the current code as to the best of my understanding the envelope used for per-pixel tests in BlendMaskGenerator?() is computed from the clipped set of lines.

[asgerpetersen]

Ok, I am not a shark in C, but the envelope used seems calculated in line 107 by the call "OGR_G_GetEnvelope( hPolygon, &sEnvelope );" And as far as I can see hPolygon is the original cutline polygon. I think it should be the envelope of poLines.

[warmerdam]

You are completely correct ... I conflated line 107 and 106 and somehow thought the envelope calculation was on poLines. I will correct...