geocodinglatitude-longitudepanoramio

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


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


Solution

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