GDAL 3.9.2 - gdalwarp raster dem degradation with VRT

[lezan]

What is the bug?

I am experience a raster degradation quality after warping a DEM file to VRT. I have the same output if I use gdalwarp from QGIS but also from prompt. With raster degradation quality I mean vertical and horizontal bands

tiff

vrt

I am warping from EPSG:4269 to a custom projection defined like that +proj=poly +lat_0=0 +lon_0=-118.375 +x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs +type=crs

Steps to reproduce the issue

1. Download the DEM from here: https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/historical/n37w119/USGS_13_n37w119_20240207.tif
2. Warp the dem with this command for the tif gdalwarp -overwrite -s_srs EPSG:4269 -t_srs "+proj=poly +lat_0=0 +lon_0=-118.375 +x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs +type=crs" -r near -of GTiff input.tif output.tif
3. Warp the dem for the VRT gdalwarp -overwrite -s_srs EPSG:4269 -t_srs "+proj=poly +lat_0=0 +lon_0=-118.375 +x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs +type=crs" -r near -of VRT input.tif output.vrt

Versions and provenance

OS: Windows 11. GDAL: 3.9.2. GDAL installation come OSGeo4W.

Additional context

No response

[jratike80]

How did you materialize (render) the vrt into image? The vrt itself is virtual, just a recipe about how the image should be processed.

I let GDAL to materialize the vrt and I can see no artifacts. gdal_translate output.vrt output_test.tif

It feels to me that the issue is in the viewer that you use for showing the vrt. What software it is?

[lezan]

Hello @jratike80 I'm sorry I didn't specify that. I imported the VRT into QGIS and then switch the render type of layer properties from Singleband gray to Hillshade. I leave all the values to default to avoid issues. Once imported the VRT, right click on the layer -> Properties -> Symbology -> Band Rendering -> Render Type.

Oh that's interesting. So from your side. once you translate the VRT to TIFF, this TIFF is the same (quality) of straight output the TIFF (commands I wrote in the first message)?

[jratike80]

Thank you for providing this essential info. I can see the same patchwork in output.tif as well when I zoom closer. This is output.tif zoomed in and warped with QGIS into some other CRS.

At this scale the .tif and .vrt looks essentially the same on my screen. Not too good, I agree, but I am not at all sure if there is an issue in GDAL and VRT. Maybe the difference comes from how QGIS is performing subsampling with hillshade with tiff vs. vrt source.

Notice that QGIS is using nearest neighbor resampling for hillshade by default. Play with the settings, which may be hidden because they are lowest on the page:

[lezan]

I am also experience the same pathwork if I warp to VRT and then output to TIFF. If I go straight to TIFF, the result is quite different, the same I shared in the first post.

I can see the same patchwork in output.tif as well when I zoom closer.

You cannot zoom to much because the DEM resolution is 10m, so I do not expect to get a good looking hillshading further 1:50000. The images attached in the first post were to the same scale, approx 1:80000.

Notice that QGIS is using nearest neighbor resampling for hillshade by default. Play with the settings, which may be hidden because they are lowest on the page

Yeah I know, but the results I shared were with the same options and at the same scale. I also used the TIFF in Blender and another 3D software to visualize the hillshading, and the VRT->TIFF is worst than the TIFF. I use this process daily and that's the first time I see something like that (different output VRT / TIFF)

I saw the same error also with 3.8.5 but before to open the issue I upgraded the latest version, but the output I got was the same.

I will try to deep into and see if I can find more info. For the moment I am experience the same issue also with, for example, EPSG:26744 and not only my custom projection.

By the way thanks for your support!

[jratike80]

Notice that GDAL offers the access to the data and QGIS hillshade rendered applies the style on-the-fly and can have its own issues.

This is 1:10000 scale rendering from QGIS when the VRT applies cubic resampling

gdalwarp -overwrite -s_srs EPSG:4269 -t_srs "+proj=poly +lat_0=0 +lon_0=-118.375 +x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs +type=crs" -r cubic -of VRT USGS_13_n37w119_20240207.tif output_cubic.vrt

and the QGIS hillshade rendered is also configured to use cubic resampling.

I see hardly any difference in the QGIS hillshade rendering with the default settings between the image warped directly into TIFF and the image warped first into VRT and then into TIFF. Anyway, the NN resampling seems to be the biggest source of artifacts so I would first make new series of tests with other resampling methods.

I suppose that you know what VRT is https://gdal.org/en/latest/drivers/raster/vrt.html and the meaning of "virtual" in its name.

[rouault]

Experimenting a bit with QGIS myself, I also get those weird artifacts at some downsampling levels when using the materialized warped dataset.

One major difference is that when using the warped VRT, given that the source GeoTIFF is a COG with several overview levels, the warped VRT also exposes virtual overviews, that can use the source overviews of the COG, depending on the subsampling factor. When gdalwarp'ing directly to a GeoTIFF, the full resolution dataset is used. To be fully honest, there might be minor differences when reading on-the-fly from a warped VRT and from a materialized warped TIFF, given how the on-the-fly reprojection works, but they lead to neglectable difference. The main issue as raised by @jratike80 is the use of nearest neighbour resampling at the GDAL and QGIS level. This is not appropriate for good hillshading result.

Amusingly, if you gdalwarp -r cubic to a VRT and a GeoTIFF, and open them in QGIS with hillshading, by default the result will be better for the VRT dataset. The reason is that it will natively apply cubic resampling for all requests that QGIS forward to it. Whereas for the TIFF, if you let QGIS default resampling settings, it will request with nearest neighbour downsampling requests to the TIFF, thus degrading its quality. As soon as you also ask QGIS to apply cubic resampling, you get identical nice results.

All in all, I don't think this is a GDAL bug