Geocentric coordinates

Back to Geodesics on the Ellipsoid. Forward to Utility Programs. Up to Contents.

The implementation of GeographicLib::Geocentric::Reverse is adapted from

• H. Vermeille, Direct transformation from geocentric coordinates to geodetic coordinates, J. Geodesy 76, 451–454 (2002).

This provides a closed-form solution but can't directly be applied close to the center of the earth. Several changes have been made to remove this restriction and to improve the numerical accuracy. Now the method is accurate for all inputs (even if h is infinite).

The problems encountered near the center of the ellipsoid are:

• There's a potential division by zero in the definition of s. The equations are easily reformulated to avoid this problem.
• t³ may be negative. This is OK; we just take the real root.
• The solution for t may be complex. However this leads to 3 real roots for u/r. It's then just a matter of picking the one that computes the geodetic result which minimizes |h| and which avoids large round-off errors.
• Some of the equations result in a large loss of accuracy due to subtracting nearly equal quantities. E.g., k= sqrt(u + v + w²) - w is inaccurate if u + v is small; we can fix this by writing k = (u + v)/(sqrt(u + v + w²) + w).

The error is computed as follows. Write a version of Geocentric::WGS84.Forward which uses long doubles (including using long doubles for the WGS84 parameters). Generate random (long double) geodetic coordinates (lat0, lon0, h0) and use the "long double" WGS84.Forward to obtain the corresponding (long double) geocentric coordinates (x0, y0, z0). [We restrict h0 so that h0 >= - a (1 - e²) / sqrt(1 - e² sin²lat0), which ensures that (lat0, lon0, h0) is the principal geodetic inverse of (x0, y0, z0).] Because the forward calculation is numerically stable and because long doubles (on Linux systems using g++) provide 11 bits additional accuracy (about 3.3 decimal digits), we regard this set of test data as exact.

Apply the double version of WGS84.Reverse to (x0, y0, z0) to compute the approximate geodetic coordinates (lat1, lon1, h1). Convert (lat1 - lat0, lon1 - lon0) to a distance, ds, on the surface of the ellipsoid and define err = hypot(ds, h1 - h0). For |h0| < 5000 km, we have err < 7 nm.

This methodology is not very useful very far from the globe, because the absolute errors in the approximate geodetic height become large, or within 50 km of the center of the earth, because of errors in computing the approximate geodetic latitude. To illustrate the second issue, the maximum value of err for h0 < 0 is about 80 mm. The error is maximum close to the circle given by geocentric coordinates satisfying hypot(x, y) = a e² (= 42.7 km), z = 0. (This is the center of meridional curvature for lat = 0.) The geodetic latitude for these points is lat = 0. However, if we move 1 nm towards the center of the earth, the geodetic latitude becomes 0.04", a distance of 1.4 m from the equator. If, instead, we move 1 nm up, the geodetic latitude becomes 7.45", a distance of 229 m from the equator. In light of this, Reverse does quite well in this vicinity.

To obtain a practical measure of the error for the general case we define

• err_h = |h1 - h0| / max(1, h0 / a)
• for h0 > 0, err_out = ds
• for h0 < 0, apply the long double version of WGS84.Forward to (lat1, lon1, h1) to give (x1, y1, z1) and compute err_in = hypot(x1 - x0, y1 - y0, z1 - z0).

We then find err_h < 8 nm, err_out < 4 nm, and err_in < 7 nm.

The testing has been confined to the WGS84 ellipsoid. The method will work for all ellipsoids used in terrestrial geodesy. However, the central region, which leads to multiple real roots for the cubic equation in Reverse, pokes outside the ellipsoid (at the poles) for ellipsoids with e > 1/sqrt(2). Reverse has not been analysed for this case. Similarly ellipsoids which are very nearly spherical near yield inaccurate results due to underflow; in the other hand, the case of the sphere, f = 0, is treated specially and gives accurate results.

Another comparable method is T. Fukushima, Fast transform from geocentric to geodetic coordinates, J. Geodesy 73, 603–610 (2003). This is an iterative method and is somewhat faster than Geocentric.Reverse. However, because of the choice of independent variable in Newton's iteration, accuracy is lost for points near the equatorial plane. As a consequence, the maximum error err near the center of meridional curvature for lat = 0 is about 50 m (as opposed to 8 mm for WGS84.Reverse).

Back to Geodesics on the Ellipsoid. Forward to Utility Programs. Up to Contents.