/usr/share/doc/geographiclib/html/geoid.html is in geographiclib-tools 1.8-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 | <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<title>GeographicLib: Geoid height</title>
<link href="tabs.css" rel="stylesheet" type="text/css"/>
<link href="doxygen.css" rel="stylesheet" type="text/css"/>
</head>
<body>
<!-- Generated by Doxygen 1.7.1 -->
<div class="navigation" id="top">
<div class="tabs">
<ul class="tablist">
<li><a href="index.html"><span>Main Page</span></a></li>
<li class="current"><a href="pages.html"><span>Related Pages</span></a></li>
<li><a href="namespaces.html"><span>Namespaces</span></a></li>
<li><a href="annotated.html"><span>Classes</span></a></li>
<li><a href="files.html"><span>Files</span></a></li>
<li><a href="dirs.html"><span>Directories</span></a></li>
</ul>
</div>
</div>
<div class="header">
<div class="headertitle">
<h1>Geoid height </h1> </div>
</div>
<div class="contents">
<center> Back to <a class="el" href="geocentric.html">Geocentric coordinates</a>. Forward to <a class="el" href="utilities.html">Utility Programs</a>. Up to <a class="el" href="index.html#contents">Contents</a>. </center><p>The gravitational equipotential surface approximately coinciding with mean sea level is called the geoid. The <a class="el" href="classGeographicLib_1_1Geoid.html" title="Looking up the height of the geoid.">GeographicLib::Geoid</a> class evaluates 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 ellipoid. 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 10km this ambiguity is only 1mm.</p>
<p>The geoid is usually approximated by an "earth gravity model". The models published by the NGA are:</p>
<ul>
<li>EGM84: <a href="http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html">http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html</a></li>
<li>EGM96: <a href="http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html">http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html</a></li>
<li>EGM2008: <a href="http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008">http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008</a></li>
</ul>
<p><a class="el" href="classGeographicLib_1_1Geoid.html" title="Looking up the height of the geoid.">GeographicLib::Geoid</a> 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.)</p>
<p>Unlike other components of GeographicLib, there is a appreciable error in the results obtained (at best, the RMS error is 1mm). However the class provides methods to report the maximum and RMS errors in the results.</p>
<p>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.</p>
<h2><a class="anchor" id="geoidinst"></a>
Installing the datasets</h2>
<p>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:</p>
<ul>
<li>EGM84, 30' grid: <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm84-30.tar.bz2/download">tar.bz2</a>, <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm84-30.zip/download">zip</a>.</li>
<li>EGM84, 15' grid: <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm84-15.tar.bz2/download">tar.bz2</a>, <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm84-15.zip/download">zip</a>.</li>
<li>EGM96, 15' grid: <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm96-15.tar.bz2/download">tar.bz2</a>, <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm96-15.zip/download">zip</a>.</li>
<li>EGM96, 5' grid: <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm96-5.tar.bz2/download">tar.bz2</a>, <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm96-5.zip/download">zip</a>.</li>
<li>EGM2008, 5' grid: <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-5.tar.bz2/download">tar.bz2</a>, <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-5.zip/download">zip</a>.</li>
<li>EGM2008, 2.5' grid: <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-2_5.tar.bz2/download">tar.bz2</a>, <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-2_5.zip/download">zip</a>.</li>
<li>EGM2008, 1' grid: <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-1.tar.bz2/download">tar.bz2</a>, <a href="http://sf.net/projects/geographiclib/files/geoids-distrib/egm2008-1.zip/download">zip</a>.</li>
</ul>
<p>Download <em>either</em> the tar.bz2 file <em>or</em> zip file (they have the same contents). If possible use the tar.bz2 format, since bzip2 compresses these files about 2 times better than zip. The data in tar and zip files both unpack into a "geoids" directory and, for convenience, all the data should be unpacked in the same parent directory. For example </p>
<div class="fragment"><pre class="fragment">
mkdir -p /usr/local/share/GeographicLib
tar xfjC egm96-5.tar.bz2 /usr/local/share/GeographicLib
tar xfjC egm2008-2_5.tar.bz2 /usr/local/share/GeographicLib
etc.
</pre></div><p> or </p>
<div class="fragment"><pre class="fragment">
mkdir -p /usr/local/share/GeographicLib
cd /usr/local/share/GeographicLib
unzip ~-/egm96-5.zip
unzip ~-/egm2008-2_5.zip
etc.
</pre></div><p> Geoid uses a compile time default (specified by the GEOID_DEFAULT_PATH macro) to locate the datasets. This is</p>
<ul>
<li>/usr/local/share/GeographicLib/geoids, for non-Windows systems</li>
<li>C:/cygwin/usr/local/share/GeographicLib/geoids, for Windows systems</li>
</ul>
<p>This may be overridden at run-time by defining the GEOID_PATH environment variable. Finally, the path may be set using the optional second argument to the GeographicLib::Geoid::Geoid constructor.</p>
<p>The sizes of these datasets are </p>
<div class="fragment"><pre class="fragment">
name geoid grid sizes (MB)
tar.bz2 zip disk
egm84-30 EGM84 30' 0.5 0.5 0.5
egm84-15 EGM84 15' 1.5 1.9 2.0
egm96-15 EGM96 15' 1.5 1.9 2.0
egm96-5 EGM96 5' 9.8 16 18
egm2008-5 EGM2008 5' 10 16 18
egm2008-2_5 EGM2008 2.5' 34 62 72
egm2008-1 EGM2008 1' 155 372 445
</pre></div><p> The final "disk" column is the size of the uncompressed data and it also gives the memory requirements for caching the entire dataset using the <a class="el" href="classGeographicLib_1_1Geoid.html#a482c6482d5ab4c5d661210327848170e">GeographicLib::Geoid::CacheAll</a> method.</p>
<h2><a class="anchor" id="geoidformat"></a>
The format of the geoid data files</h2>
<p>The gridded data used by the <a class="el" href="classGeographicLib_1_1Geoid.html" title="Looking up the height of the geoid.">GeographicLib::Geoid</a> class is stored in 16-bit PGM files. Thus the data for egm96-5 might be stored in the file</p>
<ul>
<li>/usr/local/share/GeographicLib/geoids/egm96-5.pgm</li>
</ul>
<p>PGM simple graphic format with the following properties</p>
<ul>
<li>it is well documented <a href="http://netpbm.sf.net/doc/pgm.html">here</a>;</li>
<li>there are no separate "little-endian" and "big-endian" versions;</li>
<li>it uses 1 or 2 bytes per pixel;</li>
<li>pixels can be randomly accessed;</li>
<li>comments can be added to the file header;</li>
<li>it is sufficiently simple that it can be easily read without using the <a href="http://netpbm.sf.net/doc/libnetpbm.html">libnetpbm</a> library (and thus we avoid adding a software dependency to GeographicLib).</li>
</ul>
<p>The major drawback of this format is that since there are only 65535 possible pixel values, the height must be quantized to 3mm. However, the resulting quantization error (up to 1.5mm) is typically smaller than the linear interpolation errors. The comments in the header for egm96-5 are </p>
<div class="fragment"><pre class="fragment">
# 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
</pre></div><p> Of these lines, the Scale and Offset lines are required and define the conversion from pixel value to height (in meters) using <em>height</em> = <em>offset</em> + <em>scale</em> <em>pixel</em>. The Geoid constructor also reads the Description, DataTime, and error lines (if present) and stores the resulting data so that it can be returned by <a class="el" href="classGeographicLib_1_1Geoid.html#a37d76bcfe0ddf9b84042d701c312e941">GeographicLib::Geoid::Description</a>, <a class="el" href="classGeographicLib_1_1Geoid.html#a9cd7304b5df37001f4ad3c91cdbc56f4">GeographicLib::Geoid::DateTime</a>, <a class="el" href="classGeographicLib_1_1Geoid.html#aff538da14578a02c551b411a899e567a">GeographicLib::Geoid::MaxError</a>, and <a class="el" href="classGeographicLib_1_1Geoid.html#a207f98316d20a5b9d47846e559c5e2a2">GeographicLib::Geoid::RMSError</a> 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 <a href="http://www.gdal.org">GDAL</a>. Neither of these files is read by <a class="el" href="classGeographicLib_1_1Geoid.html" title="Looking up the height of the geoid.">GeographicLib::Geoid</a>.</p>
<p>You can use gdal_translate to convert the data files to a standard GeoTiff, e.g., with </p>
<div class="fragment"><pre class="fragment">
gdal_translate -ot Float32 -scale 0 65000 -108 87 egm96-5.pgm egm96-5.tif
</pre></div><p> 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 <a href="GeoidEval.1.html">GeoidEval</a> with the "-v" option.</p>
<p>Finally, 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. </p>
<div class="fragment"><pre class="fragment">
#! /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
</pre></div><h2><a class="anchor" id="geoidinterp"></a>
Interpolating the geoid data</h2>
<p><a class="el" href="classGeographicLib_1_1Geoid.html" title="Looking up the height of the geoid.">GeographicLib::Geoid</a> evaluates the geoid height is by 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.</p>
<p>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.</p>
<p>The cubic interpolation is a least-squares fit to the values on a 12-point stencil with weights as follows: </p>
<div class="fragment"><pre class="fragment">
. 1 1 .
1 2 2 1
1 2 2 1
. 1 1 .
</pre></div><p> 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.</p>
<p>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 Geographic::Geoid; instead it is implemented as <a href="http://en.wikipedia.org/wiki/Maxima_(software)">Maxima</a> 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 <a class="el" href="Geoid_8cpp.html" title="Implementation for GeographicLib::Geoid class.">Geoid.cpp</a>.</p>
<p>The interpolations methods are quick and give good accuracy. Here is a summary of the combined quantization and interpolation errors for the heights. </p>
<div class="fragment"><pre class="fragment">
bilinear error cubic error
name geoid grid max rms max rms
egm84-30 EGM84 30' 1.546m 70mm 0.274m 14mm
egm84-15 EGM84 15' 0.413m 18mm 0.020m 1mm
egm96-15 EGM96 15' 1.152m 40mm 0.169m 7mm
egm96-5 EGM96 5' 0.140m 5mm 0.003m 1mm
egm2008-5 EGM2008 5' 0.478m 12mm 0.294m 5mm
egm2008-2_5 EGM2008 2.5' 0.135m 3mm 0.031m 1mm
egm2008-1 EGM2008 1' 0.025m 1mm 0.003m 1mm
</pre></div><p> The errors are with respect the 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 or near the centers of the original grid. (The RMS difference between EGM96 and EGM2008 is about 0.5m. The RMS difference between EGM84 and EGM96 is about 1.5m.)</p>
<h2><a class="anchor" id="geoidcache"></a>
Caching the geoid data</h2>
<p>A simple way of calling Geoid is as follows </p>
<div class="fragment"><pre class="fragment"><span class="preprocessor"> #include "<a class="code" href="Geoid_8hpp.html" title="Header for GeographicLib::Geoid class.">GeographicLib/Geoid.hpp</a>"</span>
<span class="preprocessor"> #include <iostream></span>
...
<a class="code" href="classGeographicLib_1_1Geoid.html" title="Looking up the height of the geoid.">GeographicLib::Geoid</a> g(<span class="stringliteral">"egm96-5"</span>);
<span class="keywordtype">double</span> lat, lon;
<span class="keywordflow">while</span> (std::cin >> lat >> lon)
std::cout << g(lat, lon) << <span class="stringliteral">"\n"</span>;
...
</pre></div><p>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.</p>
<p>If a large quantity of geoid calculations are needed, the calculation can be speeded up by preloading the data for a rectangular block with <a class="el" href="classGeographicLib_1_1Geoid.html#a52b5dc2d976796046aaeb8765e4a9c0f">GeographicLib::Geoid::CacheArea</a> or the entire dataset with <a class="el" href="classGeographicLib_1_1Geoid.html#a482c6482d5ab4c5d661210327848170e">GeographicLib::Geoid::CacheAll</a>. 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.5GB of RAM and should only be used on systems with sufficient memory.</p>
<p>The use of caching does not affect the values returned. Because of the caching and the random file access, this class is <em>not</em> 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:</p>
<ul>
<li>they should all construct thread-local instantiations.</li>
<li><a class="el" href="classGeographicLib_1_1Geoid.html" title="Looking up the height of the geoid.">GeographicLib::Geoid</a> should be constructed with <em>threadsafe</em> = 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.</li>
</ul>
<h2><a class="anchor" id="testgeoid"></a>
Test data for geoids</h2>
<p>A test set for the geoid models is available at</p>
<ul>
<li><a href="http://sf.net/projects/geographiclib/files/testdata/GeoidHeights.dat.gz/download">GeoidHeight.dat.gz</a></li>
</ul>
<p>This is about 10 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.</p>
<p>Each line of the test set gives 6 space delimited numbers</p>
<ul>
<li>latitude (degrees, exact)</li>
<li>longitude (degrees, exact)</li>
<li>EGM84 height (meters, accurate to 1 mm)</li>
<li>EGM96 height (meters, accurate to 0.1 mm)</li>
<li>EGM2008 height (meters, accurate to 1 mm)</li>
</ul>
<p>The latitude and longitude are all multiples of 10<sup>-6</sup> deg and should be regarded as exact. The geoid heights are computed using the harmonic NGA synthesis programs. In the case of the EGM84 and EGM96, the programs were compiled (with gfortran) and run under Linux. In the case of EGM2008, the prescription results in NaNs (possibly because of uninitialized variables in the program?), so the data was generated with the precompiled binary that NGA supplies for Windows.</p>
<center> Back to <a class="el" href="geocentric.html">Geocentric coordinates</a>. Forward to <a class="el" href="utilities.html">Utility Programs</a>. Up to <a class="el" href="index.html#contents">Contents</a>. </center> </div>
<hr class="footer"/><address class="footer"><small>Generated on Tue Feb 22 2011 for GeographicLib by
<a href="http://www.doxygen.org/index.html">
<img class="footer" src="doxygen.png" alt="doxygen"/></a> 1.7.1 </small></address>
</body>
</html>
|