Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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

share|improve this question
up vote 21 down vote accepted

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


(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 )
share|improve this answer
+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
Yes, thank you Kirk. I'll fix that. – whuber Feb 13 '12 at 22:08
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. – Dan Patterson 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

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.

share|improve this answer
Can you provide a link to the content? – Trevor Boyd Smith Feb 13 '12 at 15:43

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); 
share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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