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

[David Shean]

Hi list,
I am seeing some sub-pixel inter-tile offsets in the output of warped vrt products.  I am a big fan of vrt and use regularly, but this is the first time that I’ve noticed this issue.  I am hoping that there is a simple workaround.

I’ve isolated a simple test case from a larger project.  I’m running GDAL 2.2.4.

Test case includes:
-A single “source” DEM with UTM 43N projection (EPSG:32643), 30 m grid
-A set of 4 adjacent “reference” DEM tiles in a regional Albers equal-area projection ('+proj=aea +lat_1=25 +lat_2=47 +lat_0=36 +lon_0=85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs ‘), also 30 m grid

I need to compute elevation difference between the source and reference.  The source DEM extent spans tile boundaries of the reference DEM.

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

I then resample both ref.vrt and ref.tif to match the resolution/extent/srs of the original source DEM and compute an elevation difference.  In practice, this is done on-the-fly using https://github.com/dshean/demcoreg/blob/master/demcoreg/compute_diff.py <https://github.com/dshean/demcoreg/blob/master/demcoreg/compute_diff.py>.  For the purposes of illustration, and to ensure that these external tools are not the issue, I use gdalwarp for this resampling step:
src_extent='261016.895 3813317.750 335896.895 3885107.750'
gdalwarp $gdal_opt -r cubic -et 0 -tr 30 30  -t_srs EPSG:32643 -te $src_extent ref.tif reftif_warp.tif
gdalwarp $gdal_opt -r cubic -et 0 -tr 30 30  -t_srs EPSG:32643 -te $src_extent ref.vrt refvrt_warp.tif

The (reftif_warp.tif minus source) DEM difference map <https://drive.google.com/file/d/1VoB06LMArRmPKYxVyNkk5AHWhbNf69Kn/view?usp=sharing> looks as expected, and is consistent with the output that I obtain computing difference maps between the source DEM and each individual reference tile.

The (refvrt_warp.tif minus source) DEM difference map <https://drive.google.com/file/d/1PZmSTng8yQdC-1D0EBAVFdsTjMiQWL7C/view?usp=sharing> shows evidence of subtle horizontal shifts in the reference tiles, which introduces spatially variable vertical error in the resulting elevation differences.  I ran an image correlator on the outputs, and indeed, I see sub-pixel shifts <https://drive.google.com/file/d/1vzoUfiF9fPKO-AmzmIYOFwVNl8Y8Kom1/view?usp=sharing> between the original tiles with variable direction/magnitude (up to ~0.3 px).  The original tile boundaries are easy to identify.  

After various tests, I’ve concluded that warping the vrt is the problem.  Is this a known issue?  Is there something equivalent to -et 0 for gdalbuildvrt?  The vrt DstRect offsets are suspicious and the output vrt dimensions don’t match the sum of input tile dimensions.  Using -tap during gdalbuildvrt doesn’t help.  Maybe an issue related to my aea projection?  

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

Sample data are here: https://drive.google.com/drive/folders/15hNlUoUkCeYbgZCyfMzymmCP1__4wVgO?usp=sharing <https://drive.google.com/drive/folders/15hNlUoUkCeYbgZCyfMzymmCP1__4wVgO?usp=sharing>

Thanks for any suggestions.
-David

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20181205/928c4dc7/attachment-0001.html>