Movable Type Home Page

Movable Type Scripts


Selecting points within a bounding circle

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 (in km), 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 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.

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-2014  */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

    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);
    // 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,
                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>
<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>
</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.


Creative Commons License 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-2014 Chris Veness