How to get the cardinal direction from one location to another?

Given are two geo locations, each with latitude and longitude.

One is the current location, the other is a target location.

is there a formula for calculating the target's cardinal direction for 0 being North, 90 being East and so on?

Edit: Needed is the direction for the shortest distance. It isn't necessary to give a solution for two opposite locations (for example the poles).

-
Given two antipodal points, say on the equator, how are you going to decide in a consistent manner what the direction is from one to the other? What about if the current location is one of the poles? –  Zev Chonoles Jun 14 '11 at 13:34
Good point, added it to the question. –  Elias Jun 14 '11 at 13:44
Do you want the direction along the shortest path to the target or along the line of constant heading? –  Rahul Jun 14 '11 at 15:46
The reason I bring this up is because getting directions out of great circles can be counterintuitive: the direction from London to Anchorage, Alaska is north, for example, and so is the direction from Anchorage to London. –  Rahul Jun 14 '11 at 21:47

Consider the initial point to be $(\theta_1,\phi_1)$ and the final point to be $(\theta_2,\phi_2)$. The angle $\theta$ is chosen between $-\pi/2$ (South) and $+\pi/2$ (North), with $0$ being the equator, because that's how you wanted it. The coordinates on the unit sphere are $$V_{i} =(\cos\theta_i\cos\phi_i,\cos\theta_i\sin\phi_i,\sin\theta_i), \qquad i=1,2$$ where the angles should have the right indices $1$ or $2$ for the two points. Without a loss of generality, we may assume $\phi_1=0$ - we rotate the Earth around the pole-pole axis - and $\phi_2$ will represent what we originally called $\phi_2-\phi_1$; at the end, we will substitute $\phi_2-\phi_1$ for $\phi_2$ everywhere.

The shortest path between the points is along the maximal circle. It's useful to quantify the direction of the plane in which this maximal circle lies. Its normal direction the cross product of $V_1$ and $V_2$, i.e. $N = V_1\times V_2$. Note that this cross product is degenerate (zero) if the two points are antipodal i.e. if $\phi_2=\pi$ and $\theta_2=-\theta_1$. $$N = (-\cos\theta_2 \sin\phi_2 \sin\theta_1, \cos\phi_2 \cos\theta_2 \sin\theta_1 - \cos\theta_1 \sin\theta_2, \cos\theta_1 \cos\theta_2 \sin\phi_2 )$$ The coordinates belonging to the maximum circle are given by $N\cdot (x,y,z) = 0$ where $(x,y,z)$ is a point on the unit sphere. Needless to say, $(x,y,z)=V_i$ satisfy this condition both for $i=1$ and $i=2$.

This was really unnecessary but it could be useful for other tasks. Now, let's calculate the tangential vector $T$ along the unit sphere from the first point $V_1$ that points to the second point $V_2$. It has to be in the $V_1$-$V_2$ plane, i.e. orthogonal to $N$, but also orthogonal to $V_1$. The right solution is $T=N\times V_1$. The vector is

$$T = (\sin\theta_1 (\cos\phi_2\cos\theta_2 \sin\theta_1 - \cos\theta_1 \sin\theta_2), \cos\theta_2 \sin\phi_2, \cos\theta_1 (-\cos\phi_2 \cos\theta_2 \sin\theta_1 + \cos\theta_1 \sin\theta_2))$$

Now, the cardinal direction $\alpha$ is calculable from this vector easily if we rotate the vector above in the $xz$-plane by $\theta_1$ so that the original point would get mapped to the $(1,0,0)$ point at the equator. So we want to evaluate the vector $$U = (t_1\cos \theta_1 + t_3\sin\theta_1, t_2, t_3\cos\theta_1-t_1\sin\theta_1 )$$ and the value turns out to be simply $$U = (0, \cos\theta_2\sin\phi_2, -\cos\phi_2 \cos\theta_2 \sin\theta_1 + \cos\theta_1 \sin\theta_2 )$$ As expected, the first coordinate vanishes because the vector has to be orthogonal to $(1,0,0)$. The remaining two components determine the cardinal direction $\alpha$ via $$\alpha = \arctan \left( \frac{\cos\theta_2 \sin\phi_2}{-\cos\phi_2 \cos\theta_2 \sin\theta_1 + \cos\theta_1 \sin\theta_2} \right)$$ This arctan is obviously defined mod $\pi$ only. One has to choose a value between $\pi/2$ and $3\pi/2$ if $\theta_2 < -\theta_1$ i.e. if we go more to the South than to the North. Don't forget that $\phi_2$ above is a shortcut for $\phi_2-\phi_1$.

As a check, note that for $\phi_2=0$, we obtain $\alpha=0$ (North) or $\alpha=\pi$ (South). For $\theta_1=\theta_2=0$ (both on equator), we obtain $\alpha =\arctan(1/0) = \pm \pi/2$ (East, West). As is usual with $\arctan$, the $\pm \pi$ ambiguity is unnecessary. The more accurate condition is $$e^{i\alpha} = r_+ [(-\cos\phi_2 \cos\theta_2 \sin\theta_1 + \cos\theta_1 \sin\theta_2) + i (\cos\theta_2 \sin\phi_2)], \qquad r_+\in R^+$$

-
Comes out surprisingly clean actually. I was expecting something worse. :-) –  Jyrki Lahtonen Jun 14 '11 at 16:08

I don't think that there is a very clean formula. However, here's how you can calculate it. I assume that you are familiar with the spherical coordinates. These are much like the (longitude, latitude) coordinates except that latitude is replaced with another angle parameter that is defined to be zero at North pole, 90 degrees on the equator, and 180 degrees at the South pole.

1) Convert the points, say $A$ and $B$ longitude and latitude into $xyz$-coordinates using the spherical-to-cartesian transformation. You may assume that the units are chosen in such a way that the Earth has radius $=1$. The center of Earth $O$ is at the origin, so you now know the vectors $\vec{OA}$ and $\vec{OB}$.

2) Compute the cross product $\vec{n}=\vec{OA}\times\vec{OB}$. This vector is perpendicular to both vectors $\vec{OA}$ and $\vec{OB}$, so it gives the direction of the normal of the plane $T$ that contains all the three points $O,A,B$. Remark: if $A$ and $B$ are not antipodal points, there is only one such plane.

3) Assume that you are starting from the point $A$. The shortest path from $A$ to $B$ along the surface follows a great circle. This great circle is the intersection of Earth's surface and the plane $T$. We want to calculate a tangent vector $\vec{t}$ of this great circle at the point $A$. Here's how we do it. The great circle follows the surface of the Earth. The vector $\vec{OA}$ is normal to the surface at $A$. Therefore we know that $\vec{t}\perp\vec{OA}$. The great circle also lies on the plane $T$. Therefore we know that $\vec{t}\perp\vec{n}$. From this we can deduce that $\vec{t}$ must be parallel to (i.e. a scalar multiple of) the cross produc $\vec{r}=\vec{n}\times\vec{OA}$. Normalize this to have unit length, so set $\vec{t}=\vec{r}/|r|.$

4) We know the direction we want to go now. We want to convert it to a compass bearing. To that end we need two more unit vectors based at the starting point $A$. One should point at North and the other at East. You get this by partially differentiating the $xyz$-map of the spherical coordinates. This was probably done in a (vector) calculus course. If you don't know how to do this, ask! Let these vectors be $\vec{u}_N$ and $\vec{u}_E$.

5) Using the dot product you can then compute the angle between the vectors $\vec{t}$ and $\vec{u}_N$. If this angle is $\alpha, 0\le\alpha\le180$, then we know that the compass bearing must be either $\alpha$ or $360-\alpha$. Which is it? Well, check whether it has a positive dot product with the East pointing vector or not. If that dot product is positive, then the bearing is in the range $[0,180]$, if negative, then you should be heading into a bearing in the range $[180,360]$.

6) Caveat: We have only computed the tangent direction of the great circle. This means that you could be heading around the globe the wrong way. Shouldn't be too difficult to decide which way to go, but be prepared to switch the bearing by 180 degrees. Hopefully somebody else can say more about this.

-
If you are 'fluent' in Finnish please check out the numerical example I gave to my calculus students. The task is to compute the compass bearing heading from my hometown (22 East, 60 North) to Anchorage, Alaska (151 West, 60 North). Complete with a Mathematica picture of the plane T: users.utu.fi/lahtonen/Analyysi2011Kevat/Isoympyra.pdf –  Jyrki Lahtonen Jun 14 '11 at 16:17
Correcting the URL. –  Jyrki Lahtonen Apr 5 '13 at 16:57