Information Technology
ENContact

Great Circle on Spherical Earth

great circle

The shortest path between two points on the (spherical) Earth lies on the great circle going through these two points. The great circle is - as the name suggests - the largest possible circle on the sphere. All great circles of a given sphere have both the same radius and the same center as the sphere. The great circle's plane passes through the center of the sphere and divides it into 2 equal hemispheres.

Try the interactive applet.

Great Circle between Two Points

The equation of the great circle is computed most easily in cartesian coordinates. So, if the points were given in spherical coordinates we would transform them to cartesian coordinates first.

Let

p1 = (x1, y1, z1)T
(1a)
p2 = (x2, y2, z2)T
(1b)

be the vectors corresponding to the two points P1 and P2, respectively. The great circle is located in the plane OP1P2 spanned by these two vectors.

If we find two orthonormal vectors u and v in this plane then the equation of the great circle will be

c = r (u cos ω + v sin ω)
(2)
0° ≤ ω < 360°

The obvious choice for the vector u is the normalized vector p1:

u = p1 / R.
(3)

The vector product

w = p1 x p2
(4)

is orthogonal to the plane OP1P2. The unit vector

v = u x w / |w|.
(5)

is in the plane OP1P2 and is orthogonal to u.

Another way to compute the vector v is by using the Gram-Schmidt orthonormalization process. The angle δ between p1 and p2 can be computed using the scalar product:

p1 · p2 = R2 cos δ.
(6)

The projection of p2 to u is:

pu = R cos(δ) u.
(7)

The vector

pv = p2 - pu
(8)

is orthogonal to p1 and is in the plane OP1P2. Hence the normalized vector v is our searched vector:

v = pv / |pv| 
(9)

The above calculations break if the vectors p1 and p2 are parallel, because then they do not specify a plane. This happens if the corresponding points are either identical or antipodal. In this case there are infinite many great circles through P1 and P2. Their equations are still given by

c = r (u cos ω + v sin ω)
(10)

where v is an arbitrary unit vector orthogonal to u.

Shortest Distance between Two Points

The shortest path between two points on the earth surface is the circular arc on the great circle between these points. So, the shortest distance d is

d = R δ
(11)

where R is the earth radius and δ is the angle between the two points. From the equation (6) it follows:

cos δ = p1 · p2 / R2 = (x1x2 + y1y2 + z1z2) / R2. 
(12)

By substituting spherical coordinates and using the identity:

cos(λ1 - λ2) = cosλ1 cosλ2 + sinλ1 sinλ2
(13)

we get:

cos δ = cos(λ1 - λ2) cos φ1 cos φ2 + sin φ1 sin φ2 
(14)

Finally, the shortest distance d is:

d = R δ = R acos[cos(λ1 - λ2) cos φ1 cos φ2 + sin φ1 sin φ2]. 
(15)

Heading (Bearing, Azimuth)

Heading (also called bearing or azimuth) is the angle α between a great circle and a meridian. Heading plays an important role in navigation. When we travel from the point P1 to the point P2 on the shortest path, then the heading at point P1 specifies the starting direction of our journey. The heading changes during the voyage, of course.

Heading as angle between planes

Heading α is the angle between the great circle plane OP1P2 and the meridian plane OP1P2, or, equivalently between the respective normal vectors of these two planes. The vector orthogonal to the meridian plane is

m = p1 x (0, 0, 1)T = (y1, -x1, 0)T 
(16)

The vector orthogonal to the path plane is w = (wx, wy, wz)T = p1 x p1.

The dot product of the normalized vectors w and m equals the cosine of the heading angle α:

cosα = w/|w| · m/|m|
(17)

When we can take arccos() we get the angle α, 0 ≤ α ≤ 180°. The usual convention is that a positive angle denotes travelling eastwards while a negative angle denotes travelling westwards. Under the assumption that we travel along the shorter arc from P1 to P2 we can determine the direction by considering the vector product w. When it points upwards with respect to the positive z-axis the journey goes eastwards, when it points downwards the journey goes westwards. So, our formula for the heading angle is finally:

α =  acos[cos(α)] if wz ≥ 0 
(18a)
α = -acos[cos(α)] if wz < 0 
(18a)

Simplification of the geometry

To simplify the following calculations I will make two assumptions without loss of generality:

  1. The radius of the earth is 1. This is justified because only angles are considered.
  2. The point P1 has the longitude 0 and the point P2 has the longitude = λ2 - λ1. I just start counting the longitude at the meridian passing through P1. The vector p1 is hence (x1, 0, z1)T.

Under the simplifying assumptions we get:

m = (0, -x1, 0)T 
(19)

and

wx =  z1 y2
(20a)
wy = -x1 z2 + z1 x2
(20b)
wz =  x1 y2
(20c)

and

cosα = -wy / |w|
(21)

Heading as angle in a spherical triangle

spherical triangle

The heading is the angle at the vertex P1 of the spherical triangle NP1P2, where N is the north pole. The sides of this triangle are arcs of the respective great circles and their lengths are equal to the respective angles at the center of the Earth.

Now we can apply the cosine rule formula for the side P2N that is opposite to the angle α (see the equation (10) in Weisstein, Eric W. Spherical Trigonometry in MathWorld):

cosΦ2 = cosΦ1 cosδ + sinΦ1 sinδ cosα
(22)

Solving for cosα and taking into account that

Φ1 = 90° - φ1 
(23a)
Φ2 = 90° - φ2
(23b)

we get:

cosα = 
sinφ2 - sinφ1 cosδ 

cosφ1 sinδ
(24)

Tangens formula

Let us have a look at the geometry of the vector w:

heading tangens

In the equation (21) we computed in fact α using the cosine formula in the triangle OBC. Now, we will compute it using the tangens formula:

tanα = - wxz / wy 
(25)

Here, wxz is the projection of the vector w onto the XZ-plane:

wxz2 = wx2 + wz2 = z12 y22 + x12 y22 = y22(x12 + z12) = y22 
(26)

The last equality follows from the fact that |p1| = 1 and y1 = 0. Substituting (26) and (20b) into (25) we get:

tanα = 
y2

x1 z2 - z1 x2 
(27)

Finally, using spherical coordinates we obtain:

tanα = 
cosφ2 sin 

cosφ1 sinφ2 - sinφ1 cosφ2 cos 
(28)

Even if this formula looks more complicated than the formula (24) it is often recommended because it avoids the calculation of the angle δ between the two points P1 and P2. However, significant care must be taken by applying atan() to the equation (28). A simpler way is to use atan2() instead of atan():

  dl = lon2 - lon1
  y2 = sin(dl)*cos(lat2)
  wy = cos(lat1)*sin(lat2) - sin(lat1)*cos(lat2)*cos(dl)
  alpha = atan2(y2, wy)

References