Geodetic coordinates

Theo Martel downloaded the final coordinates of pillar CR1 from the NRCan Precise Point Positioning service and found them in two forms: latitude 49° 41′ 48.2301″ N, longitude 112° 49′ 37.5164″ W and ellipsoidal height 893.412 m, and Cartesian coordinates X, Y and Z, each in the millions of metres, two of them negative. Both describe the same point in NAD83(CSRS). Field software then expressed the vector from CR1 to CR2 in a third form, as east, north and up components. Each form has its uses, and converting between them, exactly and without loss, is one of the most common computations in modern surveying. This lesson defines geodetic coordinates and the conversions to and from Earth-centred, Earth-fixed (ECEF) Cartesian coordinates and local east–north–up (ENU) coordinates.

Latitude, longitude and ellipsoidal height

Geodetic coordinates locate a point relative to a specified ellipsoid. The geodetic latitude \(\varphi\) is the angle between the ellipsoid normal through the point and the equatorial plane, positive north. The geodetic longitude \(\lambda\) is the angle in the equatorial plane between the reference meridian and the meridian plane of the point, positive east; points in North America have negative longitudes when signed numbers are used, or carry a W. The ellipsoidal height \(h\) is the distance from the ellipsoid to the point along the normal, positive outward.

Two details matter. First, the normal is not a line to the centre of the Earth; lesson 2 showed that it meets the rotation axis below the centre in the northern hemisphere. Second, the reference meridian of modern systems is not exactly the old Greenwich meridian through the Airy transit circle. The IERS Reference Meridian, used by ITRF and therefore by GNSS, lies about 100 m east of it, a consequence of how the satellite-based frame was aligned with the Earth's centre of mass. Visitors to Greenwich who stand on the brass line with a phone find their longitude is not zero for this reason.

Coordinates are always relative to a datum, meaning an ellipsoid placed and oriented in a particular way and realized by coordinates of actual stations. The same pillar has NAD83 coordinates that differ from its ITRF2020 coordinates by more than a metre, and both change with time as the continent moves. A latitude without a datum and, for precise work, an epoch, is incomplete.

Other latitudes

Geodetic latitude is the one surveyors use, but two others appear in formulas and software. The geocentric latitude \(\psi\) is the angle at the centre of the ellipsoid between the equatorial plane and the line to the point. The reduced (or parametric) latitude \(\beta\) is defined through the auxiliary circle of radius \(a\): project the point parallel to the minor axis onto that circle and measure the angle at the centre. For points on the ellipsoid,

so that \(\psi < \beta < \varphi\) in the northern hemisphere. The difference between geodetic and geocentric latitude is largest near 45°, at about 11.5′, which is roughly 21 km on the ground. Using the wrong latitude in a formula is therefore not a small error. Reduced latitude is used in the geodesic computations of lesson 4, because on the auxiliary sphere a geodesic behaves like a great circle. A fourth kind, astronomic latitude, is defined by the actual plumb line rather than an ellipsoid normal; lessons 9 and 10 treat it.

Earth-centred, Earth-fixed coordinates

An ECEF system is a right-handed Cartesian system with its origin at the Earth's centre of mass, its Z axis along the conventional rotation axis (the IERS Reference Pole), its X axis in the equatorial plane toward the reference meridian, and its Y axis completing the system at longitude 90° E. It rotates with the Earth, so a point fixed to the ground has fixed ECEF coordinates apart from crustal motion. GNSS receivers compute positions in this form, because satellite orbits are naturally expressed in it, and then convert to geodetic coordinates for display.

The conversion from geodetic to ECEF coordinates is direct:

with \(N\) the prime vertical radius of lesson 2. The factor \((1 - e^2)\) in Z reflects the fact that the normal meets the axis at a distance \(Ne^2\sin\varphi\) below the centre.

From ECEF back to geodetic coordinates

The reverse conversion is harder, because latitude and height are tangled together. Longitude is simple, \(\lambda = \operatorname{atan2}(Y, X)\). With \(p = \sqrt{X^2 + Y^2}\), the latitude satisfies

which contains \(\varphi\) on both sides through \(N\) and \(h\). The classical solution iterates: start with \(h = 0\), compute \(\varphi\), then \(N\) and \(h\), and repeat. Convergence is fast, typically to a fraction of a millimetre in three or four iterations for points near the Earth's surface. Closed-form and one-step methods also exist. Bowring's 1976 formula, which uses the parametric latitude for its first estimate, is accurate to well below a millimetre for heights up to a few thousand kilometres in a single step, and is used in many receivers. Near the poles the height formula with \(\cos\varphi\) becomes unstable, and an alternative form \(h = Z/\sin\varphi - N(1 - e^2)\) is used there.

Local east, north and up

ECEF coordinates are awkward for local work: a 500 m baseline shows up as changes in all three coordinates, none of which is horizontal or vertical. Rotating the ECEF difference vector into a local system at the start point gives east, north and up components:

The up axis is along the ellipsoid normal, not the plumb line, and the north axis lies in the plane tangent to the ellipsoid. From these components, the horizontal distance is \(\sqrt{e^2 + n^2}\), the geodetic azimuth is \(\operatorname{atan2}(e, n)\), and the ellipsoidal zenith angle follows from \(u\). Note that the up component of a long baseline includes the curvature of the Earth: the tangent plane rises above the ellipsoid with distance, so \(u\) is smaller than the difference in ellipsoidal height. Over 561 m the difference is about 25 mm.

# Geodetic <-> ECEF (Earth-fixed geocentric) on GRS80, with a local ENU vector
import math

a, f = 6378137.0, 1 / 298.257222101
e2 = f * (2 - f); b = a * (1 - f); ep2 = e2 / (1 - e2)

def dms(d, m, s): return math.copysign(abs(d) + m / 60 + s / 3600, d)

def to_ecef(lat, lon, h):
    p, l = math.radians(lat), math.radians(lon)
    N = a / math.sqrt(1 - e2 * math.sin(p) ** 2)
    return ((N + h) * math.cos(p) * math.cos(l),
            (N + h) * math.cos(p) * math.sin(l),
            (N * (1 - e2) + h) * math.sin(p))

def from_ecef(X, Y, Z):
    lon = math.atan2(Y, X)
    p = math.hypot(X, Y)
    lat = math.atan2(Z, p * (1 - e2))              # first guess
    for i in range(10):                            # simple fixed-point iteration
        N = a / math.sqrt(1 - e2 * math.sin(lat) ** 2)
        h = p / math.cos(lat) - N
        new = math.atan2(Z, p * (1 - e2 * N / (N + h)))
        if abs(new - lat) < 1e-13:
            break
        lat = new
    N = a / math.sqrt(1 - e2 * math.sin(new) ** 2)   # height from the final latitude
    return math.degrees(new), math.degrees(lon), p / math.cos(new) - N, i + 1

def bowring(X, Y, Z):                              # one-step closed form (Bowring 1976)
    p = math.hypot(X, Y)
    t = math.atan2(Z * a, p * b)
    lat = math.atan2(Z + ep2 * b * math.sin(t) ** 3, p - e2 * a * math.cos(t) ** 3)
    N = a / math.sqrt(1 - e2 * math.sin(lat) ** 2)
    return math.degrees(lat), p / math.cos(lat) - N

def enu(dX, dY, dZ, lat, lon):                     # rotate an ECEF difference to east, north, up
    p, l = math.radians(lat), math.radians(lon)
    E = -math.sin(l) * dX + math.cos(l) * dY
    N = -math.sin(p) * math.cos(l) * dX - math.sin(p) * math.sin(l) * dY + math.cos(p) * dZ
    U = math.cos(p) * math.cos(l) * dX + math.cos(p) * math.sin(l) * dY + math.sin(p) * dZ
    return E, N, U

CR1 = (dms(49, 41, 48.2301), -dms(112, 49, 37.5164), 893.412)
CR2 = (dms(49, 41, 49.8452), -dms(112, 49, 9.6235), 893.958)

X, Y, Z = to_ecef(*CR1)
print("CR1 ECEF: X =", format(X, ".4f"), " Y =", format(Y, ".4f"), " Z =", format(Z, ".4f"))
lat, lon, h, n = from_ecef(X, Y, Z)
print("back to geodetic after", n, "iterations: lat error",
      format((lat - CR1[0]) * 3600, ".1e"), "arcsec, h error", format(h - CR1[2], ".1e"), "m")
blat, bh = bowring(X, Y, Z)
print("Bowring one step: lat error", format((blat - CR1[0]) * 3600, ".1e"), "arcsec, h error", format(bh - CR1[2], ".1e"), "m")

psi = math.degrees(math.atan2(Z, math.hypot(X, Y)))           # geocentric latitude
beta = math.degrees(math.atan((1 - f) * math.tan(math.radians(CR1[0]))))  # reduced latitude
print("geodetic", round(CR1[0], 6), " geocentric", round(psi, 6), " reduced", round(beta, 6), "deg")
print("geodetic minus geocentric =", round((CR1[0] - psi) * 60, 3), "arcmin")

X2, Y2, Z2 = to_ecef(*CR2)
E, Nn, U = enu(X2 - X, Y2 - Y, Z2 - Z, CR1[0], CR1[1])
print("CR1 -> CR2 in ENU: E =", round(E, 4), " N =", round(Nn, 4), " U =", round(U, 4), "m")
print("horizontal distance", round(math.hypot(E, Nn), 4), "m, azimuth", round(math.degrees(math.atan2(E, Nn)) % 360, 4), "deg")

CR1 ECEF: X = -1603877.9999  Y = -3810431.1941  Z = 4841719.9163
back to geodetic after 4 iterations: lat error 2.6e-11 arcsec, h error 4.8e-10 m
Bowring one step: lat error 2.3e-10 arcsec, h error 7.9e-09 m
geodetic 49.696731  geocentric 49.50681  reduced 49.601783 deg
geodetic minus geocentric = 11.395 arcmin
CR1 -> CR2 in ENU: E = 559.0602  N = 49.9349  U = 0.5214 m
horizontal distance 561.2859 m, azimuth 84.8959 deg

The round trip recovers the latitude to better than \(10^{-10}\) arcseconds and the height to a fraction of a nanometre, so the iteration loses nothing; Bowring's single step is equally good. The geocentric latitude of CR1 is 11.4′ smaller than its geodetic latitude. The vector to CR2 has an east component of 559.060 m, a north component of 49.935 m and an up component of 0.521 m. The ellipsoidal heights of the two pillars differ by 0.546 m; the up component is 25 mm smaller because the tangent plane at CR1 rises above the ellipsoid by \(s^2/2R\) over 561 m. The geodetic azimuth from CR1 to CR2 is 84.8959°.

Formats and precision

Latitude and longitude appear as degrees, minutes and seconds, as decimal degrees, and in software as radians. The table shows what each decimal place is worth on the ground. Survey coordinates in seconds are normally quoted to four or five decimals, and decimal degrees to eight or nine. GIS files that store coordinates as single-precision floating-point numbers keep only about seven significant digits, which for a longitude near −112.8° means about a metre: a hidden loss that has spoiled more than one survey deliverable.

Signs need care too. West longitudes are negative in most software but positive in some older formats, and a dropped minus sign puts a point in China. Mixing degrees and radians in trigonometric functions produces nonsense rather than an error message. Both mistakes are easy to catch with a habit: after any conversion, plot the point or compare it with a known nearby station.