Solved

PHP function to calculate bearing between latitudes and longitudes

Posted on 2004-08-15
16
3,310 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
 
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
How your wiki can always stay up-to-date

Quip doubles as a “living” wiki and a project management tool that evolves with your organization. As you finish projects in Quip, the work remains, easily accessible to all team members, new and old.
- Increase transparency
- Onboard new hires faster
- Access from mobile/offline

 
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 500 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

6 Surprising Benefits of Threat Intelligence

All sorts of threat intelligence is available on the web. Intelligence you can learn from, and use to anticipate and prepare for future attacks.

Join & Write a Comment

Suggested Solutions

Title # Comments Views Activity
scoresIncreasing challenge 10 57
sumHeights2  challenge 7 78
Replace a tag with sed 2 43
mapAB Challlenge 35 89
A short article about problems I had with the new location API and permissions in Marshmallow
Whether you’re a college noob or a soon-to-be pro, these tips are sure to help you in your journey to becoming a programming ninja and stand out from the crowd.
In this fourth video of the Xpdf series, we discuss and demonstrate the PDFinfo utility, which retrieves the contents of a PDF's Info Dictionary, as well as some other information, including the page count. We show how to isolate the page count in a…
In this fifth video of the Xpdf series, we discuss and demonstrate the PDFdetach utility, which is able to list and, more importantly, extract attachments that are embedded in PDF files. It does this via a command line interface, making it suitable …

760 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

Need Help in Real-Time?

Connect with top rated Experts

24 Experts available now in Live!

Get 1:1 Help Now