Calculating The Bearing And Compass Rose Direction Between Two Latitude / Longitude Coordinates In PHP

Does any one know how to work out the bearing from one set of coordinates to another.

For example, if the calculated bearing was 90 (degrees) the using if statements to see if a coordinate was between two parameters I could return North, North East, East and so on?

Thanks in advance.

Matt

Fortunately, the same article at Movable Type Scripts that fed the original post, also gives us formulas for calculating rhumb-line and great-circle bearings between two known points.

Rhumb Line Bearings Versus Great Circle Bearings

A rhumb line is a line of constant bearing, that crosses all meridians (longitudes) at the same angle. In other words, were we to lay down a map, with Points A and B marked on it, and draw a straight line connecting them, that’s a rhumb line.

But remember, Earth is not flat. Usually, the most efficient (i.e., the shortest) way to travel across it is an arc, especially when the distances to be traveled are long; and when we do so, the bearing of our travel isn’t constant, and we don’t start out traveling directly toward our objective.

Therefore, if we determine the great-circle bearing from Point A to Point B, we’re actually getting an initial bearing of travel toward Point B; that bearing only remains valid for some distance, until we get to some other point (usually, the next meridian), at which point we need to recalculate the bearing to Point B, based on where we are.

Fortunately, the distances we are using are not so great as to cause a significant variance in the bearings between rhumb lines and great circles. In other words, the points are so close, we’re practically traveling in a straight line, anyway. Therefore, we can use either the rhumb line bearing or the great circle bearing, because they should be the same.

For the purposes of this demo, I will provide both great-circle and rhumb line bearing calculations. If you are using very long distances (hundreds or thousands of miles), you will likely find the rhumb line bearings provide a better result than the great-circle bearing; for distances 100 miles and under, either method should produce basically the same results.

Update, July 30, 2009: My extensive thanks to tpk_ for pointing out that I did not properly amend the formula at Movable Type Scripts. Doing so has removed the wild variations I was getting between rhumb line and great circle bearings in previous versions of this code. I have amended this post to reflect the correct code.

The Rhumb Line And Great Circle Formulas

In PHP, one calculates the great-circle bearing, in degrees, from starting Point A, expressed as \$lat1 and \$lon1, to remote Point B, expressed as \$lat2 and \$lon2, thus:

```\$bearing = (rad2deg(atan2(sin(deg2rad(\$lon2) - deg2rad(\$lon1)) * cos(deg2rad(\$lat2)), cos(deg2rad(\$lat1)) * sin(deg2rad(\$lat2)) - sin(deg2rad(\$lat1)) * cos(deg2rad(\$lat2)) * cos(deg2rad(\$lon2) - deg2rad(\$lon1)))) + 360) % 360;
```

Calculating the bearing of a rhumb line is far more complicated, and pretty much requires a function. Here it is in PHP:

```function getRhumbLineBearing(\$lat1, \$lon1, \$lat2, \$lon2) {
//difference in longitudinal coordinates
\$dLon = deg2rad(\$lon2) - deg2rad(\$lon1);

//difference in the phi of latitudinal coordinates
\$dPhi = log(tan(deg2rad(\$lat2) / 2 + pi() / 4) / tan(deg2rad(\$lat1) / 2 + pi() / 4));

//we need to recalculate \$dLon if it is greater than pi
if(abs(\$dLon) > pi()) {
if(\$dLon > 0) {
\$dLon = (2 * pi() - \$dLon) * -1;
}
else {
\$dLon = 2 * pi() + \$dLon;
}
}
//return the angle, normalized
return (rad2deg(atan2(\$dLon, \$dPhi)) + 360) % 360;
}
```

I won’t get into the specifics of these calculations, both because it’s a distraction and in all honesty, it’s been a very, very long time since I learned trigonometry, and being even more brutally honest, I’m fairly certain I didn’t actually learn it. If you are interested in the methodology of these formulas, you can read full narratives at MovableType.co.uk, the Aviation Formulary and Wikipedia.

Getting A Compass Rose Direction

Both the formula and function above will return a numeric value between 0 and 359. If you want to convert that to a more meaningful compass rose direction, you would divide the result by 22.5 (for the 16-point compass rose of ENE, SSW, etc.), then use a switch to assign a string value.

Here’s a PHP function that accepts as its argument a numeric bearing and returns a 16-point compass rose direction:

```function getCompassDirection(\$bearing) {
\$tmp = round(\$bearing / 22.5);
switch(\$tmp) {
case 1:
\$direction = "NNE";
break;
case 2:
\$direction = "NE";
break;
case 3:
\$direction = "ENE";
break;
case 4:
\$direction = "E";
break;
case 5:
\$direction = "ESE";
break;
case 6:
\$direction = "SE";
break;
case 7:
\$direction = "SSE";
break;
case 8:
\$direction = "S";
break;
case 9:
\$direction = "SSW";
break;
case 10:
\$direction = "SW";
break;
case 11:
\$direction = "WSW";
break;
case 12:
\$direction = "W";
break;
case 13:
\$direction = "WNW";
break;
case 14:
\$direction = "NW";
break;
case 15:
\$direction = "NNW";
break;
default:
\$direction = "N";
}
return \$direction;
}
```

Two notes about the preceding function: First, the PHP round() function, when not provided with precision or mode arguments, rounds up fractions .5 or more, and all lesser fractions down. Therefore, this function will tend to cast directions clockwise around the compass rose. (More on this in the section about a 16-point compass rose vs. an 8-point compass rose.)

I use a switch with a default condition of North because North can be either 349-359 degrees, or 0-11 degrees. By setting the default condition to North, I simplify the function. If I didn’t do this, I would need two repetitive switch conditions: one to cover 349 degrees or greater; another to cover 0 to 11 degrees.

Again, elegance should always be the rule. The formulas are shown so complex because I want readers to be able to easily relate the source of the equation to the way it is written in PHP.

A 16-Point Compass Rose: Adequate Accuracy For Most Purposes

We could use an eight-point compass rose, (e.g., NE, SW), but the problem with that becomes accuracy; or, more specifically, the willingness of others to pick nits with your directions.

Let’s examine a compass rose:

Suppose we calculate a bearing of 225 degrees from Point A to Point B. If we divide 225 by 45 (the number of degrees between directions in an eight-point compass rose), we get 5. If we count off an eight-point compass rose’s directions clockwise, starting with NE as 1, then 5 equals SW.

But what if the bearing from Point A to Point B is 250 degrees? An eight-point compass rose would say the bearing is W [round(250 / 45) = 6], but some people would argue that 250 degrees bears a bit too far south to be considered West, and that the bearing should be SW.

That’s the primary reason for using the 16-point compass rose; unless you need absolute nautical precision, it will generally return a fair assessment of direction and prevent such arguments.

Incorporating The New Code

So, if you want to incorporate all we’ve discussed here into the canned code for “Getting All ZIP Codes In A Given Radius From A Known Point / ZIP Code Via PHP And MySQL,” you’d first add the getCompassDirection() and getRhumbLineBearing() functions, above, to your page.

Next, you would locate this block of code:

```//output all matches to screen
echo "<table class=\"bordered\" cellspacing=\"0\">\n";
echo "<tr><th>City</th><th>State</th><th>ZIP Code</th><th>Latitude</th><th>Longitude</th><th>Miles, Point A To B</th></tr>\n";
while(\$row = mysql_fetch_array(\$rs)) {
echo "<tr><td>\$row[city]</td><td>\$row[state]</td><td>\$row[zip_code]</td><td>\$row[latitude]</td><td>\$row[longitude]</td><td>";
echo acos(sin(deg2rad(\$lat1)) * sin(deg2rad(\$row['latitude'])) + cos(deg2rad(\$lat1)) * cos(deg2rad(\$row['latitude'])) * cos(deg2rad(\$row['longitude']) - deg2rad(\$lon1))) * \$r;
echo "</td></tr>\n";
}
```

You would replace it with this block of code:

```//output all matches to screen
echo "<table class=\"bordered\" cellspacing=\"0\">\n";
echo "<tr><th>City</th><th>State</th><th>ZIP Code</th><th>Latitude</th><th>Longitude</th><th>Miles, Point A To B</th><th>Rhumb Line Bearing</th><th>Great Circle Bearing</th></tr>\n";
while(\$row = mysql_fetch_array(\$rs)) {
echo "<tr><td>\$row[city]</td><td>\$row[state]</td><td>\$row[zip_code]</td><td>\$row[latitude]</td><td>\$row[longitude]</td>";
\$truedistance = acos(sin(deg2rad(\$lat1)) * sin(deg2rad(\$row['latitude'])) + cos(deg2rad(\$lat1)) * cos(deg2rad(\$row['latitude'])) * cos(deg2rad(\$row['longitude']) - deg2rad(\$lon1))) * \$r;
echo "<td>\$truedistance</td>";
echo "<td>" . getCompassDirection(getRhumbLineBearing(\$lat1, \$lon1, \$row['latitude'], \$row['longitude'])) . " [" . getRhumbLineBearing(\$lat1, \$lon1, \$row['latitude'], \$row['longitude']) . "&deg;]</td>";
\$bearing = (rad2deg(atan2(sin(deg2rad(\$row['longitude']) - deg2rad(\$lon1)) * cos(deg2rad(\$row['latitude']))), cos(deg2rad(\$lat1)) * sin(deg2rad(\$row['latitude'])) - sin(deg2rad(\$lat1)) * cos(deg2rad(\$row['latitude'])) * cos(deg2rad(\$row['longitude']) - deg2rad(\$lon1))) + 360) % 360;
echo "<td>" . getCompassDirection(\$bearing) . " [\$bearing&deg;]</td>";
echo "</tr>\n";
}
```

And that’s all there is to it. You can see a working demo here:

http://demo.dougv.com/php_zip_code_distance_bearing/

17 Comments

1. tpk_ says:

Excel’s ATAN2() function takes x, y as paramenters, as opposed to y, x of most languages (PHP included). Hence, your great circle calculations are broken. Change the order and you’ll be fine.

Great article otherwise.

1. tpk_: That explains everything. Thanks a million times for the heads-up. My code and post have been corrected to reflect your advice.

2. mustak says:

please tell me how to calculate east bearing from available co ordinates…

3. This is amazingly handy – I converted your code into C# and integrated it into my app. Worked first time, after some considerable time I’d already spent trying to hack around myself (I remember trigonometry even worse than you claim to!).

I think there’s actually a mistake in your rhumb line bearing calculation. You write:

```		if(\$dLon > 0) {
\$dLon = (2 * pi() - \$dLon) * -1;
}
else {
\$dLon = 2 * pi() * \$dLon;
}
```

but the equivalent piece of code on movable type says:

`dLon = dLon>0 ? -(2*Math.PI-dLon) : (2*Math.PI+dLon);`

The first option is equivalent, but in the second you are multiplying by dLon and the code at movable type is adding. I have to admit to not having the faintest idea which is correct, but I *think* it’s the addition one.

1. Thanks for the catch, Chris. I’ve corrected the post, demo and code. (PS: Love your Web site. British idioms amuse me to no end.)

4. Garth Powell says:

Did a “cut and paste” on the bearings formula above (thank you)…took out the “\$bearing” part, left the “=”. In place of “\$lat2” and so on, I put the specific cell for the specific information. Pressed “enter” and my 2003 version of Excel didn’t like or know what to do with “% 360;”. What am I missing here? Do I need to leave the “\$” in all the time? Will do so when I run a number of these as the starting point will be a constant for a given set and I also took out spaces in the formula…mistake on my part? Thanks, Garth

5. Garth Powell says:

Okay, when I get to that point (which will be soon, I’ll (most likely) send something from your wish list…in looking at the books and media, you’re nicely diverse in what interests you…some of them I may want to get to read!

GAP

6. The following mod to the getCompassDirection function will correct the problem “this function will tend to cast directions clockwise around the compass rose”:

`\$tmp = floor((\$bearing + 12.25) / 22.5);`
7. Adam says:

AMAZING! I’ve modified one of your queries to add a calculated bearing field to query results. So I pass the lat/lon of an origin point into the query and all of the locations in the db come back with an extra bearing field (holding the calculated bearing between the origin point and each record’s lat/lon).

Thank you Thank you!

8. Jonathan says:

@Dan. Agreed about tendency to cast clockwise but your fix is going the wrong ‘direction’ ?

On a 8 point compass, I want each of the 8 directions to cover a 360/8=45deg sector but I MAY want the 45deg sector for North, as a good example, starting at -22.5deg to +22.5deg, rather than 0 to 45deg? For a 16 point compass the sectors are 22.5 degs so the offset is -11.25 to +11.25deg.

so

`\$tmp = floor((\$bearing - 11.25) / 22.5);`
9. Peter says:

Awesome thanks

10. Thank you very much Doug, helped me a lot.
Cheers,
Chris

11. You made my day. I spent 10 hours struggling with various equations from different resources trying to figure out how to get the great circle’s angle between two points in order to determine the Qibla direction. Thanks a million for your neat, simple and ready to use code.

This site uses Akismet to reduce spam. Learn how your comment data is processed.

• Check out the Commenting Guidelines before commenting, please!
• Want to share code? Please put it into a GitHub Gist, CodePen or pastebin and link to that in your comment.
• Just have a line or two of markup? Wrap them in an appropriate SyntaxHighlighter Evolved shortcode for your programming language, please!