My latitude/longitude page presents a range of geodesy formulæ, including the distance between two points. This formula can be used for selecting points from a database.

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 (in radians) R = radius of earth d = distance between the points (in same units as R) |

Using this formula, the SQL query to select points within the bounding circle is:

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 < :rad

Assuming table includes fields Id, Postcode, Lat, and Lon, with Lat/Lon stored as signed degrees (I use double(9,6) for storing degrees, giving sub-metre precision).

### First-cut bounding box

One problem with this is that it works through the entire table, using the cosine law to find the distance of every single point from the given centre point.

We can make it significantly more efficient by using a bounding box to do an initial ‘first cut’ (how much more efficient depends on the number of points and the portion to be included).

The bounding latitudes are obtained by adding/subtracting the radius from the latitude; the bounding longitudes shrink with increasing latitude, so need to be adjusted.

_{bounds}= φ ± d/R

λ

_{bounds}= λ ± (d/R) / cosφ

Actually, as the bounding circle becomes larger, this becomes slightly inaccurate: more strictly, the bounding
latitudes should be^{*}

_{bounds}= λ ± asin(d/R) / cosφ

In this example, I will use PHP, though other web scripting languages such as JavaScript, 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.

If you have parameters $lat and $lon for the centre of the bounding circle, and $rad for the radius of the bounding circle (in km), and use $R for the radius of the earth (6371km), then:

// first-cut bounding box (in degrees) $maxLat = $lat + rad2deg($rad/$R); $minLat = $lat - rad2deg($rad/$R); $maxLon = $lon + rad2deg(asin($rad/$R) / cos(deg2rad($lat))); $minLon = $lon - rad2deg(asin($rad/$R) / cos(deg2rad($lat))); $sql = "Select Id, Postcode, Lat, Lon From MyTable Where Lat Between :minLat And :maxLat And Lon Between :minLon And :maxLon";

Assuming the Lat and Lon columns are indexed, this should be very efficient.

For some situations, this will suffice – for example, selecting points to display on a square map.

If we want to sort the results by distance (or only want those points within a given radius), we can use a sub-query to combine the simpler ‘first cut’ with the accurate cosine law:

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 Between :minLat And :maxLat And Lon Between :minLon And :maxLon ) As FirstCut Where acos(sin(:lat)*sin(radians(Lat)) + cos(:lat)*cos(radians(Lat))*cos(radians(Lon)-:lon)) * :R < :rad Order By D

In this case, we are only applying the more expensive cosine law to points within the bounding box which encompasses the bounding circle.

If the query optimisation is decent, this simpler version should also be efficient (possibly better):

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 Lat Between :minLat And :maxLat And Lon Between :minLon And :maxLon And acos(sin(:lat)*sin(radians(Lat)) + cos(:lat)*cos(radians(Lat))*cos(radians(Lon)-:lon)) * :R < :rad Order By D

### Full solution

We can put this all together to get a list of points within a bounding circle. In this example, it is a single vanilla-PHP web page which takes URL query arguments lat, lon, rad and lists the points as an html page.

In a real application, this would more typically be done within a framework, and based on a geocoded address, but this vanilla-PHP should be recognisable to any developer.

<?php /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Selection of points within specified radius of given lat/lon (c) Chris Veness 2008-2016 */ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ require 'inc/dbparams.inc.php'; // defines $dsn, $username, $password $db = new PDO($dsn, $username, $password); $db->setAttribute(PDO::ATTR_DEFAULT_FETCH_MODE, PDO::FETCH_OBJ); $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 mean radius, km // first-cut bounding box (in degrees) $maxLat = $lat + rad2deg($rad/$R); $minLat = $lat - rad2deg($rad/$R); $maxLon = $lon + rad2deg(asin($rad/$R) / cos(deg2rad($lat))); $minLon = $lon - rad2deg(asin($rad/$R) / cos(deg2rad($lat))); $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 Between :minLat And :maxLat And Lon Between :minLon And :maxLon ) As FirstCut Where acos(sin(:lat)*sin(radians(Lat)) + cos(:lat)*cos(radians(Lat))*cos(radians(Lon)-:lon)) * :R < :rad Order by D"; $params = [ 'lat' => deg2rad($lat), 'lon' => deg2rad($lon), 'minLat' => $minLat, 'minLon' => $minLon, 'maxLat' => $maxLat, 'maxLon' => $maxLon, 'rad' => $rad, 'R' => $R, ]; $points = $db->prepare($sql); $points->execute($params); ?> <html> <table> <? foreach ($points as $point): ?> <tr> <td><?= $point->Postcode ?></td> <td><?= number_format($point->Lat,4) ?></td> <td><?= number_format($point->Lon,4) ?></td> <td><?= number_format($point->D,3) ?></td> </tr> <? endforeach ?> </table> </html>

Note that special treatment would be required for a bounding circle which crosses the 180° meridian – by happy coincidence of geography, not many people need to worry about that!

If you prefer to work in miles, set the earth’s radius $R to 3959.

I offer these scripts for free use and adaptation to balance my debt to the open-source info-verse. You are welcome to re-use these scripts [under an MIT licence, without any warranty express or implied] provided solely that you retain my copyright notice and a link to this page.

If you would like to show your appreciation and support continued development of these scripts, I would most gratefully accept donations.

If you would like any advice or development work done on geodesy-related projects (including databases, server-side-, and client-side-development), I am available for consultancy.

If you have any queries or find any problems, contact me at ku.oc.epyt-elbavom@oeg-stpircs.

*© 2008-2016 Chris Veness*