Back to Transverse Mercator projection. Forward to Auxiliary latitudes. Up to Contents.

The implementation of 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 changes are described in Appendix B of

C. F. F. Karney, Geodesics on an ellipsoid of revolution, Feb. 2011; preprint arxiv:1102.1215v1.

Vermeille similarly updated his method in

H. Vermeille, An analytical method to transform geocentric into geodetic coordinates, J. Geodesy 85, 105–117 (2011).

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 (7 nanometers).

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 (1 nm = 1 nanometer).

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 analyzed for this case. Similarly ellipsoids which are very nearly spherical may yield inaccurate results due to underflow; in the other hand, the case of the sphere, f = 0, is treated specially and gives accurate results.

Other comparable methods are K. M. Borkowski, Transformation of geocentric to geodetic coordinates without approximations, Astrophys. Space Sci. 139, 1–4 (1987) (erratum) and T. Fukushima, Fast transform from geocentric to geodetic coordinates, J. Geodesy 73, 603–610 (1999). However the choice of independent variables in these methods leads to a loss of accuracy for points near the equatorial plane.

Back to Transverse Mercator projection. Forward to Auxiliary latitudes. Up to Contents.