Go Premium for a chance to win a PS4. Enter to Win

x
?
Solved

PHP function to calculate bearing between latitudes and longitudes

Posted on 2004-08-15
16
Medium Priority
?
3,824 Views
Last Modified: 2012-06-21
I need to be able to calculate the bearing between two latitude and longitude pairs, in PHP. They are decimal degrees, and since we are in Australia, our points we are working with are negative latitude, positive longitude. eg. -29.54546 Latitude and 149.68756 longitude.

I have tried this:
$bearingradians = atan2(asin($long1-$long2)*cos($lat2),
cos($lat1)*sin($lat2) - sin($lat1)*cos($lat2)*cos($long1-$long2));
$bearingdegrees = rad2deg($bearing);

This seems to work, but on closer inspection, it gets east and west mixed up! ie: When I know that a point should be west of a certain point, it comes up with say, 88degrees, which is obviously east! I am putting the lats and longs around the right way, so thats not the problem. I am trying to make a function which will take the four numbers and spit out the bearing between them, like this:

function bearing($lat1, $long1, $lat2, $long2){

$bearingradians = atan2(asin($long1-$long2)*cos($lat2),
cos($lat1)*sin($lat2) - sin($lat1)*cos($lat2)*cos($long1-$long2));
$bearingdegrees = rad2deg($bearing);

return $bearingdegrees
}

Can anyone see what im doing wrong mathematically? I need this working asap...
0
Comment
Question by:Dr_Snapid
  • 10
  • 6
16 Comments
 
LVL 2

Expert Comment

by:catknows
ID: 11812247
My first impression is that the unit is wrong:


1. function bearing($lat1, $long1, $lat2, $long2){
2. $bearingradians = atan2(asin($long1-$long2)*cos($lat2),
3. cos($lat1)*sin($lat2) - sin($lat1)*cos($lat2)*cos($long1-$long2));
4. $bearingdegrees = rad2deg($bearing);

5. return $bearingdegrees
6. }

On line 3, the first term and second term are not of the same unit. You cannot do subtraction between the two term.

If you can tell me how you derive this formular, you can help you right away.

0
 
LVL 2

Expert Comment

by:catknows
ID: 11813334
Continuing above:

I don't know how you derive your formula, but if I use formula #6 on of the below link,

http://mathworld.wolfram.com/TrigonometricAdditionFormulas.html

This is what I get:

function b2($y1, $x1, $y2, $x2)
{

  $bearingradians = atan2( $y1*$x2 - $y2*$x1, $x1*$x2 + $y1*$y2 );

  $bearingdegrees = rad2deg($bearingradians);

  return $bearingdegrees;
}
0
 
LVL 2

Expert Comment

by:catknows
ID: 11813748
Continuing above:

Just a note: The function I wrote takes values that are the coordinates of the longitude and latitude pairs.  The bearing angle is the angle on the plane that is tangent to the surface of the earth, assuming the two point are close enough.  Computing the angle bearing relative to the center of the earth would involve a little more.

Cheers.

Long
0
Industry Leaders: We Want Your Opinion!

We value your feedback.

Take our survey and automatically be enter to win anyone of the following:
Yeti Cooler, Amazon eGift Card, and Movie eGift Card!

 
LVL 1

Author Comment

by:Dr_Snapid
ID: 11815935
Cat knows, heres the code I tried, it didnt work.


function getbearing($y1, $x1, $y2, $x2)
{

  $bearingradians = atan2( $y1*$x2 - $y2*$x1, $x1*$x2 + $y1*$y2 );

  $bearingdegrees = rad2deg($bearingradians);

  return $bearingdegrees;
}


#test the function
// Grafton NSW
$lat1 = -29.6844;
$long1 = 152.9325;
//Coffs Harbour NSW, about 170 degrees from grafton
$lat2 = -30.31319;
$long2 = 153.13753;

echo "The bearing from Grafton to Coffs Harbour is ".getbearing($long1,$lat1,$long2,$lat2);
-------------------------------------------------------------------------------------------------------------
This returns -0.212214774169, which is wrong.... Should be about 170 or so. Can you see anything wrong?
0
 
LVL 2

Expert Comment

by:catknows
ID: 11816515
Oh, I got it.

By the way, if you can tell me from what source you derive your formula that would confirm my finding, too.  Thanks.

My equation compute the angle formed by the vectors having the coordinates of x1,y1, and x2,y2 that is why it does not come out right.

But hang on. Will get back to you.
0
 
LVL 1

Author Comment

by:Dr_Snapid
ID: 11816726
It came from Dr. Math, and now I cant find the URL, will post as soon as i find it.
0
 
LVL 2

Expert Comment

by:catknows
ID: 11816812
Ok,
Try this...

function getbearing($y1, $x1, $y2, $x2)
{

  $bearingradians = atan2( tan($y2-$y1), tan($x2-$x1)  );

  $bearingdegrees = rad2deg($bearingradians);

  return $bearingdegrees;
}

This is the angle originated from the destination relative to the source.
0
 
LVL 1

Author Comment

by:Dr_Snapid
ID: 11816937
Jeez we're close now, worked perfectly for grafton to Coffs Harbour, but then was -108.95732909deg for Grafton to Armidale, which should be 233.2deg.

So you can check:
Grafton -29.6844,152.9325
Coffs Harbour -30.31319,153.13753 (should be 164deg from Grafton, we got 164)
Armidale -30.51194,151.6675 (Should be 233deg from Grafton, we got -109)
Ballina -28.87398,153.56237 (Should be 34deg from Grafton, we got 34)
Tenterfield -29.05639, 152.01806 (should be 309deg, we got -60)

Looks like westerly directions are out.... if I add 360 degrees to them to get a positive value it is still out by a fair way...
You are good at this, what have we missed?
0
 
LVL 1

Author Comment

by:Dr_Snapid
ID: 11816961
If it helps, I chose those towns for being on opposite sides of grafton, Ballina is NE, Coffs is SE, Armidale is SW and Tenterfield is NW.
0
 
LVL 2

Expert Comment

by:catknows
ID: 11817102
Adding 360 is the correct way to get the number rightly orientated.  But since the numbers are still out, I can't think of any otherway.  I compute this purely on mathematically point of view, that is, assuming the earth is perfectly round, and angle are measured as though you are standing to the origin, which is Grafton, of the x-y coordinate.

Where do you get the measurements from?

The reason I am asking is that the formula I used is based on the tangent plane that touch the earth the origin where the angle is located.  If your measure is on the flat map, the angle might not come out accurate.  If you get your data out of a table, can you tell me where the source is?  If I can understand what the data you have really meant what they meant.
0
 
LVL 2

Expert Comment

by:catknows
ID: 11817182
I just add a line to the function to take care of the 360 degree wrap around:

>>
function getbearing($y1, $x1, $y2, $x2)
{

  $bearingradians = atan2( tan($y2-$y1), tan($x2-$x1)  );

  $bearingdegrees = rad2deg($bearingradians);
  if( $bearingdegrees < 0 )
    $bearingdegrees += 360;

  return $bearingdegrees;
}
>>

The result:

Coffs Harbour:
x2-x1=-0.62879
y2-y1=0.20502999999999
The bearing from Grafton to Coffs Harbour: 164.04275769504

Armidale:
x2-x1=-0.82754
y2-y1=-1.265
The bearing from Grafton to Armidale: 251.04267091014 [233]

Ballina:
x2-x1=0.81042
y2-y1=0.62986999999998
The bearing from Grafton to Ballina: 34.734391266122

Tenderfield:
x2-x1=0.62801
y2-y1=-0.91444000000001
The bearing from Grafton to Tenderfield: 299.21742521454 [309]

You can see that the larger the change in the latitude in the negative direction (more south), the worse it appears in the difference between your data and the formula results.

I think this is a matter of practical values, and theoretical values.

If you want to dig more into this, just let me know.

Cheers.
Got to go. Too hungry now.

0
 
LVL 1

Author Comment

by:Dr_Snapid
ID: 11817261
Hmmm I just had lunch. Check out the Chicken Fillet Burgers ant Country Fresh... Mmmmm!

So do you think there is a way to correct the function? It just seems to me there must be another thing weve missed... we are so close...

Someone has done it... look at this:

http://www.zeuswireless.com/support/lat-long.php

It works well. Just wish I could see his code...
0
 
LVL 2

Expert Comment

by:catknows
ID: 11821775
Ok, try this. The latest and greatest.

Cheers. By the way, no chicken for me. Thanks for the suggestion.


<?
function bearing( $lat1_d, $lon1_d, $lat2_d, $lon2_d )
{

   $lat1 = deg2rad($lat1_d);
   $lon1 = deg2rad($lon1_d);
   $lat2 = deg2rad($lat2_d);
   $lon2 = deg2rad($lon2_d);

   $L    = $lon2 - $lon1;

   $cosD = sin($lat1)*sin($lat2) + cos($lat1)*cos($lat2)*cos($L);
   $D    = acos($cosD);
   $cosC = (sin($lat2) - $cosD*sin($lat1))/(sin($D)*cos($lat1));
   
    $C = 180.0*acos($cosC)/pi();

    if( sin($L) < 0.0 )
        $C = 360.0 - $C;

    return $C;

}
?>
0
 
LVL 2

Accepted Solution

by:
catknows earned 2000 total points
ID: 11824828
Continue...

Your formula need to be fixed just a little bit then it will come out just fine, too.

function bearing( $lat1_d, $lon1_d, $lat2_d, $lon2_d )
{
   $lat1 = deg2rad($lat1_d);
   $long1 = deg2rad($lon1_d);
   $lat2 = deg2rad($lat2_d);
   $long2 = deg2rad($lon2_d);

$bearingradians = atan2(asin($long2-$long1)*cos($lat2),cos($lat1)*sin($lat2) - sin($lat1)*cos($lat2)*cos($long2-$long1));
$bearingdegrees = rad2deg($bearingradians);
$bearingdegrees = $bearingdegrees < 0? 360 + $bearingdegrees : $bearingdegrees;


return $bearingdegrees;
};

Cheers!
0
 
LVL 1

Author Comment

by:Dr_Snapid
ID: 11826556
Legend! I knew you could do it!
0
 
LVL 2

Expert Comment

by:catknows
ID: 11859224
Thank you! I am very glad that that helps.
0

Featured Post

Free Tool: Port Scanner

Check which ports are open to the outside world. Helps make sure that your firewall rules are working as intended.

One of a set of tools we are providing to everyone as a way of saying thank you for being a part of the community.

Question has a verified solution.

If you are experiencing a similar issue, please ask a related question

If you’re thinking to yourself “That description sounds a lot like two people doing the work that one could accomplish,” you’re not alone.
Make the most of your online learning experience.
An introduction to basic programming syntax in Java by creating a simple program. Viewers can follow the tutorial as they create their first class in Java. Definitions and explanations about each element are given to help prepare viewers for future …
Loops Section Overview

926 members asked questions and received personalized solutions in the past 7 days.

Join the community of 500,000 technology professionals and ask your questions.

Join & Ask a Question