Then There Was Geo

[alexander]

So i am seeing more and more code online and in phone and web apps that do various geo location things, and sometimes its a bit hard to find or figure out things with geo location, so i am going to post a little code here in hopes that it might help some people with their various "geo needs"

 

       function destination($lat1, $lon1, $brng, $dist) 
{
	// WGS-84 ellipsoid params
	$semiMajorAxis = 6378137; 
	$semiMinorAxis = 6356752.3142;  
	$inverseFlattening=0.003352811; // to cut on redundant calculations, its (1/f) 1/298.257223563

	$alpha1 = deg2rad($brng);
	$sinAlpha1 = sin($alpha1);
	$cosAlpha1 = cos($alpha1);
  
	$tanU1 = (1-$inverseFlattening) * tan(deg2rad($lat1));
	$cosU1 = 1 / sqrt((1 + pow($tanU1,2)));
	$sinU1 = $tanU1*$cosU1;
	$sigma1 = atan2($tanU1, $cosAlpha1);
	$sinAlpha = $cosU1 * $sinAlpha1;
	$cosSqAlpha = 1 - pow($sinAlpha,2);
	$uSq = $cosSqAlpha * (pow($semiMajorAxis,2) - pow($semiMinorAxis,2)) / pow($semiMinorAxis,2);
	$A = 1 + $uSq/16384*(4096+$uSq*(-768+$uSq*(320-175*$uSq)));
	$B = $uSq/1024 * (256+$uSq*(-128+$uSq*(74-47*$uSq)));
  
	$sigma = $dist / ($semiMinorAxis*$A);
	$sigmaP = 2*pi();
	while (abs($sigma-$sigmaP) > 0.000000000001) {
		$cos2SigmaM = cos(2*$sigma1 + $sigma);
		$sinSigma = sin($sigma);
		$cosSigma = cos($sigma);
		$deltaSigma = $B*$sinSigma*($cos2SigmaM+$B/4*($cosSigma*(-1+2*pow($cos2SigmaM,2))-$B/6*$cos2SigmaM*(-3+4*pow($sinSigma,2))*(-3+4*pow($cos2SigmaM,2))));
		$sigmaP = $sigma;
		$sigma = $dist / ($semiMinorAxis*$A) + $deltaSigma;  
	}

	$lat2 = atan2($sinU1*$cosSigma + $cosU1*$sinSigma*$cosAlpha1, (1-$inverseFlattening)*sqrt(pow($sinAlpha,2) + pow($sinU1*$sinSigma - $cosU1*$cosSigma*$cosAlpha1,2)));
	$lambda = atan2($sinSigma*$sinAlpha1, $cosU1*$cosSigma - $sinU1*$sinSigma*$cosAlpha1);
	$C = $inverseFlattening/16*$cosSqAlpha*(4+$inverseFlattening*(4-3*$cosSqAlpha));
	$L = $lambda - (1-$C) * $inverseFlattening * $sinAlpha * ($sigma + $C*$sinSigma*($cos2SigmaM+$C*$cosSigma*(-1+2*pow($cos2SigmaM,2))));

	return array("lat"=>rad2deg($lat2), "long"=>$lon1+rad2deg($L));
}


function distance($lat1, $lon1, $lat2, $lon2, $fast=true) 
{
	// WGS-84 ellipsoid params
	$semiMajorAxis=6378137;
	$semiMinorAxis=6356752.3142;
	$inverseFlattening=0.003352811; // to cut on redundant calculations, its (1/f) 1/298.257223563
	
	if($fast) {
		// Haversine formula, spherical earth, fast
		$lat1 = deg2rad($lat1);
		$lon1 = deg2rad($lon1);
		$lat2 = deg2rad($lat2);
		$lon2 = deg2rad($lon2);
		$dlat = ($lat2-$lat1);
		$dlon = ($lon2-$lon1);
		$a = pow(sin($dlat/2),2)+cos($lat1)*cos($lat2)*pow((sin($dlon/2)),2);
		$a = 2*atan2(sqrt($a), sqrt(1-$a));
		
		return round(($semiMajorAxis*$a),2);
	} else {
		// Vincenti formula, ellipsoid earth model, pretty precise
		$L = deg2rad($lon2-$lon1);
		$U1 = atan((1-$inverseFlattening) * tan(deg2rad($lat1)));
		$U2 = atan((1-$inverseFlattening) * tan(deg2rad($lat2)));
		$sinU1 = sin($U1);
		$cosU1 = cos($U1);
		$sinU2 = sin($U2);
		$cosU2 = cos($U2);
	  
		$lambda = $L; 
		$iterLimit = 100;
		do {
			$sinLambda = sin($lambda);
			$cosLambda = cos($lambda);
			$sinSigma = sqrt(pow(($cosU2*$sinLambda),2) + pow(($cosU1*$sinU2-$sinU1*$cosU2*$cosLambda),2));
			if($sinSigma==0) return 0;  // co-incident points
			$cosSigma = $sinU1*$sinU2 + $cosU1*$cosU2*$cosLambda;
			$sigma = atan2($sinSigma, $cosSigma);
			$sinAlpha = $cosU1 * $cosU2 * $sinLambda / $sinSigma;
			$cosSqAlpha = 1 - $sinAlpha*$sinAlpha;
			$cos2SigmaM = $cosSigma - 2*$sinU1*$sinU2/$cosSqAlpha;
			if (is_nan($cos2SigmaM)) $cos2SigmaM = 0;  // equatorial line: cosSqAlpha=0 (§6)
			$C = $inverseFlattening/16*$cosSqAlpha*(4+$inverseFlattening*(4-3*$cosSqAlpha));
			$lambdaP = $lambda;
			$lambda = $L + (1-$C) * $inverseFlattening * $sinAlpha * ($sigma + $C*$sinSigma*($cos2SigmaM+$C*$cosSigma*(-1+2*pow($cos2SigmaM,2))));
		}
		while (abs($lambda-$lambdaP) > 0.000000000001 && --$iterLimit>0);

		if ($iterLimit==0) return false;  // formula failed to converge

		$uSq = $cosSqAlpha * (pow($semiMajorAxis,2) - pow($semiMinorAxis,2)) / pow($semiMinorAxis,2);
		$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*pow($cos2SigmaM,2))-$B/6*$cos2SigmaM*(-3+4*pow($sinSigma,2))*(-3+4*pow($cos2SigmaM,2))));
		  
		return round($semiMinorAxis*$A*($sigma-$deltaSigma), 3);
	}
}


function toMiles($m)  { return $m/1609.344; }

function toMeters($mi)  { return $mi*1609.344; }

 

So you can use these for various thing, for example you can use these functions to help you search your database for localized things. For example you can use the destination function to help define a square to use as a coordinate match against the database. You could call the destination thing twice with a north west and south east bearings, define your square and generate queries based on that. Actually for that i would suggest using some quick angular offset math (to speed it up (though its more imprecise))

 

       // $lat, $lon are the center point of the search
       // $offset = offset in miles;
       $ang_offset_lat = $offset/69.05;
       $ang_offset_lon = $offset/abs(cos(deg2rad($lat))*69.05);
       $query = "SELECT * FROM `table` WHERE (`latitude` BETWEEN ".($lat-$ang_offset_lat)." AND ".($lat+$ang_offset_lat)." AND `longitude` BETWEEN ".($lon-$ang_offset_lon)." AND " . ($lon+$ang_offset_lon).")";

 

This makes it pretty easy to search for "blank" in a quickly expanding search minimizing the amounts of results returned while making sure that results are always returned. For example, if you were to map all coffee shops in the country and let people search them, by location, you would be able to return a ton of results, maybe a hundred in a 5 mile search in NYC, and there would be thousands of results in the city, but really nobody is ever going to scroll through 2-3000 results, so you save on data, and on database load, which can be (and often times is) beneficial. Where as in the mid-west you could search for a coffee shop and if you expand say doubling your radius (so 4 times the area each time) it would take you two or three queries to find a hand full of coffee shops out there. So you can always return results, and minimize the amount when densely packed together.

 

Distance is always neat to play with too, always nice to return it with the results of the previous example ;) . Just doesn't work for elevation, but that would be a lot more math and data, which is beyond the scope/usability of these much more simple functions...

 

Anyways, look, use, play, hopefully you can find something cool to do with these :)

 

If you have any questions regarding how to use/do something with geo, i will try to help you out. But for now, peace :)