Overview problems when working with warped VRTs

[mojodna]

When using gdalwarp to create a warped VRT whilst changing the resolution, incorrect overview levels are used.

E.g. Adding overviews to ​http://oin-astrodigital.s3.amazonaws.com/LC81600422014266LGN00_bands_432.TIF (3857-projected image with resolution ~= z12) and creating a warped VRT:

$ gdaladdo LC81600422014266LGN00_bands_432.TIF 2 4 8 16 32
0...10...20...30...40...50...60...70...80...90...100 - done.
$ gdalwarp LC81600422014266LGN00_bands_432.TIF \
   z20.vrt \
   -t_srs EPSG:3857 \
   -overwrite \
   -of VRT \
   -te -20037508.34 -20037508.34 20037508.34 20037508.34 \
   -ts 268435456 268435456
Creating output file that is 268435456P x 268435456L.
Processing input file LC81600422014266LGN00_bands_432.TIF.
Using internal nodata values (e.g. 0) for image LC81600422014266LGN00_bands_432.TIF.
Copying nodata values from source LC81600422014266LGN00_bands_432.TIF to destination z20.vrt.

(-20037508.34 -20037508.34 20037508.34 20037508.34 is the world extent in web mercator, 268435456 is the pixel width/height at zoom 20 (assuming 256x256 tiles))

Relevant output from gdalinfo:

Band 1 Block=512x128 Type=Byte, ColorInterp=Red
  NoData Value=0
  Overviews: 134235179x134235179, 67117590x67117590, 33576246x33576246, 16788123x16788123, 8411513x8411513

This suggests (and I've confirmed its behavior via rasterio) that overviews are available from zooms 19-15 when in reality overviews are available from 11-7. The result of this is that the "zoom 15" overviews (index 4) are used when rendering data at zoom 12 (the approximate native resolution of the source image), producing extremely pixelated output.

I would expect that the warped VRT would use "dummy" overviews until the resolution of the available overviews begins to approximately match.

(Context: I'm trying to create source files that support reading (decimated) windows specified in zoom 20-level pixel coordinates for consistency, even when imagery is available (in 3857 or otherwise) at lower resolutions.)

(Further aside: the reported blocking (512x128) doesn't match the source file (7691x1) even after re-tiling to 256x256. I'm assuming it's reporting a synthetic value, but why isn't it square?)

[Even Rouault]

Hum, what you are doing is quite particular. The issue is that you use gdalwarp to do major oversampling on the most precise resolution, but you wouldn't want overviews of the underlying raster to be oversampled the same way, but used only when at compatible resolutions. I can understand the use case, but current behaviour of implicit overviews is not completely inconsistant and more common.

What you could do instead is to use explicit overviews by editing the VRT. See ​http://gdal.org/gdal_vrttut.html#gdal_vrttut_overviews and the <Overview> element. Approximate workflow would be to create a warped VRT for each of the overview of the TIF and then use that VRT as one <Overview> of the final VRT.

Note: if quality matters, you should consider using non-nearest neighbour resampling for overview building and warping.

(as far as 512x128 blocks, this is hard-coded in the driver and I've no idea why this isn't square.)

[mojodna]

Agreed, if I'm oversampling the image, my intention is also to oversample other things associated with the image, so this behavior does make sense.

I'll give the explicit overview approach a shot in the next few days, that has the added benefit of being able to leave the source images as-is.

I'm ignoring reported blocking (it could be a potential optimization, but not when the tiles are uniform), so that's no problem. I also stripped all of the extra params (creation options, resampling, etc.); I'm using cubic for the time-being.

Thank you!

[Even Rouault]

Closing as I don't think there's an action on GDAL side