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

Viewed 129649

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

16 Answers

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

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.

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

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),
    }
}
Related