Computing Altitude from an Earth-Centered Position
Date: 01/04/2002 at 13:09:26 From: Steve Subject: Computing Altitude from Earth-Centered I'm looking for a closed form algorithm which provides altitude above the reference (WGS-84) ellipsoid given an Earth-centered position. Thanks.
Date: 01/06/2002 at 14:58:17 From: Doctor Fenton Subject: Re: Computing Altitude from Earth-Centered Hi Steve, Thanks for writing to Dr. Math. To compute altitude from ECEF xyz-coordinates, there are closed-form solutions and iterative solutions. The best closed-form solution seems to be that of K. M. Borkowski, a Polish radio-astronomer. It requires finding a root of a fourth-degree polynomial, and is therefore possible using the classical fourth-degree formula for roots. The algorithm is given in detail in the (new) Explanatory Supplement for the Astronomical Almanac (P. Seidelmann, editor). There's a link to an abstract of Borkowski's paper on his personal Web page: http://www.astro.uni.torun.pl/~kb/personal.html The abstract is at http://www.astro.uni.torun.pl/~kb/abstract.html#Transf2 A library should be able to obtain the article for you: Bulletin Geodesique, vol. 63, no. 1, 1989, pp. 50-56 or perhaps obtain or borrow the Explanatory Supplement for you. Most solutions, however, are iterative. The standard algorithm from the Astronomical Almanac is the following: Given ECEF coordinates x, y, and z, compute r = sqrt(x^2 + y^2) e^2 = f*(2-f) (where 1/f = 298.257224) and phi0 = arctan(z/r) . Initialize phi to phi0; then iterate until phi is unchanged (usually no more than 4 or 5 iterations): phi1 = phi C = 1/sqrt(1 - e^2*sin^2(phi1)) phi = arctan((z + a*C*e^2*sin(phi1))/r) . When the iteration has converged, r h = -------- - aC . cos(phi) I tested this in an Excel spreadsheet, and while it takes 5 or 6 iterations to give the correct geodetic latitude, the height above the reference ellipsoid was correct to within a few centimeters on the first iteration, and exact on the second iteration. Another algorithm due to B. R. Bowring seems to be even more accurate. The longitude is found in the usual way. To compute the latitude and altitude, given ECEF coordinates x, y, and z, compute p = sqrt(x^2 + y^2) and r = sqrt(p^2 + z^2) . Then compute u = atan((z/p)*[(1-f)+e^2*a/r], and lat = atan((z*(1-f)+e^2*a*(sin(u))^3)/((1-f)*(p-e^2*a*(cos(u))^3))) and h = p*cos(lat) + z*sin(lat) - a*sqrt(1-e^2*sin(lat)^2) . On my spreadsheet, this gave correct altitudes (and latitudes) on the first iteration. The reference I have is Bowring, B.R. (1985). "The accuracy of geodetic latitude and height equations," Survey Review. 28 Oct., pp. 202-206. If you have trouble finding the Borkowski algorithm and do not find these algorithms sufficient, please write back and I will look up the details of Borkowski's method for you. - Doctor Fenton, The Math Forum http://mathforum.org/dr.math/
Date: 01/07/2002 at 14:10:58 From: Steve Subject: Computing Altitude from Earth-Centered Thanks for the great answer... concise, complete and correct. What more could I ask for? The iterative approach works just fine for me. A tweak: using phi0 = arctan(z/((1-e^2)*r) generates a more accurate first-order approximation for latitude. In general this reduced my iteration count by one or two in most cases. For performance reasons you could also include a in the computation of C and remove this multiplication from its other occurrences. I also found that I needed to recompute C based on the final value for phi which came out of the iterations. Thanks for the algorithm! I hope that my suggestions help as well. Steve O'Neill
Search the Dr. Math Library:
Ask Dr. MathTM
© 1994- The Math Forum at NCTM. All rights reserved.