Vertical Datum conversion of a raster with gdal

Question

Is there a tool available in gdal to do vertical datum conversion of a raster DEM? I'm looking for a single command line to execute the vertical datum conversion of a raster like the las2las of Liblas does for Point Cloud data.

I found out that gdalwarp were not doing this conversion even if the .gtx grid shift file exist in the PROJ_LIB folder and is correctly set in gdal.

I was able to do vertical datum conversion of a coordinate with gdaltransform, but I can't find a way to do the conversion for a complete raster file.

Answer 1

Convert from (EGM96) geoid vertical datum to (WGS84) ellipsoid vertical datum:

gdalwarp -s_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=egm96_15.gtx" -t_srs "+proj=longlat +datum=WGS84 +no_def" dem_in_egm96_geoid.tif dem_in_wgs84_ellipsoid.tif

You will need the gtx file containing vertical datum shifts from here: http://download.osgeo.org/proj/vdatum/egm96_15/

See this page for more info: https://proj.org/usage/transformation.html

You can test the result at a location using this site: https://www.unavco.org/software/geodetic-utilities/geoid-height-calculator/geoid-height-calculator.html

Answer 2

This answer isn't a single command, but I'll put it in to get the ball rolling. Use gdalwarp to resample the geoid grid, then gdal_calc.py to shift the original raster.

gdalwarp -s_srs epsg:4326 -t_srs epsg:26910 -r cubic -tr 10 10 -tap HT2_0.gtx HT2_0_resampled.tif
gdal_calc.py -A original.tif -B HT2_0_resampled.tif --calc="A+B" --outfile=shifted.tif

If original is a 10mx10m raster with grid origin at 0,0, UTM10N and contains orthometric heights, shifted.tif will contain ellipsoidal heights.

(Obviously, the geoid raster would be huge. One should crop it at the same time to speed things up/save space.)

Edit: GDAL >= 2.2 can now perform vertical datum transformations using the updated capabilities of Proj4. For example, with gdalwarp you'd use -s_srs and -t_srs with a vertical and horizontal SRID, like: gdalwarp ... -t_srs epsg:2956+3567 which gives you NAD83(CSRS) UTM12N with CGVD27 (if memory serves.)

Answer 3

See the direct function (dem_geoid) for geoid conversion available in the NASA Ames Stereo Pipeline(https://ti.arc.nasa.gov/tech/asr/groups/intelligent-robotics/ngt/stereo/).

Answer 4

As @ShashankBhushan pointed out, Ames Stereo Pipeline can be used for this task. As an example, DEM data that was originally in EGM96 can be converted to WGS84 ellipsoid using the following line in the terminal:

dem_geoid dem.tif --geoid egm96 --reverse-adjustment -o dem_wgs84

This was tested successfully using version 3.0.0.

Answer 5

We created an open source package for this purpose for a number of DEMs including the Copernicus 30 meter DEM: https://github.com/ACCESS-Cloud-Based-InSAR/dem-stitcher. An example to convert into ellipsoidal heights is:

bounds = [-119.085, 33.402, -118.984, 35.435]
X, p = stitch_dem(bounds,
                  dem_name='glo_30',
                  # This would remove geoid
                  dst_ellipsoidal_height=True,
                  dst_area_or_point='Area')
# X is an m x n numpy array
# p is a dictionary (or a rasterio profile) including relevant GIS metadata

We read geoid data from here and remove them (this is an addition operation). Note we do some book-keeping for various DEMs related to how pixels are centered as well.

Check out the readme and notebooks for more information. Hope this helps.