Great Circle Distance question – Here in this article, we will share some of the most common and frequently asked about PHP problem in programming with detailed answers and code samples. There’s nothing quite so frustrating as being faced with PHP errors and being unable to figure out what is preventing your website from functioning as it should like php and latitude-longitude . If you have an existing PHP-based website or application that is experiencing performance issues, let’s get thinking about Great Circle Distance question.
I am familiar with the formula to calculate the Great Circle Distance between two points.
i.e.
<?php
$theta = $lon1 - $lon2;
$dist = sin(deg2rad($lat1)) * sin(deg2rad($lat2)) + cos(deg2rad($lat1)) * cos(deg2rad($lat2)) * cos(deg2rad($theta));
$dist = acos($dist);
$dist = rad2deg($dist);
//convert degrees to distance depending on units desired
?>
What I need though, is the reverse of this. Given a starting point, a distance, and a simple cardinal NSEW direction, to calculate the position of the destination point. It’s been a long time since I was in a math class. 😉
Solution :
To answer my own question just so it’s here for anyone curious, a PHP class as converted from the C function provided by Chad Birch:
class GreatCircle
{
/*
* Find a point a certain distance and vector away from an initial point
* converted from c function found at: http://sam.ucsd.edu/sio210/propseawater/ppsw_c/gcdist.c
*
* @param int distance in meters
* @param double direction in degrees i.e. 0 = North, 90 = East, etc.
* @param double lon starting longitude
* @param double lat starting latitude
* @return array ('lon' => $lon, 'lat' => $lat)
*/
public static function getPositionByDistance($distance, $direction, $lon, $lat)
{
$metersPerDegree = 111120.00071117;
$degreesPerMeter = 1.0 / $metersPerDegree;
$radiansPerDegree = pi() / 180.0;
$degreesPerRadian = 180.0 / pi();
if ($distance > $metersPerDegree*180)
{
$direction -= 180.0;
if ($direction < 0.0)
{
$direction += 360.0;
}
$distance = $metersPerDegree * 360.0 - $distance;
}
if ($direction > 180.0)
{
$direction -= 360.0;
}
$c = $direction * $radiansPerDegree;
$d = $distance * $degreesPerMeter * $radiansPerDegree;
$L1 = $lat * $radiansPerDegree;
$lon *= $radiansPerDegree;
$coL1 = (90.0 - $lat) * $radiansPerDegree;
$coL2 = self::ahav(self::hav($c) / (self::sec($L1) * self::csc($d)) + self::hav($d - $coL1));
$L2 = (pi() / 2) - $coL2;
$l = $L2 - $L1;
$dLo = (cos($L1) * cos($L2));
if ($dLo != 0.0)
{
$dLo = self::ahav((self::hav($d) - self::hav($l)) / $dLo);
}
if ($c < 0.0)
{
$dLo = -$dLo;
}
$lon += $dLo;
if ($lon < -pi())
{
$lon += 2 * pi();
}
elseif ($lon > pi())
{
$lon -= 2 * pi();
}
$xlat = $L2 * $degreesPerRadian;
$xlon = $lon * $degreesPerRadian;
return array('lon' => $xlon, 'lat' => $xlat);
}
/*
* copy the sign
*/
private static function copysign($x, $y)
{
return ((($y) < 0.0) ? - abs($x) : abs($x));
}
/*
* not greater than 1
*/
private static function ngt1($x)
{
return (abs($x) > 1.0 ? self::copysign(1.0 , $x) : ($x));
}
/*
* haversine
*/
private static function hav($x)
{
return ((1.0 - cos($x)) * 0.5);
}
/*
* arc haversine
*/
private static function ahav($x)
{
return acos(self::ngt1(1.0 - ($x * 2.0)));
}
/*
* secant
*/
private static function sec($x)
{
return (1.0 / cos($x));
}
/*
* cosecant
*/
private static function csc($x)
{
return (1.0 / sin($x));
}
}
Here’s a C implementation that I found, should be fairly straightforward to translate to PHP:
#define KmPerDegree 111.12000071117
#define DegreesPerKm (1.0/KmPerDegree)
#define PI M_PI
#define TwoPI (M_PI+M_PI)
#define HalfPI M_PI_2
#define RadiansPerDegree (PI/180.0)
#define DegreesPerRadian (180.0/PI)
#define copysign(x,y) (((y)<0.0)?-fabs(x):fabs(x))
#define NGT1(x) (fabs(x)>1.0?copysign(1.0,x):(x))
#define ArcCos(x) (fabs(x)>1?quiet_nan():acos(x))
#define hav(x) ((1.0-cos(x))*0.5) /* haversine */
#define ahav(x) (ArcCos(NGT1(1.0-((x)*2.0)))) /* arc haversine */
#define sec(x) (1.0/cos(x)) /* secant */
#define csc(x) (1.0/sin(x)) /* cosecant */
/*
** GreatCirclePos() --
**
** Compute ending position from course and great-circle distance.
**
** Given a starting latitude (decimal), the initial great-circle
** course and a distance along the course track, compute the ending
** position (decimal latitude and longitude).
** This is the inverse function to GreatCircleDist).
*/
void
GreatCirclePos(dist, course, slt, slg, xlt, xlg)
double dist; /* -> great-circle distance (km) */
double course; /* -> initial great-circle course (degrees) */
double slt; /* -> starting decimal latitude (-S) */
double slg; /* -> starting decimal longitude(-W) */
double *xlt; /* <- ending decimal latitude (-S) */
double *xlg; /* <- ending decimal longitude(-W) */
{
double c, d, dLo, L1, L2, coL1, coL2, l;
if (dist > KmPerDegree*180.0) {
course -= 180.0;
if (course < 0.0) course += 360.0;
dist = KmPerDegree*360.0-dist;
}
if (course > 180.0) course -= 360.0;
c = course*RadiansPerDegree;
d = dist*DegreesPerKm*RadiansPerDegree;
L1 = slt*RadiansPerDegree;
slg *= RadiansPerDegree;
coL1 = (90.0-slt)*RadiansPerDegree;
coL2 = ahav(hav(c)/(sec(L1)*csc(d))+hav(d-coL1));
L2 = HalfPI-coL2;
l = L2-L1;
if ((dLo=(cos(L1)*cos(L2))) != 0.0)
dLo = ahav((hav(d)-hav(l))/dLo);
if (c < 0.0) dLo = -dLo;
slg += dLo;
if (slg < -PI)
slg += TwoPI;
else if (slg > PI)
slg -= TwoPI;
*xlt = L2*DegreesPerRadian;
*xlg = slg*DegreesPerRadian;
} /* GreatCirclePos() */
Source: http://sam.ucsd.edu/sio210/propseawater/ppsw_c/gcdist.c
It would be harder to back out the Haversine fomula, then generate your own, I would think.
First the angle generated from the Earth core by traveling a “straight” line on the surface (you think it is straight, but it is curving).
Angle in Radians = Arc Length / radius.
Angle = ArcLen/6371 km
Latitude should be easy, just the “vertical” (north/south) component of your angle.
Lat1 + Cos(bearing) * Angle
Longitude divisions vary by latitude. So that becomes harder. You would use:
Sin(bearing) * Angle (with East defined as negative) to find the angle in longitude direction, but converting back to actual longitude at that latitude would be more difficult.
See the section Destination point given distance and bearing from start point at this website: http://www.movable-type.co.uk/scripts/latlong.html