Pixels from the Edge

In my post yesterday I talked about calculating geographical distances. During my research I’d come across this http://www.scribd.com/doc/2569355/Geo-Distance-Search-with-MySQL. It was a super interesting read and helped me get a much better understanding of what it is I am even doing. Also it really helps you gets super optimized with the whole subject, but as the locations table I am working with has less than 3000 entries my queries were already coming back super fast in less than 1000th of a second, so it really didn’t need to be optimized more. I didn’t really have time to mess around and investigate it all, and seeing as I still didn’t really fully understand the whole subject too well it was something I’d come back to at another time if need be. However it did point me to geonames.org that has free databases with geographical data to download. This is incredibly useful and whenever I’ve needed this data in the past I’ve had my client buy it (the disk is like $30 so no biggie, but it seems this might even be the free source they are pulling it all from). They have individual country zip codes files available to export, the US file contains exactly what I needed – all the US zip codes with their respective latitude and longitude coordinates. This site has a whole range of useful features on the subject – they appear to even have their own webservice, and a list client libraries in various languages, which pointed me towards a PHP library that appears to be doing a lot of what I need (but also way more than what I need), but by this point I already have everything I needed and had already spent too much time on the subject (not really but I didn’t have as much time allocated to the subject as it really needs – geo-location is a subject unto itself but only a small piece of the website I am building).

So geonames.org offers a free webservice, but I was already using the google API which is also free so I didn’t see any benefit to switch. But they did offer for free all the zip codes with the geographical coordinates so I decided to import this into my database and cut out the 3rd party google Geocode API. It would mean quicker response time and not having to worry about if the service is unavailable or returning errors, I was actually planning to cache the results locally anyway. And perhaps the zip code table might have to be updated on an annual basis to maintain its accuracy (like do zip codes really change that much anyways?), but the web service API is subject to change and it would be less maintainable to monitor it and react to the changes.

So I downloaded the tab delimited US Zip Codes file and imported it into a table in my database called “tbl_Zipcodes_Spatial”, with the following fields (from the included readme):

I created an index on the postal_code field and updated my GeoPos:: prepareGetByDistanceSQL() method to get the zip code from the new table. I actually created a new method so that the class still has the ability to pass in the coordinates to the old method when necessary:

<?php
function prepareGetByDistanceWithZipcodeSQL($zip_code,$distance,$table_name,$fields,$units="miles"){
  $cUnit=self::convertMeasuringUnit($units);
  $sql="
    SELECT ".$table_name.".".implode($fields,",".$table_name.".").",
      tbl_zipcodes_spatial.Latitude AS Zipcode_Latitude,
      tbl_zipcodes_spatial.Longitude AS Zipcode_Longitude,
      ACOS(SIN( PI()* tbl_zipcodes_spatial.Latitude /180 )*SIN( PI()*".$table_name.".Latitude/180 ))+(COS(PI()* tbl_zipcodes_spatial.Latitude /180)*COS( PI()*".$table_name.".Latitude/180) *COS(PI()*".$table_name.".Longitude/180-PI()* tbl_zipcodes_spatial.Longitude /180))* ".$cUnit." AS Distance
    FROM tbl_zipcodes_spatial,".$table_name."
    WHERE
      tbl_zipcodes_spatial.postal_code='".$zip_code."'
      AND ".$cUnit." * ACOS( (SIN(PI()* tbl_zipcodes_spatial.Latitude /180)*SIN(PI() * ".$table_name.".Latitude/180)) +
        (COS(PI()* tbl_zipcodes_spatial.Latitude /180)*COS(PI()*".$table_name.".Latitude/180)*COS(PI() * ".$table_name.".Longitude/180-PI()* tbl_zipcodes_spatial.Longitude /180))
        ) <= ".$distance."
    ORDER BY ".$cUnit." * ACOS(
      (SIN(PI()* tbl_zipcodes_spatial.Latitude /180)*SIN(PI()*".$table_name.".Latitude/180)) +
        (COS(PI()* tbl_zipcodes_spatial.Latitude /180)*COS(PI()*".$table_name.".Latitude/180)*COS(PI() * ".$table_name.".Longitude/180-PI()* tbl_zipcodes_spatial.Longitude /180))
      )
  ";
  return $sql;
}
?>

And it worked! Barely increased the query execution time, and I’ve completely cut out the google API middle man so the overall execution time of the page is greatly increased. Now I feel like I’ve put together something pretty damn robust (the final touch would be to create a stored procedure in the database, and this is something I’ll come back to at a future point). Behold the final script:

<?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/
And the zipcode database found here:
http://www.geonames.org/
*/

class GeoLoc{

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

	/*
	queries the google api for zip code
	returns long and lat coordinates
	*/
	public 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;
	}

	/*
	returns the radius
	*/
	public function rad($v){
		return ($v*M_PI/180);
	}

	/*
	degree,minute,seconds to decimal degrees
	note: this could also be done in a single line using ternary
	*/
	public function dms2deg($D,$M,$S,$dir){
		if(strpos(' WsSs', $dir)>0){
			return(-1 * ($D + ($M + $S/60)/60));
		}else{
			return($D + ($M + $S/60)/60);
		}
	}

	/*
	degree,minute,seconds
	*/
	public function dms($rad) {
		$d = abs($rad * 180 / M_PI);
		$d += 1/7200; // add ½ second for rounding
		$deg = floor($d);
		$min = floor(($d-$deg)*60);
		$sec = floor(($d-$deg-$min/60)*3600);
		// add leading zeros if required
		if ($deg< 100) $deg = '0' + $deg;
		if ($deg< 10) $deg = '0' + $deg;
		if ($min< 10) $min = '0' + $min;
		if ($sec< 10) $sec = '0' + $sec;
		return $deg + '\u00B0' + $min + '\u2032' + $sec + '\u2033';
	}

	/*
	gets the measuring unit
	*/
	public 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;
	}

	/*
	distance between two points using sherical law of cosines
	cos c = cos a cos b + sin a sin b cos C
	*/
	public 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;
	}

	/*
	initial bearing from point 1 to point 2
	*/
	public function bearing($lat1,$lon1, $lat2, $lon2) {
		$y = sin($lon2-$lon1) * cos($lat2);
		$x = cos($lat1)*sin($lat2) - sin($lat1)*cos($lat2)*cos($lon2-$lon1);
		return atan2($y, $x);
	}

	/*
	find out what the midpoint between two points is
	*/
	public function midpoint($lat1, $lon1, $lat2, $lon2) {
		$lat1 = self::rad($lat1); $lon1 = self::rad($lon1);
		$lat2 = self::rad($lat2); $lon2 = self::rad($lon2);
		$dLon = $lon2 - $lon1;
		$Bx = (cos($lat2) * cos($dLon));
		$By = (cos($lat2) * sin($dLon));
		$lat3 = atan2( sin($lat1) + sin($lat2), sqrt( (cos($lat1)+$Bx) * (cos($lat1)+$Bx) + ($By * $By)) );
		$lon3 = $lon1 + atan2($By, cos($lat1) + $Bx);

		if (!$lat3 || !$lon3) return false;
		return array($lat3 * 180 / M_PI, $lon3 * 180 / M_PI);
	}

	/*
	returns sql to query a table for rows by distance (miles,nmiles,kilo)
	table must contain fields named latitude and longitude
	*/
	public 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;
	}

	/*
	returns sql to query a table for rows by distance (miles,nmiles,kilo)
	table must contain fields named latitude and longitude
	requires you have the tbl_zipcodes_spatial table http://download.geonames.org/export/zip/
	*/
	public function prepareGetByDistanceWithZipcodeSQL($zip_code,$distance,$table_name,$fields,$units="miles"){
		$cUnit=self::convertMeasuringUnit($units);
		$sql="
			SELECT ".$table_name.".".implode($fields,",".$table_name.".").",
				tbl_zipcodes_spatial.Latitude AS Zipcode_Latitude,
				tbl_zipcodes_spatial.Longitude AS Zipcode_Longitude,
				ACOS(SIN( PI()* tbl_zipcodes_spatial.Latitude /180 )*SIN( PI()*".$table_name.".Latitude/180 ))+(COS(PI()* tbl_zipcodes_spatial.Latitude /180)*COS( PI()*".$table_name.".Latitude/180) *COS(PI()*".$table_name.".Longitude/180-PI()* tbl_zipcodes_spatial.Longitude /180))* ".$cUnit." AS Distance
			FROM tbl_zipcodes_spatial,".$table_name."
			WHERE
				tbl_zipcodes_spatial.postal_code='".$zip_code."'
				AND ".$cUnit." * ACOS( (SIN(PI()* tbl_zipcodes_spatial.Latitude /180)*SIN(PI() * ".$table_name.".Latitude/180)) +
				(COS(PI()* tbl_zipcodes_spatial.Latitude /180)*COS(PI()*".$table_name.".Latitude/180)*COS(PI() * ".$table_name.".Longitude/180-PI()* tbl_zipcodes_spatial.Longitude /180))
				) <= ".$distance."
			ORDER BY ".$cUnit." * ACOS(
				(SIN(PI()* tbl_zipcodes_spatial.Latitude /180)*SIN(PI()*".$table_name.".Latitude/180)) +
				(COS(PI()* tbl_zipcodes_spatial.Latitude /180)*COS(PI()*".$table_name.".Latitude/180)*COS(PI() * ".$table_name.".Longitude/180-PI()* tbl_zipcodes_spatial.Longitude /180))
			)
		";
		return $sql;
	}
}
?>