This article describes how to efficiently query a database for places that are within a certain distance from a point given in spherical coordinates (latitude and longitude). The method computes bounding coordinates that can be used for a database index scan – just like we would use minimum bounding rectangles to speed up queries in Cartesian space. It is accurate for short distances as well as for long distances, and works for any geographic location, even if it is near a pole or the 180th meridian.
The shortest distance (the geodesic) between two given points P_{1}=(lat_{1}, lon_{1}) and P_{2}=(lat_{2}, lon_{2}) on the surface of a sphere with radius R is the great circle distance. It can be calculated using the formula:
For example, the distance between the Statue of Liberty at (40.6892°, -74.0444°)=(0.7102 rad, -1.2923 rad) and the Eiffel Tower at (48.8583°, 2.2945°)=(0.8527 rad, 0.0400 rad) – assuming a spherical approximation^{a} of the figure of the Earth with radius R=6371 km – is:
Suppose we want to find places within a distance d=1000 km from M=(lat, lon)=(1.3963 rad, -0.6981 rad) in a database. Given that we have a table named Places with columns Lat and Lon that hold the coordinates in radians (the trigonometric functions of SQL expect radians), then we could use this SQL query:
SELECT * FROM Places WHERE acos(sin(1.3963) * sin(Lat) + cos(1.3963) * cos(Lat) * cos(Lon - (-0.6981))) * 6371 <= 1000;
The problem with this query is that our database system will not be able to utilize indices on the Lat and Lon columns because of the complexity of the formula in the WHERE clause.
In order to utilize indices on the Lat and Lon columns we can use a method similar to the bounding rectangle method in Cartesian space. Let’s call the circle that is formed by all points which have distance d from M the query circle. The following subsections explain how to compute bounding coordinates (lat_{min}, lon_{min}) and (lat_{max}, lon_{max}) which are opposite corners of a bounding rectangle (on the sphere) that completely contains the query circle. Points outside this rectangle are therefore definitely not within distance d from M, so in order to find points that are within the query circle, we need to consider the points within the bounding rectangle only. Those points (called candidate places) can be found quickly using an index (see section 3.5 for SQL queries).
Let’s define r = d/R = (1000 km)/(6371 km) = 0.1570 as the angular radius of the query circle.
Moving on a circle with radius R=6371 km from a point A to a point B such that there is an angle r=0.1570 between A and B means covering a distance of d=1000 km. Meridians are on such great circles with radius R. Hence we can move along a meridian, i.e. keep the longitude fixed, and simply substract/add r from/to lat in order to obtain the minimum/maximum latitude of all points within the query circle with center M=(lat, lon)=(1.3963 rad, -0.6981 rad):
Note that special care must be taken if a pole is within the query circle. See section 3.4 for details.
One approach seen on several web pages for computing the minimum and maximum longitudes is keeping the latitude fixed and changing the longitude, i.e. moving along a circle of latitude. This section will show that this approach gives inaccurate results.
Moving along a circle of latitude means moving along a small circle. A circle of latitude at latitude lat=1.3963 rad has the radius R_{s} = R · cos(lat) = 1106 km, so d=1000 km now corresponds to an angular radius of r_{s} = d/R_{s} = d/(R · cos(lat)) = 0.9039. Hence, covering d=1000 km on a circle of latitude gets you to longitude lon_{S} = lon ± d/(R · cos(lat)) = -0.6981 rad ± 0.9039 rad. But this is not the minimum/maximum longitude you can get to by moving d=1000 km in any direction! This is because we moved the distance on a small circle, but small circle distances are greater than great circle distances. Although M=(lat, lon)=(1.3963 rad, -0.6981 rad) and P_{S}=(lat, lon_{S})=(1.3963 rad, -1.6020 rad) have a distance of d=1000 km on the short circle, we could take a shortcut from M to P_{S} on a great circle. According to equation 1, that shortcut has a length of 967 km. So we could move another 33 km and would still be within the query circle. In particular, we could reach even smaller respectively greater longitudes. So using lon ± d/(R · cos(lat)) as bounding values for longitude would miss some locations that are actually within distance d. For example, the point P_{3}=(1.4618 rad, -1.6021 rad) is outside the computed bounding “box”, although its distance to M is only 872 km. That is an error of more than 12 %!
Moving along a circle of latitude in order to find the minimum and maximum longitude does not work at all as you can see in figure 1: The points on the query circle having the minimum/maximum longitude, T_{1} and T_{2}, are not on the same circle of latitude as M but closer to the pole. The formulae for the coordinates of these points can be found in a good math handbook like [1]. They are:
Note that special care must be taken if the 180th meridian is within the query circle. See section 3.4 for details.
If lat_{max}, as defined in section 3.1, is greater than π/2, then the North Pole is within the query circle, which means that parts of all meridians are within the query circle as well. Thus, in that case the bounding coordinates are (lat_{min}, -π) and (π/2, π). If lat_{min} is smaller than -π/2, then the South Pole is within the query circle and the bounding coordinates are (-π/2, -π) and (lat_{max}, π).
If one of lon_{min/max}, as defined in section 3.3, is outside the range of valid longitude values [-π, π], the 180th meridian is within the query circle. In that case one could use (lat_{min}, -π) and (lat_{max}, π) as the bounding coordinates, i.e. something with a shape like a belt around the sphere. But that would include many candidate places that are not actually within the distance. Alternatively one could use two bounding “boxes” and it is sufficient for a candidate place to be in one of them. If lon_{min} < -π, the two sets of bounding coordinates are (lat_{min}, lon_{min} + 2π), (lat_{max}, π) and (lat_{min}, -π), (lat_{max}, lon_{max}). If lon_{max} > π, the two sets of bounding coordinates are (lat_{min}, lon_{min}), (lat_{max}, π) and (lat_{min}, -π), (lat_{max}, lon_{max} - 2π).
Now that we have computed the bounding coordinates (lat_{min}, lon_{min})=(1.2393 rad, -1.8184 rad) and (lat_{max}, lon_{max})=(1.5532 rad, 0.4221 rad), we can use them in an SQL query. In the following, it is assumed that latitude and longitude – both in radians – are stored in two separate columns with an index that supports range queries (e.g. a B-tree) on each of them. If your database management system supports a point data type and a spatial index (e.g. an R-tree) on such columns, you can also use that.
There are different ways to combine the filter using the bounding coordinates with the formula for the great circle distance in an SQL statement:
SELECT * FROM Places WHERE (Lat => 1.2393 AND Lat <= 1.5532) AND (Lon >= -1.8184 AND Lon <= 0.4221) AND acos(sin(1.3963) * sin(Lat) + cos(1.3963) * cos(Lat) * cos(Lon - (-0.6981))) <= 0.1570;Most query optimizers are smart enough to perform an index scan to quickly find places satisfying
(Lat >= 1.2393 AND Lat <= 1.5532) AND (Lon >= -1.8184 AND Lon <= 0.4221)
and evaluate the formula for the great circle distance for each remaining candidate result only.SELECT * FROM Places WHERE (Lat >= 1.2393 AND Lat <= 1.5532) AND (Lon >= -1.8184 AND Lon <= 0.4221) HAVING acos(sin(1.3963) * sin(Lat) + cos(1.3963) * cos(Lat) * cos(Lon - (-0.6981))) <= 0.1570;
SELECT * FROM ( SELECT * FROM Places WHERE (Lat >= 1.2393 AND Lat <= 1.5532) AND (Lon >= -1.8184 AND Lon <= 0.4221) ) WHERE acos(sin(1.3963) * sin(Lat) + cos(1.3963) * cos(Lat) * cos(Lon - (-0.6981))) <= 0.1570;
Below you will find a ready-to-use Java GeoLocation class. The boundingCoordinates() method computes the bounding coordinates using the method explained in section 3. Section 4.2 shows an example usage of the GeoLocation class.
Download GeoLocation.java/** * <p>Represents a point on the surface of a sphere. (The Earth is almost * spherical.)</p> * * <p>To create an instance, call one of the static methods fromDegrees() or * fromRadians().</p> * * <p>This code was originally published at * <a href="http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates#Java"> * http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates#Java</a>.</p> * * @author Jan Philip Matuschek * @version 22 September 2010 */ public class GeoLocation { private double radLat; // latitude in radians private double radLon; // longitude in radians
private double degLat; // latitude in degrees private double degLon; // longitude in degrees private static final double MIN_LAT = Math.toRadians(-90d); // -PI/2 private static final double MAX_LAT = Math.toRadians(90d); // PI/2 private static final double MIN_LON = Math.toRadians(-180d); // -PI private static final double MAX_LON = Math.toRadians(180d); // PI private GeoLocation () { } /** * @param latitude the latitude, in degrees. * @param longitude the longitude, in degrees. */ public static GeoLocation fromDegrees(double latitude, double longitude) { GeoLocation result = new GeoLocation(); result.radLat = Math.toRadians(latitude); result.radLon = Math.toRadians(longitude); result.degLat = latitude; result.degLon = longitude; result.checkBounds(); return result; } /** * @param latitude the latitude, in radians. * @param longitude the longitude, in radians. */ public static GeoLocation fromRadians(double latitude, double longitude) { GeoLocation result = new GeoLocation(); result.radLat = latitude; result.radLon = longitude; result.degLat = Math.toDegrees(latitude); result.degLon = Math.toDegrees(longitude); result.checkBounds(); return result; } private void checkBounds() { if (radLat < MIN_LAT || radLat > MAX_LAT || radLon < MIN_LON || radLon > MAX_LON) throw new IllegalArgumentException(); } /** * @return the latitude, in degrees. */ public double getLatitudeInDegrees() { return degLat; } /** * @return the longitude, in degrees. */ public double getLongitudeInDegrees() { return degLon; } /** * @return the latitude, in radians. */ public double getLatitudeInRadians() { return radLat; } /** * @return the longitude, in radians. */ public double getLongitudeInRadians() { return radLon; } @Override public String toString() { return "(" + degLat + "\u00B0, " + degLon + "\u00B0) = (" + radLat + " rad, " + radLon + " rad)"; }
/** * Computes the great circle distance between this GeoLocation instance * and the location argument. * @param radius the radius of the sphere, e.g. the average radius for a * spherical approximation of the figure of the Earth is approximately * 6371.01 kilometers. * @return the distance, measured in the same unit as the radius * argument. */ public double distanceTo(GeoLocation location, double radius) { return Math.acos(Math.sin(radLat) * Math.sin(location.radLat) + Math.cos(radLat) * Math.cos(location.radLat) * Math.cos(radLon - location.radLon)) * radius; } /** * <p>Computes the bounding coordinates of all points on the surface * of a sphere that have a great circle distance to the point represented * by this GeoLocation instance that is less or equal to the distance * argument.</p> * <p>For more information about the formulae used in this method visit * <a href="http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates"> * http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates</a>.</p> * @param distance the distance from the point represented by this * GeoLocation instance. Must me measured in the same unit as the radius * argument. * @param radius the radius of the sphere, e.g. the average radius for a * spherical approximation of the figure of the Earth is approximately * 6371.01 kilometers. * @return an array of two GeoLocation objects such that:<ul> * <li>The latitude of any point within the specified distance is greater * or equal to the latitude of the first array element and smaller or * equal to the latitude of the second array element.</li> * <li>If the longitude of the first array element is smaller or equal to * the longitude of the second element, then * the longitude of any point within the specified distance is greater * or equal to the longitude of the first array element and smaller or * equal to the longitude of the second array element.</li> * <li>If the longitude of the first array element is greater than the * longitude of the second element (this is the case if the 180th * meridian is within the distance), then * the longitude of any point within the specified distance is greater * or equal to the longitude of the first array element * <strong>or</strong> smaller or equal to the longitude of the second * array element.</li> * </ul> */ public GeoLocation[] boundingCoordinates(double distance, double radius) { if (radius < 0d || distance < 0d) throw new IllegalArgumentException(); // angular distance in radians on a great circle double radDist = distance / radius; double minLat = radLat - radDist; double maxLat = radLat + radDist; double minLon, maxLon; if (minLat > MIN_LAT && maxLat < MAX_LAT) { double deltaLon = Math.asin(Math.sin(radDist) / Math.cos(radLat)); minLon = radLon - deltaLon; if (minLon < MIN_LON) minLon += 2d * Math.PI; maxLon = radLon + deltaLon; if (maxLon > MAX_LON) maxLon -= 2d * Math.PI; } else { // a pole is within the distance minLat = Math.max(minLat, MIN_LAT); maxLat = Math.min(maxLat, MAX_LAT); minLon = MIN_LON; maxLon = MAX_LON; } return new GeoLocation[]{fromRadians(minLat, minLon), fromRadians(maxLat, maxLon)}; } }Download GeoLocation.java
Here is some demo code that shows how to use the GeoLocation class from section 4.1 to query a database for places within a certain distance.
Download GeoLocationDemo.java/** * <p>See * <a href="http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates#Java"> * http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates#Java</a> * for the GeoLocation class referenced from this code.</p> * * @author Jan Philip Matuschek * @version 26 May 2010 */ public class GeoLocationDemo { /** * @param radius radius of the sphere. * @param location center of the query circle. * @param distance radius of the query circle. * @param connection an SQL connection. * @return places within the specified distance from location. */ public static java.sql.ResultSet findPlacesWithinDistance( double radius, GeoLocation location, double distance, java.sql.Connection connection) throws java.sql.SQLException { GeoLocation[] boundingCoordinates = location.boundingCoordinates(distance, radius); boolean meridian180WithinDistance = boundingCoordinates[0].getLongitudeInRadians() > boundingCoordinates[1].getLongitudeInRadians(); java.sql.PreparedStatement statement = connection.prepareStatement( "SELECT * FROM Places WHERE (Lat >= ? AND Lat <= ?) AND (Lon >= ? " + (meridian180WithinDistance ? "OR" : "AND") + " Lon <= ?) AND " + "acos(sin(?) * sin(Lat) + cos(?) * cos(Lat) * cos(Lon - ?)) <= ?"); statement.setDouble(1, boundingCoordinates[0].getLatitudeInRadians()); statement.setDouble(2, boundingCoordinates[1].getLatitudeInRadians()); statement.setDouble(3, boundingCoordinates[0].getLongitudeInRadians()); statement.setDouble(4, boundingCoordinates[1].getLongitudeInRadians()); statement.setDouble(5, location.getLatitudeInRadians()); statement.setDouble(6, location.getLatitudeInRadians()); statement.setDouble(7, location.getLongitudeInRadians()); statement.setDouble(8, distance / radius); return statement.executeQuery(); } public static void main(String[] args) { double earthRadius = 6371.01; GeoLocation myLocation = GeoLocation.fromRadians(1.3963, -0.6981); double distance = 1000; java.sql.Connection connection = ...; java.sql.ResultSet resultSet = findPlacesWithinDistance( earthRadius, myLocation, distance, connection); ...; } }Download GeoLocationDemo.java
Some people have translated the Java code from section 4 into other programming languages and published their results. There are implementations in C#, Clojure, PHP and Python. Note that I have not checked these for correctness.
If you have comments or suggestions for improving this article, please send me e-mail.
This article is licensed under a Creative Commons Attribution License. If you use this article or a part of it, please add a reference to this page (http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates).