The ellipsoid

The geoid, the level surface that best fits mean sea level, is the natural reference for heights, but it is far too irregular to compute on. It rises and falls by up to about 100 m relative to any smooth figure, following the hidden distribution of mass inside the Earth. For horizontal positions and for most computations geodesists therefore use a simpler surface: an ellipsoid of revolution, formed by rotating an ellipse about its shorter axis. It is smooth, it can be described by two numbers, and its geometry can be computed to any precision. Every latitude and longitude a GNSS receiver reports, and every ellipsoidal height, is referred to such an ellipsoid. This lesson defines its parameters and the radii of curvature that govern distances on it.

Defining the ellipsoid

An ellipsoid of revolution is fixed by two parameters. The usual pair is the semi-major axis \(a\), the equatorial radius, and the flattening \(f\). The semi-minor axis \(b\), the polar radius, and the eccentricities follow:

Here \(e^2\) is the first eccentricity squared and \(e'^2\) the second. Flattening is a small number, about 1/298, so ellipsoids are usually quoted by their inverse flattening. The Earth's polar radius is about 21.4 km shorter than its equatorial radius. Drawn to scale on a page, the difference would be less than the thickness of the line, which is why every figure in this course exaggerates it.

Modern reference systems define the ellipsoid as part of a complete model of the Earth. The Geodetic Reference System 1980 (GRS80) uses four defining constants: the semi-major axis \(a = 6\,378\,137\) m, the geocentric gravitational constant \(GM = 3.986\,005 \times 10^{14}\ \text{m}^3/\text{s}^2\), the dynamic form factor \(J_2 = 0.001\,082\,63\), and the angular velocity \(\omega = 7.292\,115 \times 10^{-5}\) rad/s. The flattening is derived from these, and comes out as \(1/f = 298.257\,222\,101\). WGS 84, the system of GPS, took its \(a\) from GRS80 but fixed \(1/f = 298.257\,223\,563\) directly. The two ellipsoids differ in their semi-minor axes by about 0.1 mm, which is negligible for every practical purpose, though the systems that use them are realized differently and their frames are not identical.

The older ellipsoids were chosen to fit the geoid well over one continent, and they were placed so that their centres were offset from the Earth's centre of mass, sometimes by hundreds of metres. Clarke 1866 fits North America well but is 69 m larger at the equator than GRS80. Mixing parameters from different ellipsoids is a classic source of error: a coordinate computed on Clarke 1866 and converted with GRS80 formulas can be wrong by tens of metres.

Radii of curvature

A sphere has the same curvature in every direction; an ellipsoid does not. At any point, the curvature of a section through the ellipsoid normal depends on its azimuth. The two principal directions are north–south, along the meridian, and east–west, along the prime vertical, the section perpendicular to the meridian. Their radii of curvature are

where \(\varphi\) is the geodetic latitude. The meridian radius \(M\) is smallest at the equator, where it equals \(a(1-e^2)\), about 6335 km, and largest at the poles, where it equals \(a^2/b\), about 6400 km. The prime vertical radius \(N\) equals \(a\) at the equator and also \(a^2/b\) at the poles. Everywhere except the poles, \(N > M\). Note that \(N\) has a geometric meaning: it is the length of the normal from the surface to the rotation axis, as the figure shows. This is why it appears in the formulas for Cartesian coordinates in lesson 3.

For a section at any other azimuth \(\alpha\), Euler's theorem gives the radius

and the Gaussian mean radius, \(\sqrt{MN}\), is the radius of the sphere that best fits the ellipsoid near the point, in the sense that it has the same total curvature. It is the radius used for local spherical approximations, such as in map projection scale factors and in reducing distances to the ellipsoid. Finally, the radius of a parallel of latitude, the circle traced by the point as the Earth rotates, is \(N\cos\varphi\). A difference in longitude \(\Delta\lambda\) along a parallel spans an arc of \(N\cos\varphi\,\Delta\lambda\).

Keeping a few of these numbers in your head is useful in the field. One arcsecond of latitude is always close to 31 m; one arcsecond of longitude is 31 m times the cosine of the latitude. One second of time in longitude, 15 arcseconds, is about 230 to 400 m between latitudes 60° and 30°, which matters when longitude is found from time in lessons 9 and 11.

# Ellipsoid parameters and radii of curvature
import math

ELLIPSOIDS = [
    ("GRS80", 6378137.0, 298.257222101),
    ("WGS 84", 6378137.0, 298.257223563),
    ("Clarke 1866", 6378206.4, 294.9786982),
]

def derived(a, inv_f):
    f = 1 / inv_f
    b = a * (1 - f)
    e2 = f * (2 - f)                 # first eccentricity squared
    ep2 = e2 / (1 - e2)              # second eccentricity squared
    return f, b, e2, ep2

for name, a, inv_f in ELLIPSOIDS:
    f, b, e2, ep2 = derived(a, inv_f)
    print(name.ljust(12), "b =", format(b, ".4f"), "m  e2 =", format(e2, ".12f"), " e'2 =", format(ep2, ".12f"))

a, inv_f = ELLIPSOIDS[0][1:]
f, b, e2, ep2 = derived(a, inv_f)
print("GRS80 minus WGS 84 semi-minor axis:", round(b - derived(6378137.0, 298.257223563)[1], 4), "m")

def radii(lat_deg):
    s = math.sin(math.radians(lat_deg))
    W = math.sqrt(1 - e2 * s * s)
    M = a * (1 - e2) / W ** 3       # meridian (north-south) radius
    N = a / W                       # prime vertical (east-west) radius
    return M, N

def radius_in_azimuth(lat_deg, az_deg):
    M, N = radii(lat_deg)
    t = math.radians(az_deg)
    return 1 / (math.cos(t) ** 2 / M + math.sin(t) ** 2 / N)   # Euler's formula

LAT = 49.6967                       # site latitude in degrees
M, N = radii(LAT)
print()
print("At latitude", LAT, "deg:")
print("  M (meridian)        =", format(M, ",.3f"), "m")
print("  N (prime vertical)  =", format(N, ",.3f"), "m")
print("  Gaussian mean sqrt(MN) =", format(math.sqrt(M * N), ",.3f"), "m")
for az in (0, 45, 90, 135):
    print("  radius in azimuth", str(az).rjust(3), "deg =", format(radius_in_azimuth(LAT, az), ",.3f"), "m")
print("  radius of the parallel N cos(lat) =", format(N * math.cos(math.radians(LAT)), ",.3f"), "m")

R1 = (2 * a + b) / 3                # mean radius of the three semi-axes
print()
print("mean radius R1 = (2a + b)/3 =", format(R1, ",.4f"), "m")
print("one second of latitude here  =", round(M * math.radians(1 / 3600), 3), "m")
print("one second of longitude here =", round(N * math.cos(math.radians(LAT)) * math.radians(1 / 3600), 3), "m")

GRS80        b = 6356752.3141 m  e2 = 0.006694380023  e'2 = 0.006739496775
WGS 84       b = 6356752.3142 m  e2 = 0.006694379990  e'2 = 0.006739496742
Clarke 1866  b = 6356583.8000 m  e2 = 0.006768657998  e'2 = 0.006814784946
GRS80 minus WGS 84 semi-minor axis: -0.0001 m

At latitude 49.6967 deg:
  M (meridian)        = 6,372,620.705 m
  N (prime vertical)  = 6,390,589.991 m
  Gaussian mean sqrt(MN) = 6,381,599.023 m
  radius in azimuth   0 deg = 6,372,620.705 m
  radius in azimuth  45 deg = 6,381,592.698 m
  radius in azimuth  90 deg = 6,390,589.991 m
  radius in azimuth 135 deg = 6,381,592.698 m
  radius of the parallel N cos(lat) = 4,133,649.001 m

mean radius R1 = (2a + b)/3 = 6,371,008.7714 m
one second of latitude here  = 30.895 m
one second of longitude here = 20.04 m

The program reproduces the worked example and adds the radius in four azimuths. At 45° the radius lies between \(M\) and \(N\), close to the Gaussian mean. The difference between the GRS80 and WGS 84 semi-minor axes, 0.1 mm, confirms that the two ellipsoids can be treated as identical. The mean radius \(R_1 = (2a + b)/3 = 6\,371\,008.8\) m is the value used whenever a single Earth radius is needed, as in the curvature correction of lesson 1.

Meridian arc length and other radii

The distance along the meridian from the equator to latitude \(\varphi\) is the integral of \(M\,d\varphi\). It has no closed form, but a rapidly converging series in powers of the eccentricity, or of the third flattening \(n = (a-b)/(a+b)\), gives millimetre accuracy anywhere. The meridian arc length is the backbone of the transverse Mercator projection, because the northing of a point on the central meridian equals its arc length times the scale factor. The full meridian quadrant on GRS80 is 10,001,965.7 m, which is why the original definition of the metre came so close to the modern one.

Several mean radii are in use for different purposes. Besides \(R_1\), there is the authalic radius, of the sphere with the same surface area as the ellipsoid (about 6,371,007 m), and the volumetric radius, of the sphere with the same volume (about 6,371,001 m). For spherical computations over a region, the Gaussian mean radius at the centre of the region is the best choice.

When is a sphere good enough?

Spherical formulas are simpler, and over short distances they work. The question is how short. Over a line of length \(s\), using a sphere of radius \(R\) in place of the correct radius \(R_\alpha\) introduces relative errors in distances computed from angular coordinates of roughly \((R - R_\alpha)/R\). Using the mean radius 6371 km at the site, where the radius in the pipeline's azimuth is about 6383 km, gives a relative error of about 0.19 per cent, or 70 m in 38 km, as lesson 4 confirms. Using the local Gaussian mean radius instead cuts the error to a few parts per million for lines of a few kilometres. For anything beyond that, or for precise work, ellipsoidal formulas are required.

The same reasoning explains why surveyors can ignore the ellipsoid's shape in plane surveying over a building site but not over a pipeline route. The choice of ellipsoid also matters when old and new data meet. Coordinates from NAD27, on Clarke 1866, differ from NAD83 coordinates on GRS80 by tens of metres, partly because of the different ellipsoids and partly because the datums were positioned differently. Software handles the conversion, but only if the user identifies the datum correctly.