19

(I see there is an equation on wikipedia that does exactly what I am asking but there is no references. I have no way of confirming the validity of this equation!)

I already understand the difference between Geocentric Latitude vs Geodetic Latitude.

Assuming known semi major, a, and semi minor, b, radiuses are given. How do you compute the radius at a given geodetic latitude?

I need some sort of expert confirmation (derivation, link to derivation, confirmation from expert, explanation, etc).

25

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).

Figure

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 )
  • +1 .. except I think the final formula has a paren out of place ... shouldn't (b^4 sin(f))^2 be changed to (b^4 sin(f)^2) ? – Kirk Kuykendall Feb 13 '12 at 21:59
  • really glad there are some experts around on this subject =). – Trevor Boyd Smith Feb 14 '12 at 18:41
  • Can a Geogebra file (html) be posted to this site? I have a radius of the prime vertical that might visually demonstrate what is going on. – user681 Feb 14 '12 at 19:20
  • You can export the original in .png format, @Dan: use the File|Export dialog. I recommend using large fonts (16 or 18 points seem to work well) and zooming as far in to the image as you can first. – whuber Feb 14 '12 at 19:30
  • I assume that the interactivity will be lost then. The demo demonstrates how varying the radii and the latitude of interest changes the properties. – user681 Feb 14 '12 at 19:47
3

Interesting to find that my math illiterate solution did the job with 5 minutes of thought and coding, wouldn't the flattening factor have to be considered rather than a perfect elliptical model?

        double pRad = 6356.7523142;
        double EqRad = 6378.137;                      
        return pRad + (90 - Math.Abs(siteLatitude)) / 90 * (EqRad - pRad); 
  • 1
    Where pRad is the polar radius, and EqRad is the equatorial radius. – Stefan Steiger Oct 9 '17 at 8:23
  • this is the only answer I could read. It seems to work for me. – Sean Bradley May 17 '18 at 19:52
  • 1
    I see that you are doing a linear interpolation of radius between pole and equator. While there is no reason to believe that linear interpolation is accurate, I'll use this as "good enough" for the Earth, given its mild flattening factor. BTW I think its slightly easier to read the equivalent: return E + (P - E) * Abs(Lat) / 90, so don't have to have 90 - ... in the formula. – ToolmakerSteve Nov 24 '18 at 20:03
2

enter image description here

At least that's the formula I found at the US Data Analysis and Assessment Center (DAAC) for the Department of Defense (DoD) High Performance Computing Modernization Program (HPCMP) wiki. It does say that they borrowed heavily from Wikipedia's entry. Still, the fact that they retained that formula should count for something.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.