12

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.

tinlyx
  • 11,057
  • 18
  • 71
  • 119
Josée-Anne Langlois
  • 121
  • 1
  • 3

5 Answers5

11

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

Vitaly Olegovitch
  • 109
  • 5
B. Raskob
  • 111
  • 1
  • 2
  • 2
    I'm using GDAL 3.* and Windows 10 with Windows PowerShell and I had to add a ./ to the egm96_15.gtx file so gdal would find it. Without specifying this local directory I was getting a "Cannot open" file error. 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 Otherwise, this looks like it worked well for me while other top Google hits for this conversion didn't work. – Brandon Jun 16 '20 at 21:04
4

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.)

Rob Skelly
  • 3,101
  • 15
  • 14
  • I just saw a ticket on gdal for the release 2.1.0 https://trac.osgeo.org/gdal/ticket/6084.. I'm not sure I understand it correctly, but it seems like gdalwarp should be able to do vertical datum conversion with this release? Do you understand the same thing? – Josée-Anne Langlois Apr 28 '16 at 19:53
  • I believe that only relates to orthorectification using RPCs. – Rob Skelly Apr 28 '16 at 20:19
  • I am interested in this exercise. See https://lists.osgeo.org/pipermail/gdal-dev/2017-November/047759.html. – Nikos Alexandris Nov 21 '17 at 12:02
  • 1
    @NikosAlexandris Did you see the follow-ups? gdalwarp in versions >=2.2 can now do vertical datum transformations, so long as the grid shift file is available to proj4/GDAL. – Rob Skelly Nov 21 '17 at 18:28
1

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/).

Shashank Bhushan
  • 43
  • 6
0

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.

Charlie
  • 31
  • 4
0

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.

Brent
  • 11
  • 2