2

I am currently trying to prepare an SRTM subset including height conversion from EGM96 to WGS84 using gdalwarp. I encountered some strange shifts in part of the image.
Here is a reproducible example:
Prepare the necessary file for the height conversion:

cd /usr/share/proj
sudo wget https://download.osgeo.org/proj/vdatum/egm96_15/egm96_15.gtx
sudo chmod 644 egm96_15.gtx

Download two SRTM 3Sec tiles and combine them in a VRT:

wget http://srtm.csi.cgiar.org/wp-content/uploads/files/srtm_5x5/TIFF/srtm_39_03.zip
wget http://srtm.csi.cgiar.org/wp-content/uploads/files/srtm_5x5/TIFF/srtm_38_03.zip
gdalbuildvrt dem.vrt /vsizip/srtm_39_03.zip/srtm_39_03.tif \
/vsizip/srtm_38_03.zip/srtm_38_03.tif \
-te 9.497235 45.432812 13.183811 47.336189

Use gdalwarp to create two new files, one that is just a simple conversion of the original file content and one that contains heights converted to WGS84 ellipsoid:

gdalwarp -r bilinear -of GTiff dem.vrt dem.tif
gdalwarp -r bilinear -of GTiff -s_srs EPSG:4326+5773 -t_srs EPSG:4326 \
dem.vrt dem_ellp.tif

Compute their difference as GTiff and JPEG to get the EGM96-WGS84 difference:

gdal_calc.py -A dem.tif -B dem_ellp.tif --outfile=dem_ellp_diff.tif --calc="A-B"
gdal_translate -of JPEG -scale -outsize 50% 50% \
dem_ellp_diff.tif dem_ellp_diff.jpg

The output looks like this: enter image description here

To the West and East the image looks like expected with simple gradient lines. In the center however it looks like a DEM aspect file. On closer inspection of the two files, one can see a pixel shift to the East that is introduced by gdalwarp in the second (ellipsoid height) file.
Is this problem known to someone? It does not appear when the DEM is not subsetted during creating the vrt file, i.e. calling gdalbuildvrt without the te argument: enter image description here

Also, the images look the same when using a different spatial CRS e.g. UTM via
-t_srs EPSG:32632

The problem is also not present when alternatively setting the extent via gdalwarp instead of gdalbuildvrt:

gdalbuildvrt dem2.vrt /vsizip/srtm_39_03.zip/srtm_39_03.tif \
/vsizip/srtm_38_03.zip/srtm_38_03.tif
gdalwarp -r bilinear -of GTiff -te 9.497235 45.432812 13.183811 47.336189 \
dem2.vrt dem2.tif
gdalwarp -r bilinear -of GTiff -s_srs EPSG:4326+5773 -t_srs EPSG:4326 \
-te 9.497235 45.432812 13.183811 47.336189 dem2.vrt dem2_ellp.tif
gdal_calc.py -A dem2.tif -B dem2_ellp.tif --outfile=dem2_ellp_diff.tif \
--calc="A-B"
gdal_translate -of JPEG -scale -outsize 50% 50% \ 
dem2_ellp_diff.tif dem2_ellp_diff.jpg

enter image description here GDAL v2.2.3, Ubuntu 18.04, libproj v4.9.3-2

  • I wish all questions were as well constructed as this one. – user2856 Mar 5 at 20:38
  • thanks @user2856! – John Truckenbrodt Mar 6 at 8:35
1

I just tried this again with GDAL version 2.4.0 from ubuntugis unstable.
Here it works perfectly fine and was apparently fixed somewhere between version 2.2.3 and 2.4.0.

  • In the release notes of GDAL 2.4.0 it reads "fix source window computation that caused sub-pixel shift with non-nearest resampling", which might be the relevant bug fix – John Truckenbrodt Mar 5 at 18:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.