GeographicLib: Geoid height (original) (raw)

Back to GeographicLib 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:

EGM84: https://earth-info.nga.mil/index.php?dir=wgs84&action=wgs84#tab_egm84
EGM96: https://earth-info.nga.mil/index.php?dir=wgs84&action=wgs84#tab_egm96
EGM2008: https://earth-info.nga.mil/index.php?dir=wgs84&action=wgs84#tab_egm2008

Geoid offers a uniform way to handle all 3 geoids models at a variety of grid resolutions. (In contrast, the software tools that NGA offers are different for each geoid model, 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 format of the geoid data files
Interpolating the geoid data
Caching the geoid data
Test data for geoids

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)	Download Links (size, MB)
tar file	Windows installer	zip file
egm84-30	EGM84	30'	0.6	link (0.5)	link (0.8)	link (0.5)
egm84-15	EGM84	15'	2.1	link (1.5)	link (1.9)	link (2.0)
egm96-15	EGM96	15'	2.1	link (1.5)	link (1.9)	link (2.0)
egm96-5	EGM96	5'	19	link (11)	link (11)	link (17)
egm2008-5	EGM2008	5'	19	link (11)	link (11)	link (17)
egm2008-2_5	EGM2008	2.5'	75	link (35)	link (33)	link (65)
egm2008-1	EGM2008	1'	470	link (170)	link (130)	link (390)

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;
comments can be added to the file header;
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 https://earth-info.nga.mil/index.php?dir=wgs84&action=wgs84#tab_egm96

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 <math xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>l</mi><mi>o</mi><mi>n</mi></mrow><annotation encoding="application/x-tex">lon </annotation></semantics></math>lonlat $res | awk '{ lon = <math xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mn>1</mn><mo separator="true">;</mo><mi>l</mi><mi>a</mi><mi>t</mi><mo>=</mo></mrow><annotation encoding="application/x-tex">1; lat = </annotation></semantics></math>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)) ((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
TEMPeTEMP eTEMPelon-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
–a_offset -108 -a_scale 0.003
egmyy-g.pgm egmyy-g.tif

(The -a_offset and -a_scale options were introduced in gdal_translate version 2.3.) 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.

The Geoid class only handles world-wide geoid models. The pgm provides geoid height postings on grid of points with uniform spacing in latitude (row) and longitude (column). If the dimensions of the pgm file are h × w, then latitude (resp. longitude) spacing is −180°/(h − 1) (resp. 360°/w), with the (0,0) pixel given the value at φ = 90°, λ = 0°. Models covering a limited area will need to be "inserted" into a world-wide grid before being accessible to the Geoid class.

Interpolating the geoid data

Geoid evaluates the geoid height using bilinear or cubic interpolation.

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_5	EGM2008	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_5	EGM2008	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

...

GeographicLib::Geoid g("egm96-5");

double lat, lon;

while (std::cin >> lat >> lon)

std::cout << g(lat, lon) << "\n";

...

Header for GeographicLib::Geoid class.

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

GeoidHeights.dat.gz

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 5 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 GeographicLib in other languages. Forward to Gravity models. Up to Contents.