For a project I’ve been working on, I needed to be able to calculate the great-cirlce distance between two points on the globe given their latitudes and longitudes. There has been much discussion on the topic (because the data is of great utility), and I’ve certainly not been able to add anything to it because, for the most part, it deals with spherical trigonometry that I only sort of remember from some math courses I took as a freshman. I have, however (I think), been able to figure it out and write some PHP that will do the calculation.
(warning, this gets into some pretty nasty math, so you may just want to skip to the PHP at the end if that’s what you’re interested in)
We start with the Haversine formula as presented by R.W. Sinnott in Sky and Telescope (1984). For this example, we’ll calculate the distance between Orlando, Florida and Shanghai, China. The latitude/longitude data used is from the Getty Thesaurus of Geographic Names Online (which is cool as hell, I might add).
L1 = latitude of the first point, 28.5333
G1 = longitude of the first point, -81.3667
L2 = latitude of the second point, 31.1000
G2 = longitude of the second point, 121.3667
R = the Earth’s mean radius, 6371 km
Δlat = |L2 - L1|
Δlong = |G2 - G1|
a = sin2(Δlat / 2) + cos(L1) * cos(L2) * sin2(Δlong / 2)
If a is less than 1 then
c = 2 * asin(√a)
otherwise,
c = 2 * asin(1);
Finally, distance = R * c, or 12960.57 kilometers.
Follow all that? Good. But there’s a problem. The above formula assumes that the Earth is a perfect sphere, which it is not. Earth is ellipsoidal (or, more precisely, oblate spheroidal). Thus its radius does not remain constant but instead varies between 6378km at the equator and 6357km at the poles. By assuming that Earth is a sphere, we’re introducing a worst-case error of about 0.5% into our calculations. Five meters per kilometer may not seem like much (and really, it isn’t, but I’m anal about this sort of stuff), but it can be rather easily reduced so we might as well do it.
What we’re going to do is calculate the approximate radius at the latitudes that we’re dealing with by calculating the local radius of curvature. There are still errors, of course, but they’re reduced significantly.
L = the mean latitude, (L1 + L2) / 2
A = equatorial radius, 6378
B = polar radius, 6357
E = the eccentricity of the ellipsoid, √(1 - B2 / A2)
Now we’ll calculate two radii, R’ for travelling directly north/south between those two latitudes and R” for travelling perpendicular to R’:
R’ = A * (1 - E2) / (1 - E2 * sin2(L))3/2
R” = A / √(1 - E2 * sin2(L))
OK, we’re almost done. Now we’ve got to take a weighted average of the two radii, to account for the actual direction of travel. This part is a little fudged, and could probably be done more precisely, but this method seems to work and isn’t particularly computationally expensive.
R avg = R’ * (Δlat / (Δlat + Δlong)) + R” * (Δlong / (Δlat + Δlong))
Now just plug R avg as the radius of the Earth into the Haversine formula we started with and you’ll be good to go. In our example, the distance between Orlando and Shanghai comes back as 12984.56 kilometers, a difference of about 24 kilometers. Not particularly significant in most applications, but when error is that easy to fix there’s little reason not to.
Now, as promised, the code to do this in PHP. Note that PHP’s trig functions take their arguments in radians, so everything is converted beforehand.
$lat1 = deg2rad(28.5333);
$lat2 = deg2rad(31.1000);
$long1 = deg2rad(-81.3667);
$long2 = deg2rad(121.3667);
$dlat = abs($lat2 - $lat1);
$dlong = abs($long2 - $long1);
$l = ($lat1 + $lat2) / 2;
$a = 6378;
$b = 6357;
$e = sqrt(1 - ($b * $b)/($a * $a));
$r1 = ($a * (1 - ($e * $e))) / pow((1 - ($e * $e) * (sin($l) * sin($l))), 3/2);
$r2 = $a / sqrt(1 - ($e * $e) * (sin($l) * sin($l)));
$ravg = ($r1 * ($dlat / ($dlat + $dlong))) + ($r2 * ($dlong / ($dlat + $dlong)));
$sinlat = sin($dlat / 2);
$sinlon = sin($dlong / 2);
$a = pow($sinlat, 2) + cos($lat1) * cos($lat2) * pow($sinlon, 2);
$c = 2 * asin(min(1, sqrt($a)));
$d = $ravg * $c;
echo "<p>Distance: $d kilometers</p>";
Your code does not work properly for some values. eg for
$lat1 = deg2rad(-23.85);
$lat2 = deg2rad(-24.7166667);
$long1 = deg2rad(151.25);
$long2 = deg2rad(152.1166667);
It returns -32437474355223 kilometers which is obviously wrong.
Left by Ken Allan on August 14th, 2003
You need to make:
$dlat = abs($lat2 - $lat1);
$dlong = abs($long2 - $long1);
Left by Ken Allan on August 14th, 2003
Ah, indeed! Thanks much, Ken. I’m not sure how I overlooked that.
Left by John on August 14th, 2003
Hi John,
I’d like to use your php code above in a college project I’m working on. Could you mail me or post here to say if you would kindly give permission to re-use it? I tried emailing you but the message bounced.
Left by Eoin Dubsky on February 23rd, 2004
What’s the difference between the local radius of curvature and the local radius ?
Left by Fabio Ponti on June 24th, 2004
I’m too old *85″ to figure out the following:
At the 26.1 latitude, what is the distance between the 82 degree and 83 degree longitudes.
I’m trying to guess the distance of a hurricand projection from the shorline.
Left by Ian Naylor on September 11th, 2004
99.9656 Kms.
Left by David Rubin on October 4th, 2004
This doesn’t appear to take into account the equatorial bulge / flattening of the earth problem. From my research, the best I can find to do this is using Vincenty’s theorems, but I have not been able to get a good PHP script for it.
Left by David on October 14th, 2004
Actually, it does take a shot at the bulge/flattening problem, finding the average radius at the latitudes in questions, but it is by no means perfect (especially if the line between your points pass through a wide range of latitudes). It’s just quick and simple, and useful for most things that don’t require real precision.
Vincenty’s theorems do indeed appear to be the best solution, but my math skills (and available free time) aren’t quite up to coding that.
Left by John on October 15th, 2004
Can any help me to find Maximum longitude and latitude difference between two points
separated by a given distance
its very urgent ……
thanx and regards
satya
Left by Satyabir Singh on April 20th, 2005
Here is a PHP version of Vincenty’s theorem, converted from Javascript @ http://www.gpsvisualizer.com/calculators.html
It is more accurate, but the code by John executes in 1/3rd of the time.
function getVincenty($lat1,$lon1,$lat2,$lon2){
if(abs($lat1)>90 || abs($lon1) >180 || abs($lat2) >90 || abs($lon2) > 180){
return NULL;
}
$lat1 = deg2rad($lat1);
$lat2 = deg2rad($lat2);
$lon1 = deg2rad($lon1);
$lon2 = deg2rad($lon2);
$a = 6378137;
$b = 6356752.3142;
$f = 1.0 / 298.257223563;
$L = $lon2 - $lon1;
$U1 = atan((1 - $f) * tan($lat1));
$U2 = atan((1 - $f) * tan($lat2));
$sinU1 = sin($U1);
$sinU2 = sin($U2);
$cosU1 = cos($U1);
$cosU2 = cos($U2);
$lambda = $L;
$lambdaP = 2 * M_PI;
$iterLimit = 20;
while( abs($lambda - $lambdaP) > 0.000000000001 && –$iterLimit > 0){
$sinLambda = sin($lambda);
$cosLambda = cos($lambda);
$sinSigma = sqrt(($cosU2 * $sinLambda) * ($cosU2 * $sinLambda) +
($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda) * ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda));
$cosSigma = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosLambda;
$sigma = atan2($sinSigma,$cosSigma);
$alpha = asin($cosU1 * $cosU2 * $sinLambda / $sinSigma);
$cosSqAlpha = cos($alpha) * cos($alpha);
$cos2SigmaM = $cosSigma - 2*$sinU1*$sinU2 / $cosSqAlpha;
$C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
$lambdaP = $lambda;
$lambda = $L + (1 - $C) * $f * sin($alpha) * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma*(-1 + 2 * $cos2SigmaM * $cos2SigmaM)));
}
if($iterLimit==0){
return acos(1.01);//Return NaN. Couldn’t find where it is defined in PHP
}
$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);
$dist = $s / 1000;
return $dist;
}
Left by proxyissues on August 16th, 2005
There is no special checking for 0 distances.. If the same point is entered in twice you get a divide by zero error.
In the original something like:
if(($dlat + $dlong) == 0){
return 0;
}
Needs to be added to avoid this.
Left by proxyissues on August 17th, 2005
Hello! I need to find out the linear mileage (how many miles further west) between the longitudes of Reno, Nevada and Los Angeles, California.
Reno, NV = 119 49 00
Los Angeles = 118 15 00
I need to take into account the curve of the longitude, obviously, but I’d like to know, how much further west IS Reno in relation to Los Angeles?
Thanks!
Left by Tre Gibbs on November 23rd, 2005
How far in miles is it from Lat. 43.59476 to 44.53361
Left by Don Toth on May 14th, 2006
asdfasdf
Left by Anonymous on February 7th, 2007