To select points from a database within a certain distance from a given latitude/longitude point, within a selection circle, you can use the spherical law of cosines formula, which gives the great-circle distance between two points on the earth’s surface. You can use a SQL query which selects points from a database table within a given distance of the centre of a bounding circle.
The cosine law is: d = acos( sin(?1).sin(?2) + cos(?1).cos(?2).cos(??) ).R
where: ? = latitude, ? = longitude
R = radius of earth
d = distance between the points
Written out as it would be coded, this is
d = acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon2-lon1))*R
In this example, I will use PHP, though other web scripting languages such as ASP or ColdFusion will work similarly, as will db procedural languages such as PL/SQL or VBA. The SQL should be common to any DBMS, though the execution methods vary between systems – here I will use MySQL+PDO. I shall assume the table includes fields ID, Postcode, Lat, and Lon (I use double(9,6) for lat/lon, stored as degrees).
If you have parameters $lat and $lon for the centre of the bounding circle (converted from degrees to radians), and $rad for the radius of the bounding circle, and use $R for the radius of the earth (6371km), then the SQL query which selects the points within the circle is:
Code: Select all
Select ID, Postcode, Lat, Lon,
acos(sin($lat)*sin(radians(Lat)) + cos($lat)*cos(radians(Lat))*cos(radians(Lon)-$lon))*$R As D
From MyTable
Where acos(sin($lat)*sin(radians(Lat)) + cos($lat)*cos(radians(Lat))*cos(radians(Lon)-$lon))*$R < $radWe can make it significantly more efficient by doing an initial ‘first cut’ (how much more efficient depends on the number of points and the size of the bounding circle). Dividing the radius of the bounding circle by the radius of the earth ($rad/$R) gives us the angular distance between the centre and the circumference of the bounding circle, in radians, so we just convert this to degrees and extend the centre point into a bounding box:
Code: Select all
// first-cut bounding box (in degrees)
$maxLat = $lat + rad2deg($rad/$R);
$minLat = $lat - rad2deg($rad/$R);
// compensate for degrees longitude getting smaller with increasing latitude
$maxLon = $lon + rad2deg($rad/$R/cos(deg2rad($lat)));
$minLon = $lon - rad2deg($rad/$R/cos(deg2rad($lat)));
$sql = "Select ID, Postcode, Lat, Lon
From MyTable
Where Lat>$minLat And Lat<$maxLat
And Lon>$minLon And Lon<$maxLon";
For some situations, this will suffice – for example, selecting points to display on a square map.
If we want only those points within a given radius, we can use a sub-query to combine the simpler ‘first cut’ with the accurate cosine law:
Code: Select all
Select ID, Postcode, Lat, Lon,
acos(sin($lat)*sin(radians(Lat)) + cos($lat)*cos(radians(Lat))*cos(radians(Lon)-$lon))*$R As D
From (
Select ID, Postcode, Lat, Lon
From MyTable
Where Lat>$minLat And Lat<$maxLat
And Lon>$minLon And Lon<$maxLon
) As FirstCut
Where acos(sin($lat)*sin(radians(Lat)) + cos($lat)*cos(radians(Lat))*cos(radians(Lon)-$lon))*$R < $rad We can put this all together to get a list of points within a bounding circle. In this example, it is a PHP web page which takes URL parameters lat, lon, rad and lists the points as an html page. The example here is simplifed and doesn’t include any error checking etc.
Code: Select all
<?php
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Selection of points within specified radius of given lat/lon (c) Chris Veness 2008-2009 */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
require 'inc/dbparams.inc.php'; // defines $dsn, $username, $password
$lat = $_GET['lat']; // latitude of centre of bounding circle in degrees
$lon = $_GET['lon']; // longitude of centre of bounding circle in degrees
$rad = $_GET['rad']; // radius of bounding circle in kilometers
$R = 6371; // earth's radius, km
// first-cut bounding box (in degrees)
$maxLat = $lat + rad2deg($rad/$R);
$minLat = $lat - rad2deg($rad/$R);
// compensate for degrees longitude getting smaller with increasing latitude
$maxLon = $lon + rad2deg($rad/$R/cos(deg2rad($lat)));
$minLon = $lon - rad2deg($rad/$R/cos(deg2rad($lat)));
// convert origin of filter circle to radians
$lat = deg2rad($lat);
$lon = deg2rad($lon);
$db = new PDO($dsn, $username, $password);
$sql = "
Select ID, Postcode, Lat, Lon,
acos(sin($lat)*sin(radians(Lat)) + cos($lat)*cos(radians(Lat))*cos(radians(Lon)-$lon))*$R As D
From (
Select ID, Postcode, Lat, Lon
From MyTable
Where Lat>$minLat And Lat<$maxLat
And Lon>$minLon And Lon<$maxLon
) As FirstCut
Where acos(sin($lat)*sin(radians(Lat)) + cos($lat)*cos(radians(Lat))*cos(radians(Lon)-$lon))*$R < $rad
Order by D";
$points = $db->query($sql, PDO::FETCH_OBJ);
?>
<html>
<body>
<table>
<? foreach ($points as $point): ?>
<tr>
<td><?= $point->Postcode ?></td>
<td><?= number_format($point->D,1) ?></td>
<td><?= number_format($point->Lat,3) ?></td>
<td><?= number_format($point->Lon,3) ?></td>
</tr>
<? endforeach ?>
</table>
</body>
</html>If you prefer to work in miles, set the earth’s radius $R to 3959.
Courtesy of movable-type.co.uk