MySQL Forums
Forum List  »  GIS

Re: WGS84 Data & distance queries
Posted by: Pierre Herman
Date: September 01, 2007 02:33PM

Here is a function in PHP that returns the distance between two points on earth based on the WGS84 datum. This is nothing new. It is called the Vincenty's folmula and the theory can be found at http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf

This is maybe not what you're looking for but it can be easily translated in any other language and is accurate to less than a meter. You just need to send the latitude and longitude in decimal degrees.

function distInvVincenty($lat1, $lon1, $lat2, $lon2) {
/* WGS84 stuff */
$a = 6378137;
$b = 6356752.3142;
$f = 1/298.257223563;
/* end of WGS84 stuff */
$L = deg2rad($lon2-$lon1);
$U1 = atan((1-$f) * tan(deg2rad($lat1)));
$U2 = atan((1-$f) * tan(deg2rad($lat2)));
$sinU1 = sin($U1);
$cosU1 = cos($U1);
$sinU2 = sin($U2);
$cosU2 = cos($U2);

$lambda = $L;
$lambdaP = 2*pi();
$iterLimit = 20;
while ((abs($lambda-$lambdaP) > pow(10, -12)) && ($iterLimit-- > 0)) {
$sinLambda = sin($lambda);
$cosLambda = cos($lambda);
$sinSigma = sqrt(($cosU2*$sinLambda) * ($cosU2*$sinLambda) + ($cosU1*$sinU2-$sinU1*$cosU2*$cosLambda) * ($cosU1*$sinU2-$sinU1*$cosU2*$cosLambda));
if ($sinSigma==0) {return 0;}
$cosSigma = $sinU1*$sinU2 + $cosU1*$cosU2*$cosLambda;
$sigma = atan2($sinSigma, $cosSigma);
$sinAlpha = $cosU1 * $cosU2 * $sinLambda / $sinSigma;
$cosSqAlpha = 1 - $sinAlpha*$sinAlpha;
$cos2SigmaM = $cosSigma - 2*$sinU1*$sinU2/$cosSqAlpha;
if (is_nan($cos2SigmaM)) {$cos2SigmaM = 0;}
$C = $f/16*$cosSqAlpha*(4+$f*(4-3*$cosSqAlpha));
$lambdaP = $lambda;
$lambda = $L + (1-$C) * $f * $sinAlpha *($sigma + $C*$sinSigma*($cos2SigmaM+$C*$cosSigma*(-1+2*$cos2SigmaM*$cos2SigmaM)));
}
if ($iterLimit==0) {return NaN;} // formula failed to converge

$uSq = $cosSqAlpha * ($a*$a - $b*$b) / ($b*$b);
$A = 1 + $uSq/16384*(4096+$uSq*(-768+$uSq*(320-175*$uSq)));
$B = $uSq/1024 * (256+$uSq*(-128+$uSq*(74-47*$uSq)));
$deltaSigma = $B*$sinSigma*($cos2SigmaM+$B/4*($cosSigma*(-1+2*$cos2SigmaM*$cos2SigmaM)- $B/6*$cos2SigmaM*(-3+4*$sinSigma*$sinSigma)*(-3+4*$cos2SigmaM*$cos2SigmaM)));

$s = $b*$A*($sigma-$deltaSigma);
return "Distance: ".$s;
}

Options: ReplyQuote


Subject
Views
Written By
Posted
8807
April 25, 2006 07:26PM
6376
December 19, 2006 03:42AM
5576
April 18, 2007 01:07AM
Re: WGS84 Data & distance queries
8823
September 01, 2007 02:33PM


Sorry, you can't reply to this topic. It has been closed.

Content reproduced on this site is the property of the respective copyright holders. It is not reviewed in advance by Oracle and does not necessarily represent the opinion of Oracle or any other party.