When approximate is better than exact

The subject of this post is a computational problem that is solved thousands of times each day (not just as part of my day job, but in almost every GPS receiver in the world).  It is a problem with an interesting history, and dozens of different solutions… but the arguably “best” solution seems to be relatively unknown.

The problem is that of converting Earth-centered Cartesian (x,y,z) coordinates to geodetic latitude, longitude, and height.  The transformation from geodetic to Cartesian coordinates is straightforward, but the reverse transformation is significantly more involved, and as a result has been addressed repeatedly in the literature.  Featherstone (see Reference (4) at the end of this post) provides a nearly comprehensive list of over 70 papers describing solutions or comparisons of solutions, spanning five decades of study.

I think this problem is interesting for three reasons.  First, it still seems to be not universally known that closed-form solutions even exist, usually obtained by finding the roots of an appropriately transformed quartic equation.  Even the first non-Wikipedia entry in a quick Google search claims that “there is no closed-form solution… if the altitude is not zero.”

Second, although all of these various closed-form solutions are “equal,” some are more equal than others.  That is, the accuracies of many of these “exact” solutions suffer when executed on a computer with limited precision.  Or worse, many papers only provide equations, with actual implementations requiring additional details not made explicit, such as handling singularities at the poles, or taking square roots of expressions that “should” be non-negative, but in double precision may yield negative values.  Borkowski’s formula in Reference (1), for example, breaks down for a spherical earth (i.e., zero eccentricity), when the transformation is in fact simplest.  And Bowring’s formula in Reference (2), which I see most often in others’ code, has a usually-neglected but easily-remedied danger of a singularity at the poles.

Which brings us to the third interesting feature of this problem: there is a solution provided by Olson in Reference (6), that, mathematically, is actually an approximation (essentially involving a truncated series expansion), but when implemented in code, and executed with limited double-precision, is not only more accurate than all of the “exact” closed-form solutions that I have evaluated in the references, but is also faster.

(It is worth emphasizing the distinction between this approximation, which although inexact is still essentially “closed-form,” and a whole class of other solutions that involve iterative approximation (e.g., Newton-Raphson), where execution time is more difficult to predict.)

Olson also actually provides working C code as a bonus; following is a port of his implementation to Python:

import math

# WGS-84 ellipsoid parameters
a = 6378137
f = 1 / 298.257223563

# Derived parameters
e2 = f * (2 - f)
a1 = a * e2
a2 = a1 * a1
a3 = a1 * e2 / 2
a4 = 2.5 * a2
a5 = a1 + a3
a6 = 1 - e2

def ecef_to_lla(ecef):
    """Convert ECEF (meters) to LLA (radians and meters).
    """
    # Olson, D. K., Converting Earth-Centered, Earth-Fixed Coordinates to
    # Geodetic Coordinates, IEEE Transactions on Aerospace and Electronic
    # Systems, 32 (1996) 473-476.
    w = math.sqrt(ecef[0] * ecef[0] + ecef[1] * ecef[1])
    z = ecef[2]
    zp = abs(z)
    w2 = w * w
    r2 = z * z + w2
    r  = math.sqrt(r2)
    s2 = z * z / r2
    c2 = w2 / r2
    u = a2 / r
    v = a3 - a4 / r
    if c2 > 0.3:
        s = (zp / r) * (1 + c2 * (a1 + u + s2 * v) / r)
        lat = math.asin(s)
        ss = s * s
        c = math.sqrt(1 - ss)
    else:
        c = (w / r) * (1 - s2 * (a5 - u - c2 * v) / r)
        lat = math.acos(c)
        ss = 1 - c * c
        s = math.sqrt(ss)
    g = 1 - e2 * ss
    rg = a / math.sqrt(g)
    rf = a6 * rg
    u = w - rg * c
    v = zp - rf * s
    f = c * u + s * v
    m = c * v - s * u
    p = m / (rf / g + f)
    lat = lat + p
    if z < 0:
        lat = -lat
    return (lat, math.atan2(ecef[1], ecef[0]), f + m * p / 2)

References:

1. K. M. Borkowski, “Accurate Algorithms to Transform Geocentric to Geodetic Coordinates,” Bulletin Geodesique, 63 (1989) 50-56 [HTML]
2. B. R. Bowring, “Transformation from Spatial to Geographical Coordinates,” Survey Review, 23:181 (1976) 323-327
3. R. Burtch, “A Comparison of Methods Used in Rectangular to Geodetic Coordinate Transformations,” Proceedings of the American Congress on Surveying and Mapping (ACSM) Annual Conference and Technology Exhibition, Orlando, FL, 2006
4. W. E. Featherstone and S. J. Claessens, “Closed-Form Transformation Between Geodetic and Ellipsoidal Coordinates,” Studia Geophysica et Geodaetica, 52 (2008) 1-18
5. M. Heikkinen, “Geschlossene Formeln zur Berechnung raumlicher geodatischer Koordinaten aus rechtwinkligen Koordinaten,” Zeitschrift Vermess. (in German), 107 (1982) 207-211
6. D. K. Olson, “Converting Earth-Centered, Earth-Fixed Coordinates to Geodetic Coordinates,” IEEE Transactions on Aerospace and Electronic Systems, 32 (1996) 473-476
7. I. Sofair, “Improved Method for Calculating Exact Geodetic Latitude and Altitude Revisited,” Journal of Guidance, Control, and Dynamics, 23:2 (2000) 369
8. J. Zhu, “Conversion of Earth-Centered Earth-Fixed Coordinates to Geodetic Coordinates,” IEEE Transactions on Aerospace and Electronic Systems, 30:3 (1994) 957-962
9. J. Zhu, “Exact Conversion of Earth-Centered, Earth-Fixed Coordinates to Geodetic Coordinates,” Journal of Guidance, Control, and Dynamics, 16:2 (1993) 389-391