How to calculate the bounding box for a given lat/lng location?
I suggest to approximate locally the Earth surface as a sphere with radius given by the WGS84 ellipsoid at the given latitude. I suspect that the exact computation of latMin and latMax would require elliptic functions and would not yield an appreciable increase in accuracy (WGS84 is itself an approximation).
My implementation follows (It's written in Python; I have not tested it):
# degrees to radians
def deg2rad(degrees):
return math.pi*degrees/180.0
# radians to degrees
def rad2deg(radians):
return 180.0*radians/math.pi
# Semi-axes of WGS-84 geoidal reference
WGS84_a = 6378137.0 # Major semiaxis [m]
WGS84_b = 6356752.3 # Minor semiaxis [m]
# Earth radius at a given latitude, according to the WGS-84 ellipsoid [m]
def WGS84EarthRadius(lat):
# http://en.wikipedia.org/wiki/Earth_radius
An = WGS84_a*WGS84_a * math.cos(lat)
Bn = WGS84_b*WGS84_b * math.sin(lat)
Ad = WGS84_a * math.cos(lat)
Bd = WGS84_b * math.sin(lat)
return math.sqrt( (An*An + Bn*Bn)/(Ad*Ad + Bd*Bd) )
# Bounding box surrounding the point at given coordinates,
# assuming local approximation of Earth surface as a sphere
# of radius given by WGS84
def boundingBox(latitudeInDegrees, longitudeInDegrees, halfSideInKm):
lat = deg2rad(latitudeInDegrees)
lon = deg2rad(longitudeInDegrees)
halfSide = 1000*halfSideInKm
# Radius of Earth at given latitude
radius = WGS84EarthRadius(lat)
# Radius of the parallel at given latitude
pradius = radius*math.cos(lat)
latMin = lat - halfSide/radius
latMax = lat + halfSide/radius
lonMin = lon - halfSide/pradius
lonMax = lon + halfSide/pradius
return (rad2deg(latMin), rad2deg(lonMin), rad2deg(latMax), rad2deg(lonMax))
EDIT: The following code converts (degrees, primes, seconds) to degrees + fractions of a degree, and vice versa (not tested):
def dps2deg(degrees, primes, seconds):
return degrees + primes/60.0 + seconds/3600.0
def deg2dps(degrees):
intdeg = math.floor(degrees)
primes = (degrees - intdeg)*60.0
intpri = math.floor(primes)
seconds = (primes - intpri)*60.0
intsec = round(seconds)
return (int(intdeg), int(intpri), int(intsec))
I wrote an article about finding the bounding coordinates:
http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates
The article explains the formulae and also provides a Java implementation. (It also shows why Federico's formula for the min/max longitude is inaccurate.)
Here I have converted Federico A. Ramponi's answer to C# for anybody interested:
public class MapPoint
{
public double Longitude { get; set; } // In Degrees
public double Latitude { get; set; } // In Degrees
}
public class BoundingBox
{
public MapPoint MinPoint { get; set; }
public MapPoint MaxPoint { get; set; }
}
// Semi-axes of WGS-84 geoidal reference
private const double WGS84_a = 6378137.0; // Major semiaxis [m]
private const double WGS84_b = 6356752.3; // Minor semiaxis [m]
// 'halfSideInKm' is the half length of the bounding box you want in kilometers.
public static BoundingBox GetBoundingBox(MapPoint point, double halfSideInKm)
{
// Bounding box surrounding the point at given coordinates,
// assuming local approximation of Earth surface as a sphere
// of radius given by WGS84
var lat = Deg2rad(point.Latitude);
var lon = Deg2rad(point.Longitude);
var halfSide = 1000 * halfSideInKm;
// Radius of Earth at given latitude
var radius = WGS84EarthRadius(lat);
// Radius of the parallel at given latitude
var pradius = radius * Math.Cos(lat);
var latMin = lat - halfSide / radius;
var latMax = lat + halfSide / radius;
var lonMin = lon - halfSide / pradius;
var lonMax = lon + halfSide / pradius;
return new BoundingBox {
MinPoint = new MapPoint { Latitude = Rad2deg(latMin), Longitude = Rad2deg(lonMin) },
MaxPoint = new MapPoint { Latitude = Rad2deg(latMax), Longitude = Rad2deg(lonMax) }
};
}
// degrees to radians
private static double Deg2rad(double degrees)
{
return Math.PI * degrees / 180.0;
}
// radians to degrees
private static double Rad2deg(double radians)
{
return 180.0 * radians / Math.PI;
}
// Earth radius at a given latitude, according to the WGS-84 ellipsoid [m]
private static double WGS84EarthRadius(double lat)
{
// http://en.wikipedia.org/wiki/Earth_radius
var An = WGS84_a * WGS84_a * Math.Cos(lat);
var Bn = WGS84_b * WGS84_b * Math.Sin(lat);
var Ad = WGS84_a * Math.Cos(lat);
var Bd = WGS84_b * Math.Sin(lat);
return Math.Sqrt((An*An + Bn*Bn) / (Ad*Ad + Bd*Bd));
}
Since I needed a very rough estimate, so to filter out some needless documents in an elasticsearch query, I employed the below formula:
Min.lat = Given.Lat - (0.009 x N)
Max.lat = Given.Lat + (0.009 x N)
Min.lon = Given.lon - (0.009 x N)
Max.lon = Given.lon + (0.009 x N)
N = kms required form the given location. For your case N=10
Not accurate but handy.
I wrote a JavaScript function that returns the four coordinates of a square bounding box, given a distance and a pair of coordinates:
'use strict';
/**
* @param {number} distance - distance (km) from the point represented by centerPoint
* @param {array} centerPoint - two-dimensional array containing center coords [latitude, longitude]
* @description
* Computes the bounding coordinates of all points on the surface of a sphere
* that has a great circle distance to the point represented by the centerPoint
* argument that is less or equal to the distance argument.
* Technique from: Jan Matuschek <http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates>
* @author Alex Salisbury
*/
getBoundingBox = function (centerPoint, distance) {
var MIN_LAT, MAX_LAT, MIN_LON, MAX_LON, R, radDist, degLat, degLon, radLat, radLon, minLat, maxLat, minLon, maxLon, deltaLon;
if (distance < 0) {
return 'Illegal arguments';
}
// helper functions (degrees<–>radians)
Number.prototype.degToRad = function () {
return this * (Math.PI / 180);
};
Number.prototype.radToDeg = function () {
return (180 * this) / Math.PI;
};
// coordinate limits
MIN_LAT = (-90).degToRad();
MAX_LAT = (90).degToRad();
MIN_LON = (-180).degToRad();
MAX_LON = (180).degToRad();
// Earth's radius (km)
R = 6378.1;
// angular distance in radians on a great circle
radDist = distance / R;
// center point coordinates (deg)
degLat = centerPoint[0];
degLon = centerPoint[1];
// center point coordinates (rad)
radLat = degLat.degToRad();
radLon = degLon.degToRad();
// minimum and maximum latitudes for given distance
minLat = radLat - radDist;
maxLat = radLat + radDist;
// minimum and maximum longitudes for given distance
minLon = void 0;
maxLon = void 0;
// define deltaLon to help determine min and max longitudes
deltaLon = Math.asin(Math.sin(radDist) / Math.cos(radLat));
if (minLat > MIN_LAT && maxLat < MAX_LAT) {
minLon = radLon - deltaLon;
maxLon = radLon + deltaLon;
if (minLon < MIN_LON) {
minLon = minLon + 2 * Math.PI;
}
if (maxLon > MAX_LON) {
maxLon = maxLon - 2 * Math.PI;
}
}
// a pole is within the given distance
else {
minLat = Math.max(minLat, MIN_LAT);
maxLat = Math.min(maxLat, MAX_LAT);
minLon = MIN_LON;
maxLon = MAX_LON;
}
return [
minLon.radToDeg(),
minLat.radToDeg(),
maxLon.radToDeg(),
maxLat.radToDeg()
];
};