Pixels from the Edge

Okay this is complicated stuff. Working on putting together a new site for a client, I have a table of locations with longitude and latitude fields for every row (I actually had a tab delimited file, but I sucked it into the db). For the functionality on the web front end a user needs to search by distance, e.g. find all locations within 10 miles of a zip code. Um okay, so how the hell do I do this?

Well first of I need to get the latitude and longitude for a zip code. For that I very quickly found the Google Geocode API and threw together the following:

<?php
$result=file_get_contents("http://maps.google.com/maps/geo?48220&sensor=false&gl=us&output=xml&key=yourgooglekey);
$data=simplexml_load_string($result);
$coord=$data->Response->Placemark->Point->coordinates;
$coord=explode(",",$coord);
$longitude=$coord[1];
$latitude=$coord[0];
?>

It took like 10 minutes, and that included signing up for the API key and reading the documentation, shit this geocoding stuff is easy, I should be done by lunch…

Yeah okay so I have two sets of coordinates, great but errr now what? So I find the Haversine formula, but alas not really being any good at math I don’t understand it, but it doesn’t really matter because I know how to convert it to do the math in PHP. Cool so I got the Haversine function in PHP now but how do I query it against the data in my database? I’d read some stuff about geographical libraries for databases and it seemed super complicated, like the GIS and Spatial extensions for MySQL. So I figure I’ll use a two step approach and do a simple query to pull out a square block of locations and then apply the math. i.e. if I am searching for within a 10 mile radius of a location I’ll do a search on -10 miles and +10 miles, but as the distance is circular I would be getting results in the four corners of that square block that are greater than the 10 miles, and I would then apply the Haversine function and be left only with those results that where truly within the 10 mile radius.

By sheer luck I found this post. It was insane, I didn’t realize you could apply the Haversine formula directly in SQL, I ran that query in my db and bam it worked instantly (well actually there’s a typo syntax error ‘< =’ -blank space between the less than and equals character - in the where clause). Obviously I had to switch out a couple of things so it played nice with my database, but it was instant success, I tried out some different latitude and longitude coordinates and it worked every time:

Wow so this was it, everything I needed. Oh wait, except that I needed to display the distance in miles. But the guy who blogged this was some kind of superman and he talks about this on another post. Really I only needed the rad() and distance() functions (because distance is dependent on rad), but I decided that they’ll probably all come in useful one day, so I took them wrapped them in a class called GeoLoc. The class doesn’t need to be instantiated into an object so all methods can be called as follows

<?php
GeoLoc::distance($lat1, $lon1, $lat2, $lon2, $units);
?>

Then decided I could abstract in the SQL statement by passing in all necessary variables, and I added it to the class:

<?php
GeoLoc::prepareGetByDistanceSQL($lat,$long,$distance_in_miles,$table_name,$fields);
?>

I then understood that the measuring unit the query is using is the same as in the distance() method so I abstracted it out into its own method. At which point I noticed that the numbers didn’t match exactly – for miles the SQL had 3963.191 whereas the distance function only had 3963.1. I decided that the SQL query had the more granular number and when I reviewed its original blog post it was across the board for all three measuring units so I updated all the numbers. I then figured I might as well abstract in the Google Geocode zip code lookup.

This is what I was left with (functions I’m not using removed); I give all credit to joeldg:

<?php
/*
Geographical Location functions
ALL CREDIT GOES TO http://blog.peoplesdns.com/
Thanks to these two blog posts:
http://blog.peoplesdns.com/archives/24
http://blog.peoplesdns.com/archives/34
This API:
http://code.google.com/apis/maps/documentation/geocoding/
*/
class GeoLoc{

const GOOGLE_API_KEY="yourkeygoeshere";
const GOOGLE_API_URL="http://maps.google.com/maps/geo";

/*
queries the google api for zip code
returns long and lat coordinates
*/
function getZipcodeCoordinates($zip_code){
  $result=file_get_contents(self::GOOGLE_API_URL."?q=".$zip_code."&sensor=false&gl=us&output=xml&key=".self::GOOGLE_API_KEY);
  $data=simplexml_load_string($result);
  $coord=$data->Response->Placemark->Point->coordinates;
  $coord=explode(",",$coord);
  $r->longitude=$coord[0];
  $r->latitude=$coord[1];
  return $r;
}
function rad($v){
  return ($v*M_PI/180);
}
/*
distance between two points using sherical law of cosines
cos c = cos a cos b + sin a sin b cos C
*/
function distance($lat1, $lon1, $lat2, $lon2, $units = 'miles'){
  $lat1 = self::rad($lat1); $lon1 = self::rad($lon1);
  $lat2 = self::rad($lat2); $lon2 = self::rad($lon2);
  $r=self::convertMeasuringUnit($units);
  return acos(sin($lat1)*sin($lat2) + cos($lat1)*cos($lat2)*cos($lon2-$lon1)) * $r;
}
/*
gets the measuring unit
*/
function convertMeasuringUnit($units){
  switch ($units){
    case "miles": $r = 3963.191; break;;
    case "nmiles": $r = 3441.596; break;;
    case "kilo": $r = 6378.137; break;;
  }
  return $r;
}
/*
returns sql to query a table for rows by distance (miles,nmiles,kilo)
table must contain fields named latitude and longitude
*/
function prepareGetByDistanceSQL($long,$lat,$distance,$table_name,$fields,$units="miles"){
  $cUnit=self::convertMeasuringUnit($units);
  $sql="
    SELECT ".implode($fields,",").",
    acos(SIN( PI()* ".$lat." /180 )*SIN( PI()*latitude/180 ))+(cos(PI()* ".$lat." /180)*COS( PI()*latitude/180) *COS(PI()*longitude/180-PI()* ".$long." /180))* ".$cUnit." AS distance
    FROM ".$table_name."
    WHERE active=1
      AND ".$cUnit." * ACOS( (SIN(PI()* ".$lat." /180)*SIN(PI() * latitude/180)) +
    (COS(PI()* ".$lat." /180)*cos(PI()*latitude/180)*COS(PI() * longitude/180-PI()* ".$long." /180))
    ) <= ".$distance."
    ORDER BY ".$cUnit." * ACOS(
    (SIN(PI()* ".$lat." /180)*SIN(PI()*latitude/180)) +
    (COS(PI()* ".$lat." /180)*Cos(PI()*latitude/180)*COS(PI() * longitude/180-PI()* ".$long." /180))
   )
  ";
  return $sql;
}
}
?>

Tags: Geocoding, Google API, PHP

This entry was posted on Thursday, April 23rd, 2009 at 3:01 pm. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.