Global proximity box calculations

Posted on September 7, 2009 in Blog, Social networking, Software Development | 0 comments

By Dwight Irving

Note:  Updated on Dec 6, 2009 – Found a bug in the code.  I’ve updated the code below and have a notation where the bug was.

A while ago, I needed to figure out how to search for database entries that were within a specified distance from a central point.  The code was needed for the search function at CrossroadsAngel.com where you can search for members by location.  I found many pages that had a similar goal of calculating great circle distance between any two points on a sphere, but that’s not quite the same problem.  I found a few web pages that described how to convert linear distance to delta latitude and longitude, but I didn’t find any complete code examples that worked across all border conditions.

The key equation for translating linear distance to latitude and longitude on a sphere is the Haversine formula.  Latitude boundary conditions for this formula occur at the  north and south poles and the equator.  For longitude, the boundary conditions occur at the prime and 180th degree meridians. For distance, the boundary condition is the circumference of the Earth.

I made a simplifying assumption that the earth is a perfect sphere.  To make the database search easier I used a rectangular bounding box rather than a circular one.  The rectangle used encloses the circle and therefore increases the distance from the centerpoint at any angles other than 0, 90, 180 and 270 degrees (0, PI/2, PI, 3PI/2 radians).

I thought it would also be handy to be able to use either kilometers or miles in the function, so I’ve included a flag where miles are true and kilometers are false

In PHP the code is as follows:

class diGeoBox
{

public $LatN;
public  $LatS;
public  $LonE;
public  $LonW;

public function __construct($Lat1, $Lon1, $Distance, $boolMiles)
{

$RadCalcConst=M_PI/180;
$InvRadCalcConst=180/M_PI;

// Convert degrees to radians
$Lat1 = $Lat1 * $RadCalcConst;
$Lon1Rad = $Lon1 * $RadCalcConst;

// Miles or kilometers?

if ($boolMiles)
$DistRad = $Distance * 2.52780586e-4; // 1/3956 radius of the earth in miles
else
$DistRad = $Distance * 1.57070574e-4; // inverse radius of the earth in km

if ($DistRad >= 2*M_PI)
{

$this->LatN = 90;
$this->LatS = -90;
$this->LonE = 180;
$this->LonW = -180;
return;

}

//     North Latitude boundary
$this->LatN = $Lat1 + $DistRad;
// Convert back to degrees
$this->LatN = $this->LatN * $InvRadCalcConst;
// check for wrapping
if ($this->LatN > 90.0)
$this->LatN = 90.0 – fmod($LatN,90.0);

// South Latitude
$this->LatS  = $Lat1 – $DistRad;
// Convert back to degrees
$this->LatS = $this->LatS * $InvRadCalcConst;
// check for wrapping
if ($this->LatS < -90.0)
$this->LatS = -90.0;

// if at a pole, Longititude goes from -180 to 180
if ( abs($Lat1) == M_PI/2)
{

$this->LonE = 180;
$this->LonW = -180;

}
else
{

$tmp=sin($Lat1);
$tmp=$tmp*$tmp;

// calculate East Longitude boundary
$dLon = abs(atan2(sin($DistRad) * cos($Lat1), cos($DistRad)- $tmp));
$this->LonE = fmod( $Lon1Rad + $dLon +M_PI,2*M_PI )-M_PI;
// convert back to degrees
$this->LonE = $this->LonE * $InvRadCalcConst;

//check to see if wrapped around 180 longitude
//**  There was a bug here in the original version where
//**  Lon1 used to be Lon1Rad
if ($this->LonE < $Lon1)
{

// LonW is an equal amount on the other side of $Lon1 as $LonE
$this->LonW = $Lon1-($this->LonE + 180);

}
else
{

// LonW is an equal amount on the other side of $Lon1 as $LonE
$this->LonW = 2*$Lon1 – $this->LonE;
// See if it wrapped
if ($this->LonW < -180)
$this->LonW = $this->LonW + 360;

}

}

}     //__construct()

}         //diGeoBox

Example php code using the diGeoBox function and tests the 180 degree meridian and the polar conditions is below.

$BoundingBox = new diGeoBox( 0, -170, 15000, 1);
echo ‘LonE = ‘.$BoundingBox->LonE. “\n”;
echo ‘LonW = ‘.$BoundingBox->LonW. “\n”;
echo ‘LatN = ‘.$BoundingBox->LatN. “\n”;
echo ‘LatS = ‘.$BoundingBox->LatS. “\n”;

Using this function at CrossroadsAngel.com requires that street address and zip codes be converted into latitude and longitude coordinates (geocoding).  Because CrossroadsAngel.com is a commercial website, and because I need to store converted locations in a local database for performance reasons, I am using the NAC Geographic Products webservice for stored conversion, and Google Maps API for geocoding the “search center” location.  The Google Maps API is accessed using JavaScript from the client browser and the converted coordinates are sent to the server for query execution.