[gdal-dev] Error in gdalwarp

[Even Rouault]

On Saturday 22 October 2016 18:30:27 Paolo Cavallini wrote:
> Hi all,
> I'm running gdalwarp on data from
> 
> https://github.com/qgis/QGIS-Training-Data/tree/master/training_manual_data/
> processing/iterative
> 
> as
> 
> gdalwarp -ot Float32 -q -of GTiff -cutline watersheds.shp -co
> COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 -wo OPTIMIZE_SIZE=TRUE
> dem25.tif OUTPUT.tif
> 
> and I get
> 
> GDAL command output:
> Warning 1: Self-intersection at or near point 583.00000000000182 326
> ERROR 1: Cutline is not valid after transformation
> 
> Previously the same command was running fine. I also checked geom
> validity of the vector, and I coud not spot any error.
> 
> Is there an explanation?

Hi Paolo,

Since GDAL 2.1, validation of the cutline geometry once reprojected into 
source coordinates is done, and defaults into an error if the validation 
fails.

Here the issue is linked to the fact that when the cutline datasource is made 
of several polygons, a multipolygon is made with them by adding them as parts 
of the multipolygon. And if you do that on watersheds.shp, it results in a 
invalid multipolygon likely because the parts share common boundaries

You can revert to the previous behaviour and avoid the error to be fatal by 
adding
--config GDALWARP_IGNORE_BAD_CUTLINE YES
to the gdalwarp command line.

The clipping to the cutline might work or not depending on the nature of the 
invalidity (in your case that shouldn't be an issue). The check was added 
because there were use cases where unnoticed invalidities resulted in rather 
critical issues (for example when you need to clip out some regions because of 
regulations on forbidden zones, and an invalidity causes an area that should 
have been clipped out to not being clipped out).

Probably that in the use case you point, the way the multipolygon is built 
should be revised to use a proper unioning operation instead of the current 
rather naive way. Could be worth a ticket.

You can do that step manually with:

ogr2ogr cutline.shp -sql "select st_union(geometry) from watersheds" 
watersheds.shp -dialect sqlite


Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com