Selecting points from a database by latitude/longitude 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(φ₁)⋅sin(φ₂) + cos(φ₁)⋅cos(φ₂)⋅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.

Node.js JavaScript version

I’ve now (2016+) moved from Apache/PHP to Node.js/JavaScript as my preferred platform.

A full solution written in Node.js / JavaScript might be as follows.

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
/*  Selection of points within specified radius of given lat/lon - Node.js (c) Chris Veness 2020  */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

import mysql from 'mysql2/promise.js';

const R = 6371e3; // earth's mean radius in metres
const sin = Math.sin, cos=Math.cos, acos = Math.acos;
const π = Math.PI;

// query selection parameters: latitude, longitude & radius of bounding circle
const lat = Number(process.argv[2]);    // or e.g. req.query.lat (degrees)
const lon = Number(process.argv[3]);    // or e.g. req.query.lon (degrees)
const radius = Number(process.argv[4]); // or e.g. req.query.radius; (metres)

// set up database connection
const Db = await mysql.createConnection({
    host:              process.env.DB_HOST,
    database:          process.env.DB_DATABASE,
    user:              process.env.DB_USER,
    password:          process.env.DB_PASSWORD,
    namedPlaceholders: true,
});

// query points within first-cut bounding box (Lat & Lon should be indexed for fast query)
const sql = `
    Select Id, Postcode, Lat, Lon
    From MyTable
    Where Lat Between :minLat And :maxLat
      And Lon Between :minLon And :maxLon`;
const params = {
    minLat: lat - radius/R*180/π,
    maxLat: lat + radius/R*180/π,
    minLon: lon - radius/R*180/π / cos(lat*π/180),
    maxLon: lon + radius/R*180/π / cos(lat*π/180),
};
const [ pointsBoundingBox ] = await Db.execute(sql, params);
Db.end(); // close connection

// add in distance d = acos( sinφ₁⋅sinφ₂ + cosφ₁⋅cosφ₂⋅cosΔλ ) ⋅ R
pointsBoundingBox.forEach(p => { p.d = acos(sin(p.Lat*π/180)*sin(lat*π/180) +
    cos(p.Lat*π/180)*cos(lat*π/180)*cos(p.Lon*π/180-lon*π/180)) * R })

// filter for points with distance from bounding circle centre less than radius, and sort
const pointsWithinCircle = pointsBoundingBox.filter(p => p.d < radius).sort((a, b) => a.d - b.d);

console.log(pointsWithinCircle); // or e.g. res.render('points', { points: pointsWithinCircle });

I have moved the distance calculation & radius filtering from a database task to a JavaScript task; I've found it easier to scale / load-balance a web server than a database server, and now prefer to keep database queries as simple as possible.

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.

Movable Type Scripts

Selecting points within a bounding circle

First-cut bounding box

Full solution

Node.js JavaScript version