Determine compass direction from one lat/lon to the other

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

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.

-
Where do you want the cut-off point to be for choosing between cardinal and intermediate directions? –  Adam P Jul 9 '10 at 4:48

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

return(bearings[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).

-
Not to be negative, but this isn't even remotely working. –  Cat Man Do Jul 10 '10 at 17:23
Check answer posted below for PHP solution –  soulkphp Apr 10 '11 at 6:47
The latitude and longitude inputs must be converted from degrees to radians before inputing them into these calculations. –  pwnsauce Oct 27 '13 at 22:13

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.

-
<sup>-1</sup> –  erickson Jul 9 '10 at 5:39
Note this fails the "50 miles East and one mile North is simply East and not Northeast" requirement. It also assumes you're talking about a flat plane, rather than a sphere (but that might be OK if the distances are small). –  Dean Harding Jul 9 '10 at 5:42
Ah, I didn't see that; will change the code. –  Lèse majesté Jul 9 '10 at 5:51
@erickson: unfortunately, that doesn't work, or at least not in code blocks. –  Lèse majesté Jul 9 '10 at 6:02
You probably want atan2 since -1/-1 is the opposite direction to 1/1, and this won't give anything like the correct angle except near the equator - at 60 north, north east would be about two degrees east for each degree north. –  Pete Kirkham Jul 9 '10 at 8:42

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.

-

Here's a post with functions for bearing and compass direction in PHP: http://www.thetechnote.com/2011/02/3-useful-geolocation-functions-in-php/

-
{"found_math": false, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0}