GeographicLib  1.46
Geoid height
Back to Implementations in other languages. Forward to Gravity models. Up to Contents.

The gravitational equipotential surface approximately coinciding with mean sea level is called the geoid. The Geoid class and the GeoidEval utility evaluate the height of the geoid above the WGS84 ellipsoid. This can be used to convert heights above mean sea level to heights above the WGS84 ellipsoid. Because the normal to the ellipsoid differs from the normal to the geoid (the direction of a plumb line) there is a slight ambiguity in the measurement of heights; however for heights up to 10 km this ambiguity is only 1 mm.

The geoid is usually approximated by an "earth gravity model" (EGM). The models published by the NGA are:

Geoid offers a uniform way to handle all 3 geoids at a variety of grid resolutions. (In contrast, the software tools that NGA offers are different for each geoid, and the interpolation programs are different for each grid resolution. In addition these tools are written in Fortran with is nowadays difficult to integrate with other software.)

The height of the geoid above the ellipsoid, N, is sometimes called the geoid undulation. It can be used to convert a height above the ellipsoid, h, to the corresponding height above the geoid (the orthometric height, roughly the height above mean sea level), H, using the relations

h = N + H;   H = −N + h.

Unlike other components of GeographicLib, there is a appreciable error in the results obtained (at best, the RMS error is 1 mm). However the class provides methods to report the maximum and RMS errors in the results. The class also returns the gradient of the geoid. This can be used to estimate the direction of a plumb line relative to the WGS84 ellipsoid.

The GravityModel class calculates geoid heights using the underlying gravity model. This is slower then Geoid but considerably more accurate. This class also can accurately compute all the components of the acceleration due to gravity (and hence the direction of plumb line).

Go to

# Installing the geoid datasets

The geoid heights are computed using interpolation into a rectangular grid. The grids are read from data files which have been are computed using the NGA synthesis programs in the case of the EGM84 and EGM96 models and using the NGA binary gridded data files in the case of EGM2008. These data files are available for download:

Available geoid data files
name geoid grid size
(MB)
tar file Windows
installer
zip file
egm84-30 EGM84
30'
0.6
egm84-15 EGM84
15'
2.1
egm96-15 EGM96
15'
2.1
egm96-5 EGM96
5'
19
egm2008-5 EGM2008
5'
19
egm2008-2_5 EGM2008
2.5'
75
egm2008-1 EGM2008
1'
470

The "size" column is the size of the uncompressed data and it also gives the memory requirements for caching the entire dataset using the Geoid::CacheAll method. At a minimum you should install egm96-5 and either egm2008-1 or egm2008-2_5. Many applications use the EGM96 geoid, however the use of EGM2008 is growing. (EGM84 is rarely used now.)

For Linux and Unix systems, GeographicLib provides a shell script geographiclib-get-geoids (typically installed in /usr/local/sbin) which automates the process of downloading and installing the geoid data. For example

   geographiclib-get-geoids best # to install egm84-15, egm96-5, egm2008-1
geographiclib-get-geoids -h   # for help


This script should be run as a user with write access to the installation directory, which is typically /usr/local/share/GeographicLib (this can be overridden with the -p flag), and the data will then be placed in the "geoids" subdirectory.

Windows users should download and run the Windows installers. These will prompt for an installation directory with the default being

   C:/ProgramData/GeographicLib


(which you probably should not change) and the data is installed in the "geoids" sub-directory. (The second directory name is an alternate name that Windows 7 uses for the "Application Data" directory.)

Otherwise download either the tar.bz2 file or the zip file (they have the same contents). If possible use the tar.bz2 files, since these are compressed about 2 times better than the zip file. To unpack these, run, for example

   mkdir -p /usr/local/share/GeographicLib
tar xofjC egm96-5.tar.bz2 /usr/local/share/GeographicLib
tar xofjC egm2008-2_5.tar.bz2 /usr/local/share/GeographicLib
etc.


and, again, the data will be placed in the "geoids" subdirectory.

However you install the geoid data, all the datasets should be installed in the same directory. Geoid and GeoidEval uses a compile time default to locate the datasets. This is

• /usr/local/share/GeographicLib/geoids, for non-Windows systems
• C:/ProgramData/GeographicLib/geoids, for Windows systems

consistent with the examples above. This may be overridden at run-time by defining the GEOGRAPHICLIB_GEOID_PATH or the GEOGRAPHICLIB_DATA environment variables; see Geoid::DefaultGeoidPath() for details. Finally, the path may be set using the optional second argument to the Geoid constructor or with the "-d" flag to GeoidEval. Supplying the "-h" flag to GeoidEval reports the default geoid path for that utility. The "-v" flag causes GeoidEval to report the full path name of the data file it uses.

# The format of the geoid data files

The gridded data used by the Geoid class is stored in 16-bit PGM files. Thus the data for egm96-5 might be stored in the file

• /usr/local/share/GeographicLib/geoids/egm96-5.pgm

PGM simple graphic format with the following properties

• it is well documented here;
• there are no separate "little-endian" and "big-endian" versions;
• it uses 1 or 2 bytes per pixel;
• pixels can be randomly accessed;
• it is sufficiently simple that it can be easily read without using the libnetpbm library (and thus we avoid adding a software dependency to GeographicLib).

The major drawback of this format is that since there are only 65535 possible pixel values, the height must be quantized to 3 mm. However, the resulting quantization error (up to 1.5 mm) is typically smaller than the linear interpolation errors. The comments in the header for egm96-5 are

   # Geoid file in PGM format for the GeographicLib::Geoid class
# Description WGS84 EGM96, 5-minute grid
# URL http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html
# DateTime 2009-08-29 18:45:03
# MaxBilinearError 0.140
# RMSBilinearError 0.005
# MaxCubicError 0.003
# RMSCubicError 0.001
# Offset -108
# Scale 0.003
# Origin 90N 0E
# AREA_OR_POINT Point
# Vertical_Datum WGS84


Of these lines, the Scale and Offset lines are required and define the conversion from pixel value to height (in meters) using height = offset + scale pixel. The Geoid constructor also reads the Description, DateTime, and error lines (if present) and stores the resulting data so that it can be returned by Geoid::Description, Geoid::DateTime, Geoid::MaxError, and Geoid::RMSError methods. The other lines serve as additional documentation but are not used by this class. Accompanying egm96-5.pgm (and similarly with the other geoid data files) are two files egm96-5.wld and egm96-5.pgm.aux.xml. The first is an ESRI "world" file and the second supplies complete projection metadata for use by GDAL. Neither of these files is read by Geoid.

You can use gdal_translate to convert the data files to a standard GeoTiff, e.g., with

   gdal_translate -ot Float32 -scale 0 65000 -108 87 egm96-5.pgm egm96-5.tif


The arguments to -scale here are specific to the Offset and Scale parameters used in the pgm file (note 65000 * 0.003 − 108 = 87). You can check these by running GeoidEval with the "-v" option.

Here is a sample script which uses GDAL to create a 1-degree squared grid of geoid heights at 3" resolution (matching DTED1) by bilinear interpolation.

   #! /bin/sh
lat=37
lon=067
res=3                           # resolution in seconds
TEMP=mktemp junkXXXXXXXXXX    # temporary file for GDAL
gdalwarp -q -te echo $lon$lat $res | awk '{ lon =$1; lat = $2; res =$3;
printf "%.14f %.14f %.14f %.14f",
lon  -0.5*res/3600, lat  -0.5*res/3600,
lon+1+0.5*res/3600, lat+1+0.5*res/3600;
}' -ts $((3600/res+1))$((3600/res+1)) -r bilinear egm96-5.tif $TEMP gdal_translate -quiet \ -mo AREA_OR_POINT=Point \ -mo Description="WGS84 EGM96,$res-second grid" \
-mo Vertical_Datum=WGS84 \
-mo Tie_Point_Location=pixel_corner \
$TEMP e$lon-n$lat.tif rm -f$TEMP


Because the pgm files are uncompressed, they can take up a lot of disk space. Some compressed formats compress in tiles and so might be compatible with the requirement that the data can be randomly accessed. In particular gdal_translate can be used to convert the pgm files to compressed tiff files with

gdal_translate -co COMPRESS=LZW -co PREDICTOR=2 \
-co TILED=YES -co BLOCKXSIZE=256 -co BLOCKYSIZE=256 \
egmyy-g.pgm egmyy-g.tif


The resulting files sizes are

                  pgm      tif
egm84-30      0.6 MB   0.5 MB
egm84-15      2.1 MB   1.4 MB
egm96-15      2.1 MB   1.5 MB
egm96-5        19 MB   8.5 MB
egm2008-5      19 MB   9.8 MB
egm2008-2_5    75 MB    28 MB
egm2008-1     470 MB    97 MB


Currently, there are no plans for GeographicLib to support this compressed format.

# Interpolating the geoid data

Geoid evaluates the geoid height using bilinear or cubic interpolation. The gradient of the geoid height is obtained by differentiating the interpolated height and referencing the result to distance on the WGS84 ellipsoid.

WARNING: Although Geoid computes the gradient of the geoid height, this is subject to large quantization errors as a result of the way that the geoid data is stored. This is particularly acute for fine grids, at high latitudes, and for the easterly gradient. If you need to compute the direction of the acceleration due to gravity accurately, you should use GravityModel::Gravity.

The bilinear interpolation is based on the values at the 4 corners of the enclosing cell. The interpolated height is a continuous function of position; however the gradient has discontinuities are cell boundaries. The quantization of the data files exacerbates the errors in the gradients.

The cubic interpolation is a least-squares fit to the values on a 12-point stencil with weights as follows:

   . 1 1 .
1 2 2 1
1 2 2 1
. 1 1 .


The cubic is constrained to be independent of longitude when evaluating the height at one of the poles. Cubic interpolation is considerably more accurate than bilinear interpolation; however, in this implementation there are small discontinuities in the heights are cell boundaries. The over-constrained cubic fit slightly reduces the quantization errors on average.

The algorithm for the least squares fit is taken from, F. H. Lesh, Multi-dimensional least-squares polynomial curve fitting, CACM 2, 29–30 (1959). This algorithm is not part of Geoid; instead it is implemented as Maxima code which is used to precompute the matrices to convert the function values on the stencil into the coefficients from the cubic polynomial. This code is included as a comment in Geoid.cpp.

The interpolation methods are quick and give good accuracy. Here is a summary of the combined quantization and interpolation errors for the heights.

Interpolation and quantization errors for geoid heights
name geoid grid
bilinear error
cubic error
max
rms
max
rms
egm84-30 EGM84
30'
1.546 m
70 mm
0.274 m
14 mm
egm84-15 EGM84
15'
0.413 m
18 mm
0.021 m
1.2 mm
egm96-15 EGM96
15'
1.152 m
40 mm
0.169 m
7.0 mm
egm96-5 EGM96
5'
0.140 m
4.6 mm
0.0032 m
0.7 mm
egm2008-5 EGM2008
5'
0.478 m
12 mm
0.294 m
4.5 mm
egm2008-2_5EGM2008
2.5'
0.135 m
3.2 mm
0.031 m
0.8 mm
egm2008-1 EGM2008
1'
0.025 m
0.8 mm
0.0022 m
0.7 mm

The errors are with respect to the specific NGA earth gravity model (not to any "real" geoid). The RMS values are obtained using 5 million uniformly distributed random points. The maximum values are obtained by evaluating the errors using a different grid with points at the centers of the original grid. (The RMS difference between EGM96 and EGM2008 is about 0.5 m. The RMS difference between EGM84 and EGM96 is about 1.5 m.)

The errors in the table above include the quantization errors that arise because the height data that Geoid uses are quantized to 3 mm. If, instead, Geoid were to use data files without such quantization artifacts, the overall error would be reduced but only modestly as shown in the following table, where only the changed rows are included and where the changed entries are given in bold:

Interpolation (only!) errors for geoid heights
name geoid grid
bilinear error
cubic error
max
rms
max
rms
egm96-5 EGM96
5'
0.140 m
4.6 mm
0.0026 m
0.1 mm
egm2008-2_5EGM2008
2.5'
0.135 m
3.2 mm
0.031 m
0.4 mm
egm2008-1 EGM2008
1'
0.025 m
0.6 mm
0.0010 m
0.011 mm

# Caching the geoid data

A simple way of calling Geoid is as follows

#include <iostream>
...
GeographicLib::Geoid g("egm96-5");
double lat, lon;
while (std::cin >> lat >> lon)
std::cout << g(lat, lon) << "\n";
...

The first call to g(lat, lon) causes the data for the stencil points (4 points for bilinear interpolation and 12 for cubic interpolation) to be read and the interpolated value returned. A simple 0th-order caching scheme is provided by default, so that, if the second and subsequent points falls within the same grid cell, the data values are not reread from the file. This is adequate for most interactive use and also minimizes disk accesses for the case when a continuous track is being followed.

If a large quantity of geoid calculations are needed, the calculation can be sped up by preloading the data for a rectangular block with Geoid::CacheArea or the entire dataset with Geoid::CacheAll. If the requested points lie within the cached area, the cached data values are used; otherwise the data is read from disk as before. Caching all the data is a reasonable choice for the 5' grids and coarser. Caching all the data for the 1' grid will require 0.5 GB of RAM and should only be used on systems with sufficient memory.

The use of caching does not affect the values returned. Because of the caching and the random file access, this class is not normally thread safe; i.e., a single instantiation cannot be safely used by multiple threads. If multiple threads need to calculate geoid heights, there are two alternatives:

• they should all construct thread-local instantiations.
• Geoid should be constructed with threadsafe = true. This causes all the data to be read at the time of construction (and if this fails, an exception is thrown), the data file to be closed and the single-cell caching to be turned off. The resulting object may then be shared safely between threads.

# Test data for geoids

A test set for the geoid models is available at

This is about 11 MB (compressed). This test set consists of a set of 500000 geographic coordinates together with the corresponding geoid heights according to various earth gravity models.

Each line of the test set gives 6 space delimited numbers

• latitude (degrees, exact)
• longitude (degrees, exact)
• EGM84 height (meters, accurate to 1 μm)
• EGM96 height (meters, accurate to 1 μm)
• EGM2008 height (meters, accurate to 1 μm)

The latitude and longitude are all multiples of 10−6 deg and should be regarded as exact. The geoid heights are computed using the harmonic NGA synthesis programs, where the programs were compiled (with gfortran) and run under Linux. In the case of the EGM2008 program, a SAVE statement needed to be added to subroutine latf, in order for the program to compile correctly with a stack-based compiler. Similarly the EGM96 code requires a SAVE statement in subroutine legfdn.

Back to Implementations in other languages. Forward to Gravity models. Up to Contents.