Determine compass direction from one lat/lon to the other

Question

Does anyone have an algorithm to determine the direction from one lat/lon to another (pseudo-code):

CalculateHeading( lat1, lon1, lat2, long2 ) returns string heading

Where heading is e.g. NW, SW, E, etc.

Basically, I have two points on a map and I want to get a general idea of the direction taking into account that 50 miles East and one mile North is simply East and not Northeast.

Answer 1

Do you remember your trig functions? I.e. SOHCAHTOA:

• SOH: Sin(θ) = Opposite over Hypotenuse
• CAH: Cos(θ) = Adjacent over Hypotenuse
• TOA: Tan(θ) = Opposite over Adjacent

In pseudo-code:

function getDir(lat1, long1, lat2, long2) {
    margin = π/90; // 2 degree tolerance for cardinal directions
    o = lat1 - lat2;
    a = long1 - long2;
    angle = atan2(o,a);

    if (angle > -margin && angle < margin):
            return "E";
    elseif (angle > π/2 - margin && angle < π/2 + margin):
            return "N";
    elseif (angle > π - margin && angle < -π + margin):
            return "W";
    elseif (angle > -π/2 - margin && angle < -π/2 + margin):
            return "S";
    }
    if (angle > 0 && angle < π/2) {
        return "NE";
    } elseif (angle > π/2 && angle < π) {
        return "NW";
    } elseif (angle > -π/2 && angle < 0) {
        return "SE";
    } else {
        return "SW";
    }
}

Edit 1: As Pete and Dean pointed out, this does not take into account the curvature of the Earth. For more accurate calculations for points away from the equator, you'll need to use spherical triangle formulas, as used in Dean's answer.

Edit 2: Another correction; as Pete noted, arctan() does not give the correct angles, as -1/-1 and 1/1 are the same (as are -1/1 and 1/-1). arctan2(y, x) is a two argument variation of arctan() that is designed to compensate for this. arctan() has a range of (-π, π], positive for y >= 0 and negative for y < 0.

Answer 2

Convert to a numeric angle and use the result to look up the text. For example, -22.5..+22.5 = N. +22.5..67.5 = NE, 67.5..112.5 = E, etc. Of course, that's assuming you're using only N, NE, E, SE, S, SW, W, NW -- if you decide (for example) to go with the old "32 points of the compass", each text string obviously represents a smaller range.

Answer 3

This site has the basic algorithm:

// in javascript, not hard to translate...
var y = Math.sin(dLon) * Math.cos(lat2);
var x = Math.cos(lat1)*Math.sin(lat2) -
        Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
var brng = Math.atan2(y, x).toDeg();

UPDATED: See here for complete algorith http://stackoverflow.com/questions/1971585/mapping-math-and-javascript

That'll give you a number between 0 and 360 then it's just a matter of having a simple lookup:

var bearings = ["NE", "E", "SE", "S", "SW", "W", "NW", "N"];

var index = brng - 22.5;
if (index < 0)
    index += 360;
index = parseInt(index / 45);

alert(bearnings[index]);

It's important to note that your bearing actually changes as you move around the earth. The algorithm above shows you initial bearing, but if you're traveling a long distance, your bearing to be significantly different when you reach the destination (if you're only traveling a short distance [< a few hundred kms] then it probably won't change enough to be a concern).