gdalwarp cutline with GCP raster input problem

[klokan]

Hi,

when I use gdalwarp cutline on a raster which is georeferenced by GCPs the result is not correct:

What I did:

gdalwarp -cutline 13218_1.shp -cl 13218_1 13218_1.KAP 13218_1.tif

the cutline location differs from:

gdalwarp 13218_1.KAP 13218_1.warp.tif
gdal_rasterize -burn 0 -l 13218_1 13218_1.shp 13218_1.warp.tif

The zero data area should be the same in both result files, but it is not. In the burnded warp.tif it is correct and in the directly warped it is not.

In stable 1.6.0 the cutline was shifted but polygon was at least visible. When I use latest svn-trunk the result is just a black window.

To test this bug I am attaching the related shapefile. The input raster (13218_1.KAP) is possible to download here:  http://www.charts.noaa.gov/RNCs/13218.zip.

All shapefiles and raster files has the same srs (epsg:4326) in this example. I intend to use reprojection during the cutline warping from gcps in 4326 into 900913 (then the shapefile will be in 900913 too). I hope I can use VRT with cutline - because I am going to read it just once in the native resolution.

[klokan]

Maybe better test for cutline and GCP reference is this file:  http://www.charts.noaa.gov/RNCs/14500.zip

which is deformed by GCPs polynomial transformation more. I am attaching also the shape file with cutline for testing...

[warmerdam]

I have done a substantial analysis.

First, with the 14500 file the BSB driver was returning an approximate geotransform even though GDALGetGeoTransform() returned FALSE indicating it was not valid. Gdalwarp was using it anyways giving an approximate affine transformation for the cutline to image coordinates.

The 14500 file is in polyconic, but this is not recognised by the driver so the plain WGS84 GCPs are utilized. Mercator files (like 13218) were working closer to correct because the affine transformation was nearly the same as a higher order GCP based transformation.

I have modified gdalwarp in trunk (r16589) to support transforming the cutline to image coordinates using essentially the same algorithm used for transforming the image. This includes the ability to use GCPs, and even to reproject if the cutline projection differs from the image projection.

With these changes results with the 13218 image are now very good due to it being mercator (rectangles remain rectangular in mercator).

With the polyconic image (14500) the results are still not that great, primarily because the GCP transformer does not give a good approximation of polyconic reprojection. However, it is better than it was.

The next step would be to implement better projection recognition for Polyconic in the BSB driver, and also to reduce the requirement for precision in the GCPs to GeoTransform? converter as used by the BSB driver which is preventing even the Mercator dataset from being treated as a "plain" mercator image.

[hamish]

see also #3409, GDAL is setting +lat_ts wrong for Mercator charts. It's probably close enough not to matter for harbour scale charts, but rather important for 1:10million scale ocean charts.

Hamish