How to calculate the bounding box for a given lat/lng location?

Question

I have given a location defined by latitude and longitude. Now i want to calculate a bounding box within e.g. 10 kilometers of that point.

The bounding box should be defined as latmin, lngmin and latmax, lngmax.

I need this stuff in order to use the panoramio API.

Does someone know the formula of how to get thos points?

Edit: Guys i am looking for a formula/function which takes lat & lng as input and returns a bounding box as latmin & lngmin and latmax & latmin. Mysql, php, c#, javascript is fine but also pseudocode should be okay.

Edit: I am not looking for a solution which shows me the distance of 2 points

Exosmosis

If you are using a geodatabase somewhere, they surely have a bounding box calculation integrated. You could even go check the source of PostGIS/GEOS, for example. – Ledesma

Answer 1

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

Tohubohu

As pointed out in the documentation of the suggested CPAN library, this makes sense only for halfSide <= 10km. – Tohubohu

Does this work near the poles? It doesn't seem to, because it looks like it ends up with latMin < -pi (for the south pole) or latMax > pi (for the north pole)? I think when you are within halfSide of a pole you need to return a bounding box that includes all longitudes and the latitudes computed normally for the side away from the pole and at the pole on the side near the pole. – Xerox

This does not work around the poles either. JanMatuschek.de/LatitudeLongitudeBoundingCoordinates gives the "correct algorithm" – Pleasance

Here is a PHP implementation from the specification found at JanMatuschek.de: github.com/anthonymartin/GeoLocation.class.php – Chic

I have added a C# implementation of this answer down below. – Papiamento

@FedericoA.Ramponi what is the haldSideinKm here? don't understand... what I must to pass in this argyments, the radius between two points in map or what? – Bushnell

Answer 2

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.)

Sylvestersylvia

I've made a PHP port of your GeoLocation class. It can be found here: pastie.org/5416584 – Chic

I've uploaded it on github now: github.com/anthonymartin/GeoLocation.class.php – Chic

Does this even answer the question? If we only have 1 starting point we can't calculate the great circle distance as is done in this code, that requires two lat,long locations. – Maniemanifest

there is a bad code in your C# variant, for e.g.: public override string ToString(), it's very bad to override such a global method only for one purpose, better just to add another method, then overriding standard method, which can be used in other parts of application, not for the gis exact... – Bushnell

Here's an updated link to the PHP port of Jan's GeoLocaiton class: github.com/anthonymartin/GeoLocation.php – Chic

Thanks. Can you also explain what would be the error percentage if we only use the bounding coordinates and not the distance formula ? – Sealey

Your articles has been rather helpful to me in more than one occasion, thanks a lot for writing and sharing it. – Execution

Answer 3

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));
}

Affiliate

Thanks, this work for me. Had to test the code by hand, didn't know how to write a unit test for this but it generates accurate results to the degree of accuracy i need – Maniemanifest

what is the haldSideinKm here? don't understand... what I must to pass in this argyments, the radius between two points in map or what? – Bushnell

@GeloVolro: That's the half length of the bounding box you want. – Papiamento

You don't necessarily have to write your own MapPoint class. There's a GeoCoordinate class in System.Device.Location that takes Latitude and Longitude as parameters. – Orbital

@Ε Г И І И О could you have a look at my question please, you might be able to help: #35390334 – Rouble

This works beautifully. I really appreciate the C# port. – Aeroplane

Answer 4

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.

Endow

Indeed, not accurate but still useful and very easy to implement. – Chinaman

what does 0.009 mean in this formula? – Irradiance

Answer 5

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()
  ];
};

Foxhound

This code does not work, at all. I mean, even after fixing up the obvious errors like minLon = void 0; and maxLon = MAX_LON; it still doesn't work. – Male

@aroth, I just tested it and had no issue. Remember that centerPoint argument is an array consisting of two coordinates. E.g., getBoundingBox([42.2, 34.5], 50) — void 0 is the CoffeeScript output for "undefined" and won't affect the codes ability to run. – Foxhound

this code doesn't work. degLat.degToRad is not a function – Calcite

Code worked in Node and Chrome as-is, until I put it inside a project I'm working on and started getting "degToRad is not a function" errors. Never found out why but Number.prototype. isn't a good idea for a utility function like this so I converted those to normal local functions. Also it's important to note that the box returned is [LNG,LAT,LNG,LAT] instead of [LAT,LNG,LAT,LNG]. I modified the return function when I used this to avoid confusion. – Tlingit

Answer 6

Here is an simple implementation using javascript which is based on the conversion of latitude degree to kms where 1 degree latitude ~ 111.2 km.

I am calculating bounds of the map from a given latitude, longitude and radius in kilometers.

function getBoundsFromLatLng(lat, lng, radiusInKm){
     var lat_change = radiusInKm/111.2;
     var lon_change = Math.abs(Math.cos(lat*(Math.PI/180)));
     var bounds = { 
         lat_min : lat - lat_change,
         lon_min : lng - lon_change,
         lat_max : lat + lat_change,
         lon_max : lng + lon_change
     };
     return bounds;
}

Ozellaozen

Answer 7

Illustration of @Jan Philip Matuschek excellent explanation.(Please up-vote his answer, not this; I am adding this as I took a little time in understanding the original answer)

The bounding box technique of optimizing of finding nearest neighbors would need to derive the minimum and maximum latitude,longitude pairs, for a point P at distance d . All points that fall outside these are definitely at a distance greater than d from the point. One thing to note here is the calculation of latitude of intersection as is highlighted in Jan Philip Matuschek explanation. The latitude of intersection is not at the latitude of point P but slightly offset from it. This is a often missed but important part in determining the correct minimum and maximum bounding longitude for point P for the distance d.This is also useful in verification.

The haversine distance between (latitude of intersection,longitude high) to (latitude,longitude) of P is equal to distance d.

Python gist here https://gist.github.com/alexcpn/f95ae83a7ee0293a5225

Steersman

Answer 8

You're looking for an ellipsoid formula.

The best place I've found to start coding is based on the Geo::Ellipsoid library from CPAN. It gives you a baseline to create your tests off of and to compare your results with its results. I used it as the basis for a similar library for PHP at my previous employer.

Geo::Ellipsoid

Take a look at the location method. Call it twice and you've got your bbox.

You didn't post what language you were using. There may already be a geocoding library available for you.

Oh, and if you haven't figured it out by now, Google maps uses the WGS84 ellipsoid.

Wamsley

Answer 9

I adapted a PHP script I found to do just this. You can use it to find the corners of a box around a point (say, 20 km out). My specific example is for Google Maps API:

http://www.richardpeacock.com/blog/2011/11/draw-box-around-coordinate-google-maps-based-miles-or-kilometers

Murther

-1 What the OP is looking for is: given a reference point (lat, lon) and a distance, find the smallest box such that all points that are <= "distance" away from the reference point are not outside the box. Your box has its corners "distance" away from the reference point and is thus too small. Example: the point that is "distance" due north is well outside your box. – Gauche

Well, by chance, it's exactly what I just needed. So thank you, even if it doesn't quite answer this question :) – Ortrud

Well, I'm glad it could help somebody! – Murther

Answer 10

This is javascript code for getting bounding box co-ordinates based on lat/long and distance. tested and working fine.

Number.prototype.degreeToRadius = function () {
    return this * (Math.PI / 180);
};

Number.prototype.radiusToDegree = function () {
    return (180 * this) / Math.PI;
};

function getBoundingBox(fsLatitude, fsLongitude, fiDistanceInKM) {

    if (fiDistanceInKM == null || fiDistanceInKM == undefined || fiDistanceInKM == 0)
        fiDistanceInKM = 1;
    
    var MIN_LAT, MAX_LAT, MIN_LON, MAX_LON, ldEarthRadius, ldDistanceInRadius, lsLatitudeInDegree, lsLongitudeInDegree,
        lsLatitudeInRadius, lsLongitudeInRadius, lsMinLatitude, lsMaxLatitude, lsMinLongitude, lsMaxLongitude, deltaLon;
    
    // coordinate limits
    MIN_LAT = (-90).degreeToRadius();
    MAX_LAT = (90).degreeToRadius();
    MIN_LON = (-180).degreeToRadius();
    MAX_LON = (180).degreeToRadius();

    // Earth's radius (km)
    ldEarthRadius = 6378.1;

    // angular distance in radians on a great circle
    ldDistanceInRadius = fiDistanceInKM / ldEarthRadius;

    // center point coordinates (deg)
    lsLatitudeInDegree = fsLatitude;
    lsLongitudeInDegree = fsLongitude;

    // center point coordinates (rad)
    lsLatitudeInRadius = lsLatitudeInDegree.degreeToRadius();
    lsLongitudeInRadius = lsLongitudeInDegree.degreeToRadius();

    // minimum and maximum latitudes for given distance
    lsMinLatitude = lsLatitudeInRadius - ldDistanceInRadius;
    lsMaxLatitude = lsLatitudeInRadius + ldDistanceInRadius;

    // minimum and maximum longitudes for given distance
    lsMinLongitude = void 0;
    lsMaxLongitude = void 0;

    // define deltaLon to help determine min and max longitudes
    deltaLon = Math.asin(Math.sin(ldDistanceInRadius) / Math.cos(lsLatitudeInRadius));

    if (lsMinLatitude > MIN_LAT && lsMaxLatitude < MAX_LAT) {
        lsMinLongitude = lsLongitudeInRadius - deltaLon;
        lsMaxLongitude = lsLongitudeInRadius + deltaLon;
        if (lsMinLongitude < MIN_LON) {
            lsMinLongitude = lsMinLongitude + 2 * Math.PI;
        }
        if (lsMaxLongitude > MAX_LON) {
            lsMaxLongitude = lsMaxLongitude - 2 * Math.PI;
        }
    }

    // a pole is within the given distance
    else {
        lsMinLatitude = Math.max(lsMinLatitude, MIN_LAT);
        lsMaxLatitude = Math.min(lsMaxLatitude, MAX_LAT);
        lsMinLongitude = MIN_LON;
        lsMaxLongitude = MAX_LON;
    }

    return [
        lsMinLatitude.radiusToDegree(),
        lsMinLongitude.radiusToDegree(),
        lsMaxLatitude.radiusToDegree(),
        lsMaxLongitude.radiusToDegree()
    ];
};

Use getBoundingBox function like below to draw a bounding box.

var lsRectangleLatLong = getBoundingBox(parseFloat(latitude), parseFloat(longitude), lsDistance);
if (lsRectangleLatLong != null && lsRectangleLatLong != undefined) {
    latLngArr.push({ lat: lsRectangleLatLong[0], lng: lsRectangleLatLong[1] });
    latLngArr.push({ lat: lsRectangleLatLong[0], lng: lsRectangleLatLong[3] });
    latLngArr.push({ lat: lsRectangleLatLong[2], lng: lsRectangleLatLong[3] });
    latLngArr.push({ lat: lsRectangleLatLong[2], lng: lsRectangleLatLong[1] });
}

Maximomaximum

Answer 11

I was working on the bounding box problem as a side issue to finding all the points within SrcRad radius of a static LAT, LONG point. There have been quite a few calculations that use

maxLon = $lon + rad2deg($rad/$R/cos(deg2rad($lat)));
minLon = $lon - rad2deg($rad/$R/cos(deg2rad($lat)));

to calculate the longitude bounds, but I found this to not give all the answers that were needed. Because what you really want to do is

(SrcRad/RadEarth)/cos(deg2rad(lat))

I know, I know the answer should be the same, but I found that it wasn't. It appeared that by not making sure I was doing the (SRCrad/RadEarth) First and then dividing by the Cos part I was leaving out some location points.

After you get all your bounding box points, if you have a function that calculates the Point to Point Distance given lat, long it is easy to only get those points that are a certain distance radius from the fixed point. Here is what I did. I know it took a few extra steps but it helped me

-- GLOBAL Constants
gc_pi CONSTANT REAL := 3.14159265359;  -- Pi

-- Conversion Factor Constants
gc_rad_to_degs          CONSTANT NUMBER := 180/gc_pi; -- Conversion for Radians to Degrees 180/pi
gc_deg_to_rads          CONSTANT NUMBER := gc_pi/180; --Conversion of Degrees to Radians

lv_stat_lat    -- The static latitude point that I am searching from 
lv_stat_long   -- The static longitude point that I am searching from 

-- Angular radius ratio in radians
lv_ang_radius := lv_search_radius / lv_earth_radius;
lv_bb_maxlat := lv_stat_lat + (gc_rad_to_deg * lv_ang_radius);
lv_bb_minlat := lv_stat_lat - (gc_rad_to_deg * lv_ang_radius);

--Here's the tricky part, accounting for the Longitude getting smaller as we move up the latitiude scale
-- I seperated the parts of the equation to make it easier to debug and understand
-- I may not be a smart man but I know what the right answer is... :-)

lv_int_calc := gc_deg_to_rads * lv_stat_lat;
lv_int_calc := COS(lv_int_calc);
lv_int_calc := lv_ang_radius/lv_int_calc;
lv_int_calc := gc_rad_to_degs*lv_int_calc;

lv_bb_maxlong := lv_stat_long + lv_int_calc;
lv_bb_minlong := lv_stat_long - lv_int_calc;

-- Now select the values from your location datatable 
SELECT *  FROM (
SELECT cityaliasname, city, state, zipcode, latitude, longitude, 
-- The actual distance in miles
spherecos_pnttopntdist(lv_stat_lat, lv_stat_long, latitude, longitude, 'M') as miles_dist    
FROM Location_Table 
WHERE latitude between lv_bb_minlat AND lv_bb_maxlat
AND   longitude between lv_bb_minlong and lv_bb_maxlong)
WHERE miles_dist <= lv_limit_distance_miles
order by miles_dist
;

Millennial

Answer 12

Here I have converted Federico A. Ramponi's answer to PHP if anybody is interested:

<?php
# deg2rad and rad2deg are already within PHP

# 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]
function WGS84EarthRadius($lat)
{
    global $WGS84_a, $WGS84_b;
    
    $an = $WGS84_a * $WGS84_a * cos($lat);
    $bn = $WGS84_b * $WGS84_b * sin($lat);
    $ad = $WGS84_a * cos($lat);
    $bd = $WGS84_b * sin($lat);
    
    return 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
function 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*cos($lat);
    
    $latMin = $lat - $halfSide / $radius;
    $latMax = $lat + $halfSide / $radius;
    $lonMin = $lon - $halfSide / $pradius;
    $lonMax = $lon + $halfSide / $pradius;
    
    return array(rad2deg($latMin), rad2deg($lonMin), rad2deg($latMax), rad2deg($lonMax));
}
?>

Clubhouse

For PHP 7.4, I had to slightly mod the code to get it to work: ` # Semi-axes of WGS-84 geoidal reference $WGS84_a = floatval( 6378137.0 ); # Major semiaxis [m] $WGS84_b = floatval( 6356752.3 ); # Minor semiaxis [m] ` – Cliquish

Answer 13

Thanks @Fedrico A. for the Phyton implementation, I have ported it into a Objective C category class. Here is:

#import "LocationService+Bounds.h"

//Semi-axes of WGS-84 geoidal reference
const double WGS84_a = 6378137.0; //Major semiaxis [m]
const double WGS84_b = 6356752.3; //Minor semiaxis [m]

@implementation LocationService (Bounds)

struct BoundsLocation {
    double maxLatitude;
    double minLatitude;
    double maxLongitude;
    double minLongitude;
};

+ (struct BoundsLocation)locationBoundsWithLatitude:(double)aLatitude longitude:(double)aLongitude maxDistanceKm:(NSInteger)aMaxKmDistance {
    return [self boundingBoxWithLatitude:aLatitude longitude:aLongitude halfDistanceKm:aMaxKmDistance/2];
}

#pragma mark - Algorithm 

+ (struct BoundsLocation)boundingBoxWithLatitude:(double)aLatitude longitude:(double)aLongitude halfDistanceKm:(double)aDistanceKm {
    double radianLatitude = [self degreesToRadians:aLatitude];
    double radianLongitude = [self degreesToRadians:aLongitude];
    double halfDistanceMeters = aDistanceKm*1000;
    
    
    double earthRadius = [self earthRadiusAtLatitude:radianLatitude];
    double parallelRadius = earthRadius*cosl(radianLatitude);
    
    double radianMinLatitude = radianLatitude - halfDistanceMeters/earthRadius;
    double radianMaxLatitude = radianLatitude + halfDistanceMeters/earthRadius;
    double radianMinLongitude = radianLongitude - halfDistanceMeters/parallelRadius;
    double radianMaxLongitude = radianLongitude + halfDistanceMeters/parallelRadius;
    
    struct BoundsLocation bounds;
    bounds.minLatitude = [self radiansToDegrees:radianMinLatitude];
    bounds.maxLatitude = [self radiansToDegrees:radianMaxLatitude];
    bounds.minLongitude = [self radiansToDegrees:radianMinLongitude];
    bounds.maxLongitude = [self radiansToDegrees:radianMaxLongitude];

    return bounds;
}

+ (double)earthRadiusAtLatitude:(double)aRadianLatitude {
    double An = WGS84_a * WGS84_a * cosl(aRadianLatitude);
    double Bn = WGS84_b * WGS84_b * sinl(aRadianLatitude);
    double Ad = WGS84_a * cosl(aRadianLatitude);
    double Bd = WGS84_b * sinl(aRadianLatitude);
    return sqrtl( ((An * An) + (Bn * Bn))/((Ad * Ad) + (Bd * Bd)) );
}

+ (double)degreesToRadians:(double)aDegrees {
    return M_PI*aDegrees/180.0;
}

+ (double)radiansToDegrees:(double)aRadians {
    return 180.0*aRadians/M_PI;
}

@end

I have tested it and seems be working nice. Struct BoundsLocation should be replaced by a class, I have used it just to share it here.

Vey

Answer 14

Here is the c++ implementation.

#include <array>
#include <cmath>

double degrees_to_radian(double deg)
{
    return deg * M_PI / 180.0;
}

double radian_to_degrees(double rad)
{
    return rad * 180.0 / M_PI ;
}

double WGS84EarthRadius(double lat_radians){
    static const double WGS84_a = 6378137.0;  // Major semiaxis [m]
    static const double WGS84_b = 6356752.3;  // Minor semiaxis [m]
    double An = WGS84_a*WGS84_a * cos(lat_radians);
    double Bn = WGS84_b*WGS84_b * sin(lat_radians);
    double Ad = WGS84_a * cos(lat_radians);
    double Bd = WGS84_b * sin(lat_radians);

    return std::sqrt( (An*An + Bn*Bn)/(Ad*Ad + Bd*Bd) );
}

std::array<double, 4> BoundingBox(double latitude, double longitude, uint32_t range_meters)
{
    double lat = degrees_to_radian(latitude);
    double lon = degrees_to_radian(longitude);
    
    uint32_t halfSide = range_meters;

    double radius =  WGS84EarthRadius(lat);
    double pradius = radius*cos(lat);

    double latMin = lat - halfSide/radius;
    double latMax = lat + halfSide/radius;
    double lonMin = lon - halfSide/pradius;
    double lonMax = lon + halfSide/pradius;

    return {radian_to_degrees(latMin), radian_to_degrees(lonMin), radian_to_degrees(latMax) , radian_to_degrees(lonMax)};
}

Whittle

Answer 15

It is very simple just go to panoramio website and then open World Map from panoramio website.Then go to specified location whichs latitude and longitude required.

Then you found latitude and longitude in address bar for example in this address.

http://www.panoramio.com/map#lt=32.739485&ln=70.491211&z=9&k=1&a=1&tab=1&pl=all

lt=32.739485 =>latitude ln=70.491211 =>longitude

this Panoramio JavaScript API widget create a bounding box around a lat/long pair and then returning all photos with in those bounds.

Another type of Panoramio JavaScript API widget in which you can also change background color with example and code is here.

It does not show in composing mood.It show after publishing.

<div dir="ltr" style="text-align: center;" trbidi="on">
<script src="https://ssl.panoramio.com/wapi/wapi.js?v=1&amp;hl=en"></script>
<div id="wapiblock" style="float: right; margin: 10px 15px"></div>
<script type="text/javascript">
var myRequest = {
  'tag': 'kahna',
  'rect': {'sw': {'lat': -30, 'lng': 10.5}, 'ne': {'lat': 50.5, 'lng': 30}}
};
  var myOptions = {
  'width': 300,
  'height': 200
};
var wapiblock = document.getElementById('wapiblock');
var photo_widget = new panoramio.PhotoWidget('wapiblock', myRequest, myOptions);
photo_widget.setPosition(0);
</script>
</div>

Jaquenette

Answer 16

All of the above answer are only partially correct. Specially in region like Australia, they always include pole and calculate a very large rectangle even for 10kms.

Specially the algorithm by Jan Philip Matuschek at http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates#UsingIndex included a very large rectangle from (-37, -90, -180, 180) for almost every point in Australia. This hits a large users in database and distance have to be calculated for all of the users in almost half the country.

I found that the Drupal API Earth Algorithm by Rochester Institute of Technology works better around pole as well as elsewhere and is much easier to implement.

https://www.rit.edu/drupal/api/drupal/sites%21all%21modules%21location%21earth.inc/7.54

Use earth_latitude_range and earth_longitude_range from the above algorithm for calculating bounding rectangle

And use the distance calculation formula documented by google maps to calculate distance

https://developers.google.com/maps/solutions/store-locator/clothing-store-locator#outputting-data-as-xml-using-php

To search by kilometers instead of miles, replace 3959 with 6371. For (Lat, Lng) = (37, -122) and a Markers table with columns lat and lng, the formula is:

SELECT id, ( 3959 * acos( cos( radians(37) ) * cos( radians( lat ) ) * cos( radians( lng ) - radians(-122) ) + sin( radians(37) ) * sin( radians( lat ) ) ) ) AS distance FROM markers HAVING distance < 25 ORDER BY distance LIMIT 0 , 20;

Read my detailed answer at https://mcmap.net/q/48546/-calculating-bounding-box-a-certain-distance-away-from-a-lat-long-coordinate-in-java

Needs

Answer 17

Here is Federico Ramponi's answer in Go. Note: no error-checking :(

import (
    "math"
)

// Semi-axes of WGS-84 geoidal reference
const (
    // Major semiaxis (meters)
    WGS84A = 6378137.0
    // Minor semiaxis (meters)
    WGS84B = 6356752.3
)

// BoundingBox represents the geo-polygon that encompasses the given point and radius
type BoundingBox struct {
    LatMin float64
    LatMax float64
    LonMin float64
    LonMax float64
}

// Convert a degree value to radians
func deg2Rad(deg float64) float64 {
    return math.Pi * deg / 180.0
}

// Convert a radian value to degrees
func rad2Deg(rad float64) float64 {
    return 180.0 * rad / math.Pi
}

// Get the Earth's radius in meters at a given latitude based on the WGS84 ellipsoid
func getWgs84EarthRadius(lat float64) float64 {
    an := WGS84A * WGS84A * math.Cos(lat)
    bn := WGS84B * WGS84B * math.Sin(lat)

    ad := WGS84A * math.Cos(lat)
    bd := WGS84B * math.Sin(lat)

    return math.Sqrt((an*an + bn*bn) / (ad*ad + bd*bd))
}

// GetBoundingBox returns a BoundingBox encompassing the given lat/long point and radius
func GetBoundingBox(latDeg float64, longDeg float64, radiusKm float64) BoundingBox {
    lat := deg2Rad(latDeg)
    lon := deg2Rad(longDeg)
    halfSide := 1000 * radiusKm

    // Radius of Earth at given latitude
    radius := getWgs84EarthRadius(lat)

    pradius := radius * math.Cos(lat)

    latMin := lat - halfSide/radius
    latMax := lat + halfSide/radius
    lonMin := lon - halfSide/pradius
    lonMax := lon + halfSide/pradius

    return BoundingBox{
        LatMin: rad2Deg(latMin),
        LatMax: rad2Deg(latMax),
        LonMin: rad2Deg(lonMin),
        LonMax: rad2Deg(lonMax),
    }
}

Athalia