11

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.

5 Answers 5

10

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

1
  • 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, 2020 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.)

4
  • I just saw a ticket on gdal for the release 2.1.0 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? Apr 28, 2016 at 19:53
  • I believe that only relates to orthorectification using RPCs.
    – Rob Skelly
    Apr 28, 2016 at 20:19
  • I am interested in this exercise. See lists.osgeo.org/pipermail/gdal-dev/2017-November/047759.html. Nov 21, 2017 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, 2017 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/).

0

We created a package specifically 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. Note we do some book-keeping for various DEMs related to how pixels are centered.

Check out the readme/notebooks/code for more information. Functions in the repo could be used on their own to remove the geoid for regional, high-resolution DEMs that are not specified.

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.

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.