views:

5192

answers:

7

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: http://www.panoramio.com/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

+3  A: 

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.

http://search.cpan.org/~jgibson/Geo-Ellipsoid-1.12/lib/Geo/Ellipsoid.pm

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.

jcoby
A: 

What You are looking for is called "great circle distance".

Black
you can't actually use a sphere to model the earth without introducing inaccuracy. you have to use an ellipsoid- and there are a dozen or so to choose from. wgs84 seeming to be the most common.
jcoby
But for a 10km distance you can probably get away with a sphere
Martin Beckett
a formula would be great ;)
Michal
This probably happened in an edit after the original question was posted, but the question is after a function that converts latitude, longitude, radius to a bounding box.
mjs
+6  A: 

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))
Federico Ramponi
As pointed out in the documentation of the suggested CPAN library, this makes sense only for halfSide <= 10km.
Federico Ramponi
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.
Doug McClean
A: 

a formula would be great ;)

susheel
+1  A: 

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

Jan Philip Matuschek
A: 

Can you please tell me what does latMin, lonMin, latMax and lonMax represent? According to me this is the representation: latMax: North latMin: South lonMax: East and lonMin: West. Please correct me if I am wrong.

A: 

Is this the following representation:

latMax: North latMin: South lonMax: East lonMin: West