When working with elevation data, sometimes you may discover that 2 datasets from different providers have very different elevation values for the same location. A common reason for this being each dataset being referenced to a different surface.
Orthometric (or geoid) elevation measures height of a point from mean-sea-level (MSL). Ground based surveys will produce orthometric elevations
Ellipsoidal elevation measures height of a point from a reference ellipsoid (such as WGS84). Elevations derived from space-borne platforms are usually ellipsoidal elevations.
Ellipsoids are a simple model of the earth’s surface and are defined using semi-major axis (a) and semi-minor axis (b). Geoid is a much more complex model as it attempts to model the gravitational force of the earth which is not uniform. There exists several models which define the geoid shape across the entire surface of the earth.
As both of these are mathematical models, it is possible to convert elevations from one reference to another. GDAL command-line tools are a free and easy way to do these conversions
Here I will show a real example of converting between 2 elevation datasets.
Cartosat-1 Digital Elevation Model (CartoDEM) is a National DEM developed by the Indian Space Research Organization (ISRO). It is derived from the Cartosat-1 stereo payload launched in May 2005. The Cartosat DEM is referenced to the WGS84 ellipsoid model (and hence have ellipsoidal elevation). The source data at 30m resolution can be downloaded from BHUVAN.
Shuttle Radar Topography Mission (SRTM) is an international project spearheaded by the U.S. National Geospatial-Intelligence Agency (NGA) and the U.S. National Aeronautics and Space Administration (NASA). During its 11-day mission, the space shuttle Endeavour orbited the Earth 16 times and captured Earth’s topography at 1 arc-second (30 meters) for over 80% of the Earth’s surface. An easy interface to download SRTM 30m data is the 30-Meter SRTM Tile Downloader. SRTM data is distributed with referenced to the EGM96 geoid model (and hence have orthometric elevations)
The GDAL command needed to convert the dataset is gdalwarp. We need to supply it with a source spatial reference system (s_srs) and a target spatial reference system (t_srs) values.
A previous version of this article used EPSG codes to specify a compound CRS, such as EPSG:4326+5773. Since the release of GDAL3/PROJ6, this syntax is no longer supported. The article has been updated with the new syntax.
How to run these commands?
Install GDAL utilities following these instructions for your platform.
One installed, open the shell, then
cd to the directory containing your files. For example, if your file is in C:\Users\ujaval\Downloads\, run the following command first
The commands in the sections below assume your filename is
cartosat.tif. Replace that with your filenames before running the commands.
Converting Between Vertical CRSs
Both the datasets have the same horizontal srs but a different vertical srs. So we can reference them using a proj strings as follows:
SRTM = WGS84 ellipsoid with EGM96 elevation
+proj=longlat +datum=WGS84 +no_defs +geoidgrids=egm96_15.gtx
Cartosat-1 = WGS84 ellipsoid with WGS84 elevation
+proj=longlat +datum=WGS84 +no_def
So now the conversion between the 2 is a matter of specifying the source and target srs.
Below is the command used for converting Cartosat-1 DEM to orthometric elevation.
gdalwarp -s_srs "+proj=longlat +datum=WGS84 +no_def" -t_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=egm96_15.gtx" cartosat.tif cartosat_orthometric.tif
If you load original cartosat.tif, cartosat_orthometric.tif and srtm.tif in QGIS, and inspect the elevation values, you can see the the cartosat_orthometric and srtm elevations match almost exactly now since they are referenced to the same geoid.
Above image contains a visualization of CartoDEM Version-3 R1 provided by National Remote Sensing Centre, ISRO, Government of India, Hyderabad, India.
GDAL using the PROJ library to carry out CRS transformations. PROJ uses pre-defined grids for datum transformations. Vertical datum transformations are defined using a GTX file. In the example above, when we did the transformation to EGM96 vertical CRS, it used the parameters supplied in the egm96_15.gtx file. Many such grid files are included when you download and install GDAL. But there are other grids which are very large and are not included by default. If you want to convert to a vertical datum whose grid files are not included by default, you need to download them separately and copy them to appropriate directory on your system. Learn more about PROJ Resource files.
One such vertical CRS is EGM2008 (EPSG:3855). Say you wanted to transform WGS84 heights to this new and more precise geoid model. You can run a command such as below
gdalwarp -s_srs "+proj=longlat +datum=WGS84 +no_def" -t_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=egm08_25.gtx" cartosat.tif cartosat_orthometric.tif
But if you run it, you may get an error such as below
ERROR 1: Cannot open egm08_25.gtx.
This is because the EGM2008 grid is very large and is not included in many GDAL installations. To fix this, you can download the grid file from http://download.osgeo.org/proj/vdatum/egm08_25/
Copy the egm08_25.gtx file to the PROJ resource directory on your computer. The location of this directory will depend on your platform and the installation method. Some common paths are as below
Windows: C:\OsGeo4W64\share\proj\ Mac: /Library/Frameworks/PROJ.framework/Resources/proj/ Linux: /usr/share/proj/
Once you copy the file, gdalwarp should be able to do the transformation correctly.
If this doesn’t work and GDAL can’t locate the grid file, you can also specify the full path to your grid file. The grid file can be either in the
gtx or a
tif format. For example, if you have a a German Grid (GCG2016) transform file located in C:\ drive, replace
+geoidgrids=C:\GCG2016.tif (Thanks to Christian Loeser for this tip)