Just to check that I'm on the right track:

Are all great-circles on the sphere and in equirectangular projection (i.e. latitude, longitude pairs) either:

  1. meridians (i.e. going pole to pole)
  2. of the form tan latitude = sin360(longitude + rotation) * amplitude + offset

(with additional restrictions on the offset/amplitude combinations -- obviously, all great circle paths with an amplitude of 0 also have offset 0 - the equator).

Or are there great-circle paths that do not fit into this scheme (again, only in a longitude-latitude coordinate system, not on other map projections).

Note: I added the tan above after posting the question, in a reaction to whubers excellent reply. It turns out that offset then always is 0.

  • Form #2 is not a set of (latitude, longitude) pairs. What does it mean? And why does it depend on three parameters (rotation, amplitude, and offset) when the great circles naturally have only two parameters (each one corresponds to a pair of diametrically opposite points that are "polar" to it)? – whuber Feb 26 '13 at 17:49
  • Yes, one is redundant, but I don't have the formula ready. As noted, amplitude==0 implies offset=0; these two are obviously coupled. See updated question for the missing latitude to make the relationship well-formed. – Erich Schubert Feb 26 '13 at 17:53

Although geodesics do look a little like sine waves in some projections, the formula is incorrect.

Here is one geodesic in an Equirectangular projection. Clearly it is not a sine wave:

enter image description here

(The background image is taken from http://upload.wikimedia.org/wikipedia/commons/thumb/e/ea/Equirectangular-projection.jpg/800px-Equirectangular-projection.jpg.)

Because all Equirectangular projections are affine transformations of this one (where the x-coordinate is the longitude and the y-coordinate is the latitude), and affine transformations of sine waves are still sine waves, we cannot expect any geodesics in any form of the Equirectangular projection to be sine waves (except for the Equator, which plots as a horizontal line). So let's start at the beginning and work out the correct formula.

Let the equation of such a geodesic be in the form

latitude = f(longitude)

for a function f to be found. (This approach has already given up on the meridians, which cannot be written in such a form, but otherwise is fully general.) Converting to 3D Cartesian coordinates (x,y,z) gives

x = cos(l) cos(f(l))
y = sin(l) cos(f(l))
z = sin(f(l))

where l is the longitude and a unit radius is assumed (without any loss of generality). Since geodesics on the sphere are intersections with planes (passing through its center), there must exist a constant vector (a,b,c)--which is directed between the poles of the geodesic--for which

a x + b y + c z = 0

no matter what the value of l might be. Solving for f(l) gives

f(l) = ArcTan(-(a cos(l) + b sin(l)) / c)

provided c is nonzero. Evidently, when c approaches 0, we obtain in the limit a pair of meridians differing by 180 degrees--precisely the geodesics we abandoned at the outset. So all is good. By the way, despite appearances, this uses just two parameters equal to a/c and b/c.

Note that all geodesics can be rotated until they cross the equator at zero degrees longitude. This indicates that f(l) can be written in terms of f0(l-l0) where l0 is the longitude of the equatorial crossing and f0 is the expression of a geodesic crossing at the Prime Meridian. From this we obtain the equivalent formula

f(l) = ArcTan(gamma * sin(l - l0))

where -180 <= l0 < 180 degrees is the longitude of the equatorial crossing (as the geodesic enters the Northern Hemisphere when traveling eastwards) and gamma is a positive real number. This does not include the meridian pairs. When gamma = 0 it designates the Equator with a starting point at longitude l0; we may always take l0 = 0 in that case if we wish a unique parameterization. There are still just two parameters, given by l0 and gamma this time.

Mathematica 8.0 was used to create the image. In fact, it created a "dynamic manipulation" in which the vector (a,b,c) can be controlled and the corresponding geodesic is instantly displayed. (That's pretty cool.) First we obtain the background image:

i = Import[

Here is the code in its entirety:

 {a, b, c} = {Cos[u] Cos[v], Sin[u] Cos[v], Sin[v]};
    Polygon[{{-\[Pi], -\[Pi]/2}, {\[Pi], -\[Pi]/2}, {\[Pi], \[Pi]/2}, {-\[Pi], \[Pi]/2}}, 
     VertexTextureCoordinates -> {{0, 0}, {1, 0}, {1, 1}, {0, 1}}]}], 
  Plot[ArcTan[(a Cos[\[Lambda]] + b Sin[\[Lambda]]) / (-c)], {\[Lambda], -\[Pi], \[Pi]}, 
   PlotRange -> {Automatic, {-\[Pi]/2, \[Pi]/2}}, PlotStyle -> {Thick, Red}]],
   {u, 0, 2 \[Pi]}, {v, -\[Pi]/2, \[Pi]}]
  • Thanks. I lost that arctan somewhere. Guess I started with arctan latitude somewhere. – Erich Schubert Feb 27 '13 at 11:16
  • I wish I could up vote this more than once! – Ian Turton Feb 27 '13 at 12:46

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.