GeographicLib: Gravity models (original) (raw)

Back to Geoid height. Forward to Normal gravity. Up to Contents.

GeographicLib can compute the earth's gravitational field with an earth gravity model using the GravityModel and GravityCircle classes and with the Gravity utility. These models expand the gravitational potential of the earth as sum of spherical harmonics. The models also specify a reference ellipsoid, relative to which geoid heights and gravity disturbances are measured. Underlying all these models is normal gravity which is the exact solution for an idealized rotating ellipsoidal body. This is implemented with the NormalGravity class and some notes on are provided in section Normal gravity

Go to

The supported models are

See

for more information.

Acknowledgment: I would like to thank Mathieu Peyréga for sharing EGM_Geoid_CalculatorClass from his Geo library with me. His implementation was the first I could easily understand and he and I together worked through some of the issues with overflow and underflow the occur while performing high-degree spherical harmonic sums.

Installing the gravity models

These gravity models are available for download:

Available gravity models

name max degree size(kB) Download Links (size, kB)
tar file Windows installer zip file
egm84 180 27 link (26) link (55) link (26)
egm96 360 2100 link (2100) link (2300) link (2100)
egm2008 2190 76000 link (75000) link (72000) link (73000)
wgs84 20 1 link (1) link (30) link (1)

The "size" column is the size of the uncompressed data.

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

geographiclib-get-gravity all # to install egm84, egm96, egm2008, wgs84 geographiclib-get-gravity -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 "gravity" 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 "gravity" 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). To unpack these, run, for example

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

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

However you install the gravity models, all the datasets should be installed in the same directory. GravityModel and Gravity uses a compile time default to locate the datasets. This is

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

The format of the gravity model files

The constructor for GravityModel reads a file called NAME.egm which specifies various properties for the gravity model. It then opens a binary file NAME.egm.cof to obtain the coefficients of the spherical harmonic sum.

The first line of the .egm file must consist of "EGMF-v" where EGMF stands for "Earth Gravity Model Format" and v is the version number of the format (currently "1").

The rest of the file is read a line at a time. A # character and everything after it are discarded. If the result is just white space it is discarded. The remaining lines are of the form "KEY WHITESPACE VALUE". In general, the KEY and the VALUE are case-sensitive.

GravityModel only pays attention to the following keywords

Other keywords are ignored.

The coefficient file NAME.egm.cof is a binary file in little endian order. The first 8 bytes of this file must match the ID given in NAME.egm. This is followed by 2 sets of spherical harmonic coefficients. The first of these gives the gravity potential and the second gives the zeta-to-N corrections to the geoid height. The format for each set of coefficients is:

Although the coefficient file is in little endian order, GeographicLib can read it on big endian machines. It can only be read on machines which store doubles in IEEE format.

As an illustration, here is egm2008.egm:

EGMF-1

An Earth Gravity Model (Format 1) file. For documentation on the

format of this file see

https://geographiclib.sourceforge.io/C++/doc/gravity.html#gravityformat

Name egm2008 Publisher National Geospatial Intelligence Agency Description Earth Gravity Model 2008 URL https://earth-info.nga.mil/index.php?dir=wgs84&action=wgs84#tab_egm2008 ReleaseDate 2008-06-01 ConversionDate 2011-11-19 DataVersion 1 ModelRadius 6378136.3 ModelMass 3986004.415e8 AngularVelocity 7292115e-11 ReferenceRadius 6378137 ReferenceMass 3986004.418e8 Flattening 1/298.257223563 HeightOffset -0.41

Gravitational and correction coefficients taken from

EGM2008_to2190_TideFree and Zeta-to-N_to2160_egm2008 from

the egm2008 distribution.

ID EGM2008A

The code to produce the coefficient files for the wgs84 and grs80 models is

#include

#include

#include

using namespace std;

try {

Utility::set_digits();

const char* filenames[] = {"wgs84.egm.cof", "grs80.egm.cof"};

const char* ids[] = {"WGS1984A", "GRS1980A"};

for (int grs80 = 0; grs80 < 2; ++grs80) {

ofstream file(filenames[grs80], ios::binary);

Utility::writearray<char, char, false>(file, ids[grs80], 8);

const int N = 20, M = 0,

cnum = (M + 1) * (2 * N - M + 2) / 2;

vector num(2);

num[0] = N; num[1] = M;

Utility::writearray<int, int, false>(file, num);

vectorMath::real c(cnum, 0);

const NormalGravity& earth(grs80 ? NormalGravity::GRS80() :

for (int n = 2; n <= N; n += 2)

Utility::writearray<double, Math::real, false>(file, c);

num[0] = num[1] = -1;

Utility::writearray<int, int, false>(file, num);

}

}

catch (const exception& e) {

cerr << "Caught exception: " << e.what() << "\n";

return 1;

}

catch (...) {

cerr << "Caught unknown exception\n";

return 1;

}

}

int main(int argc, const char *const argv[])

Header for GeographicLib::NormalGravity class.

Header for GeographicLib::Utility class.

The normal gravity of the earth.

Math::real DynamicalFormFactor(int n=2) const

Namespace for GeographicLib.

Comments on the NGA harmonic synthesis code

GravityModel attempts to reproduce the results of NGA's harmonic synthesis code for EGM2008, hsynth_WGS84.f. Listed here are issues that I encountered using the NGA code:

  1. A compiler which allocates local variables on the stack produces an executable with just returns NaNs. The problem here is a missing SAVE statement in subroutine latf.
  2. Because the model and references masses for egm2008 differ (by about 1 part in 109), there should be a 1/r contribution to the disturbing potential T. However, this term is set to zero in hsynth_WGS84 (effectively altering the normal potential). This shifts the surface W = _U_0 outward by about 5 mm. Note too that the reference ellipsoid is no longer a level surface of the effective normal potential.
  3. Subroutine radgrav computes the ellipsoidal coordinate β incorrectly. This leads to small errors in the deflection of the vertical, ξ and η, when the height above the ellipsoid, h, is non-zero (about 10−7 arcsec at h = 400 km).
  4. There are several expressions which will return inaccurate results due to cancellation. For example, subroutine grs computes the flattening using f = 1 − sqrt(1 − _e_2). Much better is to use f = _e_2/(1 + sqrt(1 − _e_2)). The expressions for q and q' in grs and radgrav suffer from similar problems. The resulting errors are tiny (about 50 pm in the geoid height); however, given that's there's essentially no cost to using more accurate expressions, it's preferable to do so.
  5. hsynth_WGS84 returns an "undefined" value for xi and eta at the poles. Better would be to return the value obtained by taking the limit lat → ± 90°.

Issues 1–4 have been reported to the authors of hsynth_WGS84. Issue 1 is peculiar to Fortran and is not encountered in C++ code and GravityModel corrects issues 3–5. On issue 2, GravityModel neglects the 1/r term in T in GravityModel::GeoidHeight and GravityModel::SphericalAnomaly in order to produce results which match NGA's for these quantities. On the other hand, GravityModel::Disturbance and GravityModel::T do include this term.

Details of the geoid height and anomaly calculations

Ideally, the geoid represents a surface of constant gravitational potential which approximates mean sea level. In reality some approximations are taken in determining this surface. The steps taking by GravityModel in computing the geoid height are described here (in the context of EGM2008). This mimics NGA's code hsynth_WGS84 closely because most users of EGM2008 use the gridded data generated by this code (e.g., Geoid) and it is desirable to use a consistent definition of the geoid height.

A useful discussion of the problems with defining a geoid is given by Dru A. Smith in There is no such thing as "The" EGM96 geoid: Subtle points on the use of a global geopotential model, IGeS Bulletin No. 8, International Geoid Service, Milan, Italy, pp. 17–28 (1998).

GravityModel::GeoidHeight reproduces the results of the several NGA codes for harmonic synthesis with the following maximum discrepancies:

The formula for the gravity anomaly vector involves computing gravity and normal gravity at two different points (with the displacement between the points unknown ab initio). Since the gravity anomaly is already a small quantity it is sometimes acceptable to employ approximations that change the quantities by O(f). The NGA code uses the spherical approximation described by Heiskanen and Moritz, Sec. 2-14 and GravityModel::SphericalAnomaly uses the same approximation for compatibility. In this approximation, the gravity disturbance delta = grad T is calculated. Here, T once again excludes the 1/r term (this is issue 2 in Comments on the NGA harmonic synthesis code and is consistent with the computation of the geoid height). Note that delta compares the gravity and the normal gravity at the same point; the gravity anomaly vector is then found by estimating the gradient of the normal gravity in the limit that the earth is spherically symmetric. delta is expressed in spherical coordinates as deltax, deltay, deltaz where, for example, deltaz is the radial component of delta (not the component perpendicular to the ellipsoid) and deltay is similarly slightly different from the usual northerly component. The components of the anomaly are then given by

NormalGravity computes the normal gravity accurately and avoids issue 3 of Comments on the NGA harmonic synthesis code. Thus while GravityModel::SphericalAnomaly reproduces the results for xi and eta at h = 0, there is a slight discrepancy if h is non-zero.

The effect of the mass of the atmosphere

All of the supported models use WGS84 for the reference ellipsoid. This has (see TR8350.2, table 3.1)

The value of GM includes the mass of the atmosphere and so strictly only applies above the earth's atmosphere. Near the surface of the earth, the value of g will be less (in absolute value) than the value predicted by these models by about δ_g_ = (4π_G_/g) A = 8.552 × 10−11 A m2/kg, where G is the gravitational constant, g is the earth's gravity, and A is the pressure of the atmosphere. At sea level we have A = 101.3 kPa, and δ_g_ = 8.7 × 10−6 m s−2, approximately. (In other words the effect is about 1 part in a million; by way of comparison, buoyancy effects are about 100 times larger.)

Geoid heights on a multi-processor system

The egm2008 model includes many terms (over 2 million spherical harmonics). For that reason computations using this model may be slow; for example it takes about 78 ms to compute the geoid height at a single point. There are two ways to speed up this computation:

Both of these techniques are illustrated by the following code, which computes a table of geoid heights on a regular grid and writes on the result in a .gtx file. On an 8-processor Intel 3.0 GHz machine using OpenMP (-DHAVE_OPENMP=1), it takes about 10 minutes of elapsed time to compute the geoid height for EGM2008 on a 1' grid. (Without these optimizations, the computation would have taken about 50 days!)

#include

#include

#include

#include

#include

#if defined(_OPENMP)

#define HAVE_OPENMP 1

#else

#define HAVE_OPENMP 0

#endif

#if HAVE_OPENMP

# include <omp.h>

#endif

using namespace std;

int main(int argc, const char* const argv[]) {

if (argc != 4) {

cerr << "Usage: " << argv[0]

<< " gravity-model intervals-per-degree output.gtx\n";

return 1;

}

try {

int ndigits = Utility::set_digits();

string model(argv[1]);

int ndeg = Utility::val(string(argv[2]));

string filename(argv[3]);

int

nlat = Math::hd * ndeg + 1,

nlon = Math::td * ndeg;

delta = 1 / Math::real(ndeg),

latorg = -Math::qd,

lonorg = -Math::hd;

ofstream file(filename.c_str(), ios::binary);

{

Math::real transform[] = {latorg, lonorg,

unsigned sizes[] = {unsigned(nlat), unsigned(nlon)};

Utility::writearray<double, Math::real, true>(file, transform, 4);

Utility::writearray<unsigned, unsigned, true>(file, sizes, 2);

}

const int nbatch = 64;

vector< vector > N(nbatch, vector(nlon));

for (int ilat0 = 0; ilat0 < nlat; ilat0 += nbatch) {

int nlat0 = min(nlat, ilat0 + nbatch);

#if HAVE_OPENMP

# pragma omp parallel for

#endif

for (int ilat = ilat0; ilat < nlat0; ++ilat) {

Utility::set_digits(ndigits);

Math::real lat = (latorg * ndeg + ilat) / ndeg, h = 0;

GravityCircle c(g.Circle(lat, h, GravityModel::GEOID_HEIGHT));

for (int ilon = 0; ilon < nlon; ++ilon) {

Math::real lon = (lonorg * ndeg + ilon) / ndeg;

N[ilat - ilat0][ilon] = float(c.GeoidHeight(lon));

}

}

for (int ilat = ilat0; ilat < nlat0; ++ilat)

Utility::writearray<float, float, true>(file, N[ilat - ilat0]);

}

}

catch (const exception& e) {

cerr << "Caught exception: " << e.what() << "\n";

return 1;

}

catch (...) {

cerr << "Caught unknown exception\n";

return 1;

}

}

Header for GeographicLib::GravityCircle class.

Header for GeographicLib::GravityModel class.

Gravity on a circle of latitude.

Model of the earth's gravity field.

cmake will add in support for OpenMP for examples/GeoidToGTX.cpp, if it is available.

Back to Geoid height. Forward to Normal gravity. Up to Contents.