Why Bowring 1976 for ecef2geodetic?

I'm not sure if this is the right place to ask this question, but I thought it would be a good first attempt. I am wondering why in the Matlab function ecef2geodetic.m which then calls cylindrical2geodetic.m, Matlab chooses to use Bowring's method from his 1976 paper when Bowring published an improved version in 1985. In the 1985 version he uses a different initial guess and gives an improved calculation for height. So I'm curious why, if you choose Bowring, you're not using his best solution. Also curious why choose Bowring in the first place over other methods like Borkowski which tend to run faster because there are fewer trig function calls.

Bowring. (1985). The Accuracy of Geodetic Latitude and Height. Survey Review, 28(218):202-206.

[James Tursa]

Maybe you could post this as a potential improvement to TMW Support.

Answer 1

[Rob Comer]

Good questions. You asked about (1) an improved calculation for height, (2) a different initial guess, and (3) performance and other (non-iterative) methods that potentially require fewer trigonometric function calls.

(1) Height Calculation

Actually, the Mapping Toolbox function ecef2geodetic has always used the height formula from Bowring (1985) rather than the height formula from the 1976 paper. Here's where the toolbox implementation came from.

The purpose of ecef2geodetic is to convert a point in Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates to a point in geodetic coordinates (geodetic latitude, longitude, and height above the reference ellipsoid). That's done using a fixed-point iteration to invert the much more straightforward geodetic-to-ECEF conversion (implemented in geodetic2ecef). The geodetic-to-ECEF conversion, from to , can be written as

where

is the distance to the point from the polar axis, and

.

(And where a and e are the semimajor axis and eccentricity of the ellipsoid, respectively.)

Bowring's 1976 height formula (applied after iterating to determine ϕ) is

.

It follows from solving for h in the expression for p above. At the time ecef2geodetic was introduced, in R2012b, it was clear that this formula is unstable if is too close to , although it's fine for ϕ close to 0. Something more had to be done.

Fortunately, one can also solve for h in the expression for z,

.

This formula is stable when is close to or equal to , but not for ϕ near 0. So, the two formulas for h are complementary. If we combine them with a pair of nonzero weights, and , such that , vanishes rapidly as ϕ approaches , and vanishes rapidly as ϕ approaches 0, we'll have a stable formula. A natural choice is and The result is

.

This stable formula is what's always been used in ecef2geodetic (more recently in a slightly refactored form, to improved the odds of a perfect "roundtrip" between the two conversion functions). It's easy to see its equivalence to the formula in Bowring (1985),

.

Rather than considering a weighted combination, however, Bowring explained his formula with a geometric argument.

(2) Initial Guess

In R2021a the initial guess was updated to use the starting values from Bowring (1985). That was a small change. It helped reduce the number of iterations needed in most cases, but by only a fraction of an iteration on average. A few cases actually required an additonal iteration after the change. But none of this is very significant; ecef2geodetic typically converges in 3 iterations or less.

(3) Performance and Non-Iterative Methods with Fewer Trigonometric Calls

Changes to the detailed implementation in R2021a provided measurable perforance improvements by eliminating nealy all trigonometric function calls (which was easier than one might think). I'm unaware of any comparison between ecef2geodetic and any of the many non-iterative methods that have been published over the years.