Determine compass direction from one lat/lon to the other
Asked Answered
K

3

4

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.

Kulturkampf answered 9/7, 2010 at 4:44 Comment(1)
Where do you want the cut-off point to be for choosing between cardinal and intermediate directions?Bastardy
E
21

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

Ehr answered 9/7, 2010 at 4:56 Comment(1)
The latitude and longitude inputs must be converted from degrees to radians before inputing them into these calculations.Prelature
B
2

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.

Belly answered 9/7, 2010 at 4:51 Comment(3)
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).Ehr
@erickson: unfortunately, that doesn't work, or at least not in code blocks.Hilaria
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.Empyrean
S
1

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.

Stendhal answered 9/7, 2010 at 4:51 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.