The following paper has been published in:
Bull. Géod., 63 (1989), pp. 50–56.

K.M. BORKOWSKI

Torun Radio Astronomy Observatory,
Nicolaus Copernicus University,
ul. Chopina 12/18, PL–87–100 Torun, Poland

ACCURATE ALGORITHMS TO TRANSFORM
GEOCENTRIC TO GEODETIC COORDINATES

 Abstract

The problem of the transformation is reduced to solving of the equation
 2 sin(ψ – Ω) = c sin 2ψ,
where Ω = arctg [bz/(ar)], c = (a2 – b2)/√{(ar)2 + (bz)2}, a and b are the semi-axes of the reference ellipsoid, and z and r are the polar and equatorial, respectively, components of the position vector in the Cartesian system of coordinates. Then, the geodetic latitude is found as φ = arctg [(a/b) tg ψ], and the height above the ellipsoid as h = (r – a cosψ) cosφ + (z – b sinψ) sinφ. Two accurate closed solutions are proposed of which one is approximate in nature and the other is exact. They are shown to be superior to others, found in literature and in practice, in both or either accuracy and/or simplicity.

The transformation of coordinates between the geodetic and Cartesian system is a matter of everyday practice in geodesy. While conversion of the geodetic coordinates, the latitude φ, the longitude λ and the height above a reference ellipsoid h, to the geocentric Cartesian x, y and z coordinates is done straightforwardly the reverse conversion poses certain difficulties. Usually it is being performed by successive approximations in an iterative process (e.g. Heiskanen and Moritz, 1967, Bartelme and Meissl, 1975, Urmaev, 1981, Hlibowicki, 1981), or using closed approximative formulae cast as series expansions (e.g. Morrison and Pines, 1961, Izotov et al., 1974, Long, 1976, Taff, 1985) or other expressions (e.g. Pick, 1985, Baranov et al., 1986). Even though there exist closed exact solutions for transformation of the Cartesian coordinates onto the geodetic coordinates (Paul, 1973, Hedgley, 1976, Heikkinen, 1982, Vaniček and Krakiwsky, 1982) many authors seem to believe to the contrary (e.g. Woolard and Clemence, 1966, Heiskanen and Moritz, 1967, Izotov et al., 1974, Taff, 1985, Baranov et al., 1986, Astronomical Almanac, 1987). The exact solutions are in principle error free, however, in practice errors arise due to rounding-off in relatively complicated algorithms. Nevertheless, in certain applications, they might be preferred to simpler approximative solutions both due to their higher accuracy and completeness (range of latitudes and heights).

In this report we offer two solutions: a very simple approximation algorithm, which gives accuracies well above any practical needs, and in fact approaching the precision of exact algorithms, and one truly exact, independent of the so far published. The derivation of the latter is so simple that it may serve as a standard reference. The simplicity comes from the fact that in our quartic (biquadratic) equation only two of the five terms have coefficients different from zero or ±1 (other authors solve a general form of the quartic equation). We shall also present a comparison of performance of a few of more popular algorithms, along with ours.

To start with, observe that the Cartesian coordinates can be obtained from the geodetic ones through the following formulae

 r = a cosψ + h cosφ     and (1a)
 z = b sinψ + h sinφ, (1b)

where a (b) is the semimajor (semiminor) axis of the reference ellipsoid and r, the equatorial component, is split into x and y by multiplying Eq.(1a) by the cosine and sine, respectively, function of the longitude. The ψ parameter, sometimes called parametric or reduced latitude, is calculated from

 ψ = arctg([b/a] tgφ). (2)

The Eqs. (1) are equivalent to the more common ones

 r = (N + h) cosφ (3a)
 z = [N(1 – e2) + h] sinφ, (3b)

where N = a/(1 – e2 sin2φ)1/2, and e2 = (a2 – b2)/a2. Now, given r and z, the problem of finding φ and h can be shown to reduce to solving of the following equation:

 2 sin(ψ – Ω) – c sin 2ψ = 0, (4)

where Ω = arctg [bz/(ar)] and c = (a2 – b2)/[(ar)2 + (bz)2]1/2. This equation obtains from geometrical considerations. Equivalently, it can be derived directly from Eqs. (1) by eliminating the h parameter. A ψ that solves Eq. (4), by virtue of Eq. (2), is used to calculate the geodetic latitude:

 φ = arctg[(a/b) tgψ]. (5)

This, in turn, allows the use of both or any one of the two equations (1) to find the height, e.g.

 h = (r – a cosψ) cosφ + (z – b sinψ) sinφ. (6)

We prefer this form for it has similar accuracy over all range of latitudes unlike the formulas met in practice and involving division by the sine or cosine of the latitude. Solving of the Eq. (4) should not present any difficulty in view of the fact that nowadays even handheld calculators have a built-in capability of finding the zeroes of transcendental functions of the type of the left side of (4). A convenient starting value for ψ can be chosen to equal arctg [az/(br)], which, as evidenced by Eqs. (1), is exact on the surface of the ellipsoid (with h = 0). Applying the Newton's method results in a rapidly converging series of approximations. This method uses the derivative of a function to be solved with respect to the variable (the ψ here), which in our case is

 w = 2[cos(ψ′ – Ω) – c cos 2ψ′], (7)

to find next approximation:

 ψ = ψ′ – [2sin(ψ′ – Ω) – c sin 2ψ′]/w, (8)

where the prime indicates the starting value or previous approximation. Computer tests have shown that just two iterations provide accuracies by far exceeding the needs of even very demanding practical user. The accuracies are of the order of 1E–9 rad in the geodetic latitude for all points further than a thousand of kilometers away from the origin of the Cartesian system. Another way of saying roughly the same is to claim that it is some six orders of magnitude more accurate than other algorithms met in literature, save the few truly exact solutions.

Table 1 presents errors of position calculated with the use of 10 different algorithms. The errors are defined to be the distances between the true position and the one calculated. To make the comparison meaningful the first three iterative algorithms (based on Heiskanen and Moritz (1967), Astronomical Almanac (1987), and Bartelme and Meissl (1975), respectively) were restricted to only two iterations (their performance for the third iteration is shown in brackets). The algorithms numbered 4 to 7 were taken from Izotov at al. (1974), Long (1976; Eq. (8) and (15) therein), Baranov et al. (1986) and Pick (1985), respectively. The Baranov et al. formula has been corrected for a misprint in the cited source (the a parameter in their expression for tgφ has been inadvertadly squared). It is clear from this table that our relatively simple algorithm, whose performance is listed in the column headed 8, compares favourably even with the exact solutions (the two rightmost columns in Table 1, of which the first corresponds to the algorithm given by Heikkinen (1982) and the second is ours, still to be presented, see also Borkowski 1987). Although our approximative solution is iterative in principle, due to its exceptionally fast convergence, all the process of iterations is reduced to two formulas which are essentially the first (Eq. (8) with ψ′ = Ω) and the second approximation (Eq. (8) itself). With such an understanding this algorithm can be safely treated also as a closed solution. It was used for calculation of entries to Table 1.

Eq.(4), when expressed in tg(π/4 – ψ/2) ≡ t, leads to the fourth-degree polynomial (biquadratic or quartic equation):

 t4 + 2Et3 + 2Ft – 1 = 0, (9)

where

 E = [bz – (a2 – b2)]/(ar)     and (10)
 F = [bz + (a2 – b2)]/(ar). (11)

Applying the Ferrari's solution (e.g. Korn and Korn, 1961) to Eq. (9) one obtains:

 t = ±√{G2 + (F – vG)/(2G – E)} – G, (12)

where

 G = [±(E2 + v)1/2 + E]/2, (13)
 v = (D1/2 – Q)1/3 – (D1/2 + Q)1/3, (14a)
 D = P3 + Q2, (15)
 P = (4/3)(EF + 1)    and (16)
 Q = 2(E2 – F2). (17)

For D < 0, to avoid the complex arithmetic, it is recommended to use the following equivalent of Eq.(14a)

 v = 2 √{–P} cos{(1/3)arccos[Q/(–P)3/2]}. (14b)

The case of D < 0 will be of marginal use in practice since D is a positive quantity anywhere further than about 45 km from the Earth center. Now, any real of the four solutions (corresponding to the four combinations of +/– signs in Eqs. (12) and (13)) represents a distinctive set of geodetic coordinates obtainable through:

 φ = arctg [a(1 – t2)/(2bt)]    and (18)
 h = (r – at) cosφ + (z – b) sinφ . (19)

As a practical matter, one does not need all the four solutions even though each of them satisfies Eqs. (1) (they are all real for D < 0 and two are real for D > 0). To obtain readily programmable algorithm yielding the single desired solution just skip out the double signs in Eqs. (12) and (13) to leave only positive square roots there. Then it will work provided a > b and φ > 0. Noteworthy, this algorithm will work properly also for southern latitudes (z < 0) if one additionally preassigns the z-coordinate sign to the minor semiaxis, b.

 Table 1 Position displacement errors (in mm) for ten algorithms. The numbers in brackets are for 3 iterations with procedures of Heiskanen and Moritz (1967), Astronomical Almanac (1987) and Bartelme and Meissl (1975), algorithms No. 1, 2 and 3, respectively. Other algorithms are based on approximate closed solutions: Izotov et al., (1974; No. 4), Long (1975; 5), Baranov et al., (1986; 6), Pick (1985; 7) and this paper (8), and exact solutions: Heikkinen (1982; 9) and this paper (10). The bottom line lists respective timing measurements in units of 0.1 ms.

 φ [deg] h [km] 1 2 3 4 5 6 7 8 9 10 89.000 100000 109(0) 0(0) 1(0) 721 190 1 925 .000000 .000015 .000000 89.000 10000 460(1) 0(0) 22(0) 721 158 0 126 .000004 .000002 .000000 89.000 1000 228(1) 0(0) 25(0) 721 857 0 0 .000000 .000000 .000000 89.000 0 0(0) 0(0) –(–) 721 1219 0 0 .000000 .000000 .000000 89.000 –1000 431(3) 0(0) 64(0) 721 1836 0 1 .000000 .000001 .000001 70.000 100000 89(0) 0(0) 0(0) 550 112 248 13722 .000021 .000017 .000000 70.000 10000 377(1) 1(0) 10(0) 635 105 105 2094 .000000 .000000 .000000 70.000 1000 187(1) 6(0) 11(0) 660 370 5 6 .000000 .000001 .000000 70.000 0 0(0) 9(0) –(–) 646 506 0 0 .000000 .000001 .000000 70.000 –1000 353(2) 12(0) 29(0) 611 734 10 34 .000000 .000000 .000001 45.000 100000 30(0) 1(0) 0(0) 167 24 452 10075 .000000 .000000 .000000 45.000 10000 128(0) 37(0) 0(0) 365 5 191 2683 .000000 .000000 .000000 45.000 1000 63(0) 181(1) 0(0) 587 215 9 14 .000000 .000001 .000001 45.000 0 0(0) 242(1) –(–) 636 327 0 0 .000000 .000000 .000000 45.000 –1000 120(0) 340(1) 0(0) 695 519 18 81 .000000 .000000 .000000 20.000 100000 2(0) 2(0) 0(0) 13 22 91 13641 .000007 .000000 .000015 20.000 10000 7(0) 73(0) 10(0) 118 25 38 2090 .000000 .000000 .000002 20.000 1000 3(0) 358(2) 11(0) 424 72 2 7 .000000 .000001 .000001 20.000 0 0(0) 478(3) –(–) 546 121 0 0 .000000 .000001 .000001 20.000 –1000 6(0) 672(5) 28(0) 736 205 4 35 .000000 .000001 .000000 1.000 100000 0(0) 0(0) 1(0) 0 1 0 922 .000015 .000015 .000009 1.000 10000 0(0) 5(0) 22(0) 5 1 0 126 .000000 .000000 .000000 1.000 1000 0(0) 25(0) 24(0) 25 10 0 0 .000000 .000001 .000000 1.000 0 0(0) 33(0) –(–) 33 17 0 0 .000000 .000000 .000000 1.000 –1000 0(0) 47(0) 63(0) 47 27 0 1 .000000 .000001 .000000 Exec. time: 56(78) 51(70) 50(63) 44 34 46 146 74 56 63

We have tested our exact algorithm, as embodied in Eqs. (10) through (19), over wide range of (r,z) plane (save the z-axis, where the quantities E and F are obviously indeterminate) without encountering any problems. As to the round-off errors, they are certainly not serious for practical applications (we have always used the DOUBLE PRECISION arithmetic in FORTRAN programming language), though we found them considerably larger in a close vicinity of the r and z axes than elsewhere. The main source of these errors belongs to evaluation of Eq. (14a), which involves the difference of cubic roots (computed as the exponent of one third of the natural logarithm). The simple way to circumvent this problem is to correct a result obtained from Eq. (14a), say v′, by using the cubic resolvent (of which v′ is one of the three roots):

 v = –(v′3 + 2Q)/(3P). (20)

It can be easily proved that v from the above equation is more accurate than v′ if v2 < P, which is the case everywhere outside the radial distance from the Earth's centre of about 70 km. The solution of Heikkinen (1982), which also involves the cubic root operation, is free from the considered weakness, however, it is not designed to handle the case of D < 0. For D > 0 and with the correction (20) in our algorithm the two solutions can be considered equally good.

To compare the effectiveness of the discussed algorithms we have measured the amount of CPU time required on an IBM PC type machine using the Lahey Computer Systems FORTRAN compiler. The results, given in units of 0.1 ms in the last row of Table 1, show that all times are comparable despite considerable differences in complexity of the algorithms. This is so, because the main contribution to the computation time is due to the trigonometric functions and they are a few in number in each algorithm. For example, the evaluation of the cosine function takes about 0.6 ms, or 6 units of Table 1. Since the computing is really cheap, no special attempt was undertaken to optimize the execution times. Also, sizes of the algorithms in our implementation are comparable and range from 7 to 15 FORTRAN cards, with the single exception of Pick's solution, which anyway seems to be the worst in every respect.

To conclude, we find that the iterative algorithms, such like these published by Heiskanen and Moritz (1967) or in Astronomical Almanac (1987), are the simplest to program and give acceptable results with 3 to 4 iterations. With one or two additional iterations they can also be used when higher accuracy is required. However, due to round-off errors attempting more than 6 iterations may not improve the result obtained so far, which may be slightly worse than the one obtainable with algorithms based on exact solutions. In this case the two algorithms presented in this paper, or the one of Heikkinen (1982), seem to be more appropriate. They guarantee accuracies a few orders below 1 mm over all practically useful area on the (r,z)-plane, from as close as a few tens of kilometers from the very centre of the Earth to many Earth's radii out into the space. For that matter the correction (20) to our mathematically exact solution may be safely omitted. For more general solutions of the transformation between the Cartesian and geodetic coordinates the only readily available algorithm is our exact solution (presented originally in full detail, unfortunately with an annoying misprint in the equation corresponding to Eq. (12), in Borkowski (1987)).

FORTRAN listings of all algorithms analysed in this report may be requested from the author.

Acknowledgements

I would like to thank dr B . Benciolini for helpful comments on an earlier version of this report. This work was financially supported through the Polish Government research problem RPBR No. RR.I.11/2.

 REFERENCES

Astronomical Almanac for the Year 1987, USNO and RGO, Washington and London, p. K11.

V.N. BARANOV et al., 1986, Space Geodesy, Nedra, Moskva (in Russian), p. 27.

N . BARTELME, P . MEISSL, 1975, Ein einfaches, rasches und numerisch stabiles Verfahren zur Bestimmung des kürzesten Abstandes eines Punktes von einen sphäroidischen Rotationsellipsoid, AVN, Heft 12, p. 436–439 (Wichmann, Karlsruhe) (in German).

W. BENNING, 1974, Der Kürzeste Abstand eines in rechtwinkeligen Koordinaten gegebenen Aussenpunktes vom Ellipsoid, AVN, Heft 11, (Wichmann, Karlsruhe) (in German).

K.M. BORKOWSKI, 1987, Transformation of Geocentric to Geodetic Coordinates without Approximations, Astrophys. Space Sci., 139, 1–4.

D.R. HEDGLEY, 1976, An Exact Transformation from Geocentric to Geodetic Coordinates for Nonzero Altitudes, NASA Techn. Rep. R–458.

M. HEIKKINEN, 1982, Geschlossene Formeln zur Berechnung raumlicher geodätischer Koordinaten aus rechtwinkligen Koordinaten, Zeitschrift Vermess. (in German), 107, 207–211 .

W.A. HEISKANEN and H. MORITZ, 1967, Physical Geodesy, W.H. Freeman and Co., San Francisco, p. 181.

R.HLIBOWICKI (ed.), 1981, Higher Geodesy and Geodetic Astronomy, PWN, Warszawa (in Polish), p. 273.

A.A. IZOTOV et al., 1974, Fundamentals of Satellite Geodesy, Nedra, Moskva (in Russian), p. 34.

G.A. KORN and T.M. KORN, 1961, Mathematical Handbook for Scientists and Engineers, McGraw–Hill, New York, p. 24.

S.A.T. LONG, 1975, General Altitude Transformation between Geocentric and Geodetic Coordinates, Celestial Mech., 12, 225–230.

J. MORRISON and S. PINES, 1961, The Reduction from Geocentric to Geodetic Coordinates, Astron. J., 66, 15–16.

M.K. PAUL, 1973, A Note on Computation of Geodetic Coordinates from Geocentric (Cartesian) Coordinates, Bull. Géod., Nouv. Ser., No 108, 135–139.

M . PICK, 1985, Closed Formulae for Transformation of the Cartesian Coordinate System into a System of Geodetic Coordinates, Studia geoph. et geod., 29, 112–119.

L.G. TAFF, 1985, Celestial Mechanics. A Computational Guide for the Practitioner, J . Wiley and Sons, New York, ch. 3.

M.S. URMAEV, 1981, Orbital Methods of Space Geodesy, Nedra, Moskva (in Russian), p . 14.

P . VANIČEK and E.J. KRAKIWSKI, 1982, Geodesy: The Concepts, North Holland Publ. Co., Amsterdam, p . 324.

E.W. WOOLARD and G.M. CLEMENCE, 1966, Spherical Astronomy, Academic Press, New York, ch. 3.