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

[David Shean]

Hi Even,
Many thanks for looking into this and pushing a fix so quickly.  Some follow up responses inline.
D

> On Dec 5, 2018, at 10:08 AM, Even Rouault <even.rouault at spatialys.com> wrote:
> 
> 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

Ha.  I hear you.  Things get messy down there, and I’ve spent too many hours working with sub-pixel issues.  Fortunately, most of the time GDAL performs as expected.

> 
>> 
>> 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)

Ah, yes.  That would have been simpler.  I incorrectly attributed to the second warp.

Your finding on subwindow requests within a single tile is consistent with what I’ve observed for other output (and why this issue went unnoticed for so long).

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

Yep. But I do enough raster calculations on sources with different dimensions that I don’t want to think about preparing new files with identical extent/res/srs to feed to gdal_calc.  

That compute_diff.py utility handles all of this on the fly with convenience options (intersection extent, mean resolution, etc) using MEM datasets and the underlying gdal.ReprojectImage functionality wrapped in the pygeotools/lib/warplib.py module.  Probably cleaner ways to do this MEM ds “matching” with latest rasterio functionality.

> 
> 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.

OK, I will test and check for any new edge artifacts.

> 
>> Is there something equivalent to -et 0 for gdalbuildvrt?  
> 
> Nope, modulo bugs, VRT reading is supposed to be ‘perfect'

:) Glad to hear it.

> 
>> 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)

I'm glad that my imperfect tiles helped reveal an underlying bug in vrt generation :).

OK.  These tiles are generated independently.  I precompute all tile extents: https://github.com/dshean/gmbtools/blob/master/gmbtools/dem_mosaic_validtiles.py#L144
then make separate calls to the Ames Stereo Pipeline dem_mosaic utility in parallel.  It’s a little clunky, but efficient for the available resources.

The origin, extent and dimensions of all tiles should be systematic at this point, always multiples of resolution.  For example:

dem_mosaic -o hma_mos_30m-tile-389-median.tif --tr 30.0 --t_projwin -1204970.0 -38010.0 -1104970.0 61990.0 …

Produces output tile with extent: -1204985.0 -38045.0 -1104935.0 62005.0

So looks like dem_mosaic is padding the tile, apparently nonuniformly.  Will look into this, thanks for the heads up.  

dem_mosaic also recomputes output tile extent based on valid data in input rasters, trimming nodata where appropriate, which can lead to non-systematic origin and tile dimensions.  While the vrt should account for this, it’s probably a better idea to create more “perfect” tiles, even if we have to store a bunch of compressed nodata pixels.  

> 
>> 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.

OK, this is good to know.  I agree it would be nice to limit to one resampling step, but not sure it will work for my needs.

If I understand correctly, I need the -r cubic in gdalbuildvrt to handle my imperfect tiles. 

My workflow requires that I create multiple vrts from large tiled datasets, each with different resolution, extent, srs.  I then need to extract the same projected extent from each vrt, with on-the-fly resampling to a common grid for subsequent analysis.  So I believe a second cubic warp is also required.

I was kind of hoping that GDAL might be able to recognize this and only perform a single warp under the hood, but I realize that is not trivial.

Thanks again for investigating, fixing and providing useful feedback.

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