Federico Cargnelutti

I’m working on a project that requires Geo proximity search. Basically, what I’m doing is plotting a radius around a point on a map, which is defined by the distance between two points on the map given their latitudes and longitudes. To achieve this I’m using the Haversine formula (spherical trigonometry). This equation is important in navigation, it gives great-circle distances between two points on a sphere from their longitudes and latitudes. You can see it in action here: Radius From UK Postcode.

This has already been covered in some PHP blogs, however, I found most of the information to be inaccurate and, in some cases, incorrect.

Of course, there are always exceptions. In this post, Kevin uses the Haversine equation to develop a postcode radius facility for his website, like the one on freemaptools.com.

He assumes you have the following values:

$lon = (float) -2.708077;
$lat = (float) 53.754842;
$radius = 20; // in miles

And calculates the latitude (min/max) and longitude (min/max) based on that:

(float)$dpmLAT = 1 / 69.1703234283616;

// Latitude calculation
(float)$usrRLAT = $dpmLAT * $radius;
(float)$latMIN = $lat - $usrRLAT;
(float)$latMAX = $lat + $usrRLAT;

// Longitude calculation
(float)$mpdLON = 69.1703234283616 * cos($userLat * (pi/180));

// degrees per mile longintude
(float)$dpmLON = 1 / $mpdLON;
$usrRLON = $dpmLON * $radius;
$lonMIN = $lon - $usrRLON;
$lonMAX = $lon + $usrRLON;

The variable names are a bit confusing, but you can tell straight away he knows what he’s doing. Unfortunately, the script doesn’t work, $userLat is not defined. If you change $userLat to $lat, it works. Also, you can use the deg2rad() function to convert the $lat number to the radian equivalent.

Here is Kevin’s implementation refactored (Codepad)

// Degrees per mile latitude
$dpm_lat = (float) 1 / 69;

// Latitude calculation
$rlat = (float) $dpm_lat * $radius;
$lat_min = (float) $latitude - $rlat;
$lat_max = (float) $latitude + $rlat;

// Longitude calculation
$mpd_lng = (float) abs(cos(deg2rad($latitude)) * 69);
$dpm_lng = (float) 1 / $mpd_lng;

$rlng = $dpm_lng * $radius;
$lng_min = $longitude - $rlng;
$lng_max = $longitude + $rlng; 

echo 'lng (min/max): ' . $lng_min . '/' . $lng_max;
echo 'lat (min/max): ' . $lat_min . '/' . $lat_max;

Outputs:

lng (min/max): -3.1983251898421/-2.2178288101579
lat (min/max): 53.464986927536/54.044697072464

Cool, it works. Now, for the sake of James Inman, lets refactor the script again and simplify the hole thing. The Haversine equation is very straight forward, so there’s no need to complicate things.

I’ve implemented the solution in PHP, Python and SQL. Use the one that suits you best.

PHP implementation (Codepad)

$longitude = (float) -2.708077;
$latitude = (float) 53.754842;
$radius = 20; // in miles

$lng_min = $longitude - $radius / abs(cos(deg2rad($latitude)) * 69);
$lng_max = $longitude + $radius / abs(cos(deg2rad($latitude)) * 69);
$lat_min = $latitude - ($radius / 69);
$lat_max = $latitude + ($radius / 69);

echo 'lng (min/max): ' . $lng_min . '/' . $lng_max . PHP_EOL;
echo 'lat (min/max): ' . $lat_min . '/' . $lat_max;

It outputs the same result:

lng (min/max): -3.1983251898421/-2.2178288101579
lat (min/max): 53.464986927536/54.044697072464

SQL implementation (Codepad)

set @latitude=53.754842;
set @longitude=-2.708077;
set @radius=20;

set @lng_min = @longitude - @radius/abs(cos(radians(@latitude))*69);
set @lng_max = @longitude + @radius/abs(cos(radians(@latitude))*69);
set @lat_min = @latitude - (@radius/69);
set @lat_max = @latitude + (@radius/69);

SELECT * FROM postcode
WHERE (longitude BETWEEN @lng_min AND @lng_max)
AND (latitude BETWEEN @lat_min and @lat_max);

Python implementation (Codepad)

from __future__ import division
import math

longitude = float(-2.708077)
latitude = float(53.754842)
radius = 20

lng_min = longitude - radius / abs(math.cos(math.radians(latitude)) * 69)
lng_max = longitude + radius / abs(math.cos(math.radians(latitude)) * 69)
lat_min = latitude - (radius / 69)
lat_max = latitude + (radius / 69)

print 'lng (min/max): %f %f' % (lng_min, lng_max)
print 'lat (min/max): %f %f' % (lat_min, lat_max)

That’s it. Happy Geolocating!