coordinate systemcoordinatesecef"latitude longitude
I wish to implement this formula from Wikipedia to convert LLA to ECEF
coordinates:
where
and a and b are the semi-major and semi-minor axes respectively.
My question is: what units will this work with? This paper says the equation is in meters. I am wondering if it can be more general than that. Can I use any unit of distance so long as I am consistent?
This question assumes an ellipsoidal model of the earth. Its reference surface is obtained by rotating an ellipse around its minor axis (plotted vertically by convention). Such an ellipse is just a circle that has been stretched horizontally by a factor of a and vertically by a factor of b. Using the standard parameterization of the unit circle,
t --> (cos(t), sin(t))
(which defines cosine and sine), we obtain a parameterization
t --> (a cos(t), b sin(t)).
(The two components of this parameterization describe a trip around the curve: they specify, in Cartesian coordinates, our location at "time" t.)
The geodetic latitude, f, of any point is the angle that "up" makes to the equatorial plane. When a differs from b, the value of f differs from that of t (except along the equator and at the poles).
In this picture, the blue curve is one quadrant of such an ellipse (greatly exaggerated compared to the earth's eccentricity). The red dot at the lower left corner is its center. The dashed line designates the radius to one point on the surface. Its "up" direction there is shown with a black segment: it is, by definition, perpendicular to the ellipse at that point. Due to the exaggerated eccentricity, it is easy to see that "up" is not parallel to the radius.
In our terminology, t is related to the angle made by the radius to the horizontal and f is the angle made by that black segment. (Note that any point on the surface can be viewed from this perspective. This allows us to limit both t and f to lie between 0 and 90 degrees; their cosines and sines will be positive, so we don't have to worry about negative square roots in the formulas.)
The trick is to convert from the t-parameterization to one in terms of f, because in terms of t the radius R is easy to compute (via the Pythagorean theorem). Its square is the sum of squares of the components of the point,
R(t)^2 = a^2 cos(t)^2 + b^2 sin(t)^2.
To make this conversion we need to relate the "up" direction f to the parameter t. This direction is perpendicular to the tangent of the ellipse. By definition, a tangent to a curve (expressed as a vector) is obtained by differentiating its parameterization:
Tangent(t) = d/dt (a cos(t), b sin(t)) = (-a sin(t), b cos(t)).
(Differentiation computes the rate of change. The rate of change of our position as we travel around the curve is, of course, our velocity, and that always points along the curve.)
Rotate this clockwise by 90 degrees to obtain the perpendicular, called the "normal" vector:
Normal(t) = (b cos(t), a sin(t)).
The slope of this normal vector, equal to (a sin(t)) / (b cos(t)) ("rise over run"), is also the tangent of the angle it makes to the horizontal, whence
tan(f) = (a sin(t)) / (b cos(t)).
Equivalently,
(b/a) tan(f) = sin(t) / cos(t) = tan(t).
(If you have good insight into Euclidean geometry, you could obtain this relationship directly from the definition of an ellipse without going through any trig or calculus, simply by recognizing that the combined horizontal and vertical expansions by a and b respectively have the effect of changing all slopes by this factor b / a.)
Look again at the formula for R(t)^2: we know a and b -- they determine the shape and size of the ellipse -- so we only need to find cos(t)^2 and sin(t)^2 in terms of f, which the preceding equation lets us do easily:
cos(t)^2 = 1/(1 + tan(t)^2)
= 1 / (1 + (b/a)^2 tan(f)^2)
= a^2 / (a^2 + b^2 tan(f)^2);
sin(t)^2 = 1 - cos(t)^2
= b^2 tan(f)^2 / (a^2 + b^2 tan(f)^2).
(When tan(f) is infinite, we're at the pole, so just set f = t in that case.)
This is the connection we need. Substitute these values for cos(t)^2 and sin(t)^2 into the expression for R(t)^2 and simplify to get
R(f)^2 = ( a^4 cos(f)^2 + b^4 sin(f)^2 ) / ( a^2 cos(f)^2 + b^2 sin(f)^2 ).
A simple transformation shows that this equation is the same as the one found on Wikipedia. Because a^2 b^2 = (ab)^2 and (a^2)^2 = a^4,
R(f)^2 = ( (a^2 cos(f))^2 + (b^2 sin(f))^2 ) / ( (a cos(f))^2 + (b sin(f))^2 )
You can compare the two. In most applications I suspect the second (direct) method will be the one to choose.
Accuracy of the first (iterative) method depends on the accuracy with which you do the computations and when you decide to stop iterating. It therefore can be made as accurate as the second method for all inputs where both are valid (the first method works only for terrestrial heights, not for astronomical ones).
Which is faster depends on the programming environment, the computing architecture, and how much accuracy is needed. (In my tests with Mathematica, the second--direct--method is actually twice as fast as the first, practically no matter how much inaccuracy was tolerated in the iterative method.) Because both methods take about the same amount of computation but the first has to iterate at least once, it actually can be slower. If you're only doing the conversions at sea level (h=0), the iterative method might be a little faster, but the difference isn't going to be huge (I would be surprised by a twofold advantage).
BTW, note that the "closed formula" (for the second method) is a little deceptive: when you get to computing the height h, you need to obtain the radius of curvature N in terms of the latitude you have just calculated (phi). To do that, use the formula for N found in the preceding section.
Best Answer