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 |

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:

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

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

// 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 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 only those points within a given radius (or want to sort the results by distance), 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.

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.

<?php /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Selection of points within specified radius of given lat/lon (c) Chris Veness 2008-2013 */ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 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 mean 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))); $db = new PDO($dsn, $username, $password); $db->setAttribute(PDO::ATTR_DEFAULT_FETCH_MODE, PDO::FETCH_OBJ); $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 = array( '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> <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>

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 formulæ & scripts for free use and adaptation as my contribution to the open-source info-sphere from which I have received so much. You are welcome to re-use these scripts [under a simple attribution license, without any warranty express or implied] provided solely that you retain my copyright notice and a link to this page.

If you would like any development work done on geodesy-related projects (including databases, server-side-, and client-side-development), please contact me at ku.oc.epyt-elbavom@ofni.

If you would like to show your appreciation, I would most gratefully accept donations.

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

*© 2008-2013 Chris Veness*