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