[gdal-dev] vrt warping precision issues and tile offsets

[Even Rouault]

David,

Ah, shouldn't have looked at this, as with all sub-pixel issues, it took me 
hours to figure it out, but the good news is I managed to fix it:
https://github.com/OSGeo/gdal/pull/1132

> 
> I build a vrt
> <https://drive.google.com/file/d/1AsXk5KFF-tQNaFfWZsAQLz0wQW_lsoOA/view?usp
> =sharing> of the reference tiles: gdalbuildvrt -r cubic ref.vrt tile1.tif
> tile2.tif tile3.tif tile4.tif
> 
> For comparison, I create a mosaic of the reference tiles using gdalwarp:
> gdal_opt='-co COMPRESS=LZW -co TILED=YES -co BIGTIFF=IF_SAFER'
> gdalwarp $gdal_opt -r cubic -et 0  tile1.tif tile2.tif tile3.tif tile4.tif
> ref.tif

At that stage, you could have just compared ref.vrt and ref.tif, and you would 
have already seen that, except the top-left tile, there were shift in the 
areas covered by the 3 other ones (the issue was somehow subtle, in that, it 
occured only if doing RasterIO requests that go beyond the tile extent, as 
done with gdal_translate/gdalwarp when operating on the whole raster, but if 
you picked up only a single pixel, or smaller subwindows, the result was 
correct)

You can use gdal_calc -A x.tif -B y.tif --outfile=diff.tif --calc "A-B"
to compute differences.

The issue was indeed in reading the VRT (indepedantly of the later warping you 
did it). With my fix, ref.vrt and ref.tif are now identical except a few 
pixels at tile boundaries and raster edges due to subtle difference in how VRT 
compositing and warp compositing deal with edge effects, but there is no 
longer a systematic offset.

> Is there something equivalent to -et 0 for gdalbuildvrt?  

Nope, modulo bugs, VRT reading is supposed to be 'perfect'

> The vrt DstRect offsets are suspicious and the output vrt
> dimensions don’t match the sum of input tile dimensions. 

That's expected with your data. If you look at the georeferencing of your 
tiles, they have the same resolution, but the shift between their origin is 
not a multiple of the resolution, and they have a few pixels overlap, so it is 
expected that gdalbuildvrt keeps that overlap and there are non-integer 
coordinates (hence the need for -r cubic with gdalbuildvrt. Otherwise if you 
had perfect mosaicing, it would be useless)

> Also, does the -r cubic option for gdalbuildvrt matter if I specify a
> different resampling algorithm during a later warp operation?

Yes, the -r cubic in gdalbuildvrt matters, and is independant of the 
resampling algorithm of a later warp operation. When you cubic-warp a cubic-
VRT, you have in fact two cubic interpolations occuring.

I would not that in your workflow (not sure to have completely understood it), 
you could directly do gdalwarp of your 4 ref tiles and reproject in the same 
step: this would avoid those two cubic interpolation stages that introduce 
some loss.

Even

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