Bagaimana cara menghitung kotak pembatas untuk lokasi lintang / bujur tertentu?

Saya telah memberikan lokasi yang ditentukan oleh lintang dan bujur. Sekarang saya ingin menghitung kotak pembatas dalam jarak misalnya 10 kilometer dari titik itu.

Kotak pembatas harus didefinisikan sebagai latmin, lngmin dan latmax, lngmax.

Saya membutuhkan barang ini untuk menggunakan API panoramio .

Apakah seseorang mengetahui rumus cara mendapatkan poin itu?

Sunting: Teman-teman saya mencari rumus / fungsi yang mengambil lat & lng sebagai masukan dan mengembalikan kotak pembatas sebagai latmin & lngmin dan latmax & latmin. Mysql, php, c #, javascript baik-baik saja tetapi juga pseudocode harus baik-baik saja.

Sunting: Saya tidak mencari solusi yang menunjukkan jarak 2 poin

Saya menyarankan untuk memperkirakan secara lokal permukaan bumi sebagai bola dengan radius yang diberikan oleh ellipsoid WGS84 pada garis lintang yang diberikan. Saya menduga bahwa penghitungan yang tepat dari latMin dan latMax akan membutuhkan fungsi elips dan tidak akan menghasilkan peningkatan akurasi yang cukup (WGS84 sendiri merupakan perkiraan).

Implementasi saya mengikuti (Ini ditulis dengan Python; Saya belum mengujinya):

# 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: Kode berikut mengubah (derajat, bilangan prima, detik) menjadi derajat + pecahan derajat, dan sebaliknya (tidak diuji):

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

Saya menulis artikel tentang menemukan koordinat pembatas:

http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates

Artikel tersebut menjelaskan rumus dan juga menyediakan implementasi Java. (Ini juga menunjukkan mengapa rumus Federico untuk bujur min / maks tidak akurat.)

Di sini saya telah mengubah jawaban Federico A. Ramponi menjadi C # untuk siapa pun yang tertarik:

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

Saya menulis fungsi JavaScript yang mengembalikan empat koordinat kotak pembatas persegi, diberi jarak dan sepasang koordinat:

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

Karena saya membutuhkan perkiraan yang sangat kasar, jadi untuk menyaring beberapa dokumen yang tidak perlu dalam kueri elasticsearch, saya menggunakan rumus di bawah ini:

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 diperlukan dari lokasi yang diberikan. Untuk kasus Anda N = 10

Tidak akurat tapi berguna.

Anda sedang mencari rumus ellipsoid.

Tempat terbaik yang saya temukan untuk memulai pengkodean didasarkan pada pustaka Geo :: Ellipsoid dari CPAN. Ini memberi Anda dasar untuk membuat pengujian Anda dari dan untuk membandingkan hasil Anda dengan hasilnya. Saya menggunakannya sebagai dasar untuk pustaka serupa untuk PHP di tempat kerja saya sebelumnya.

Geo :: Ellipsoid

Lihatlah locationmetodenya. Sebut dua kali dan Anda mendapatkan bbox Anda.

Anda tidak memposting bahasa yang Anda gunakan. Mungkin sudah ada pustaka geocoding yang tersedia untuk Anda.

Oh, dan jika Anda belum mengetahuinya sekarang, Google maps menggunakan ellipsoid WGS84.

Ilustrasi penjelasan yang sangat baik dari @Jan Philip Matuschek. (Harap pilih jawabannya, bukan ini; Saya menambahkan ini karena saya meluangkan sedikit waktu untuk memahami jawaban aslinya)

Teknik kotak pembatas untuk mengoptimalkan menemukan tetangga terdekat perlu mendapatkan garis lintang minimum dan maksimum, pasangan bujur, untuk titik P pada jarak d. Semua titik yang berada di luar ini pasti berada pada jarak yang lebih besar dari d dari titik tersebut. Satu hal yang perlu diperhatikan di sini adalah perhitungan garis lintang persimpangan seperti yang disorot dalam penjelasan Jan Philip Matuschek. Garis lintang persimpangan tidak berada pada garis lintang titik P tetapi sedikit diimbangi darinya. Ini adalah bagian yang sering terlewat tetapi penting dalam menentukan garis bujur batas minimum dan maksimum yang benar untuk titik P untuk jarak d. Ini juga berguna dalam verifikasi.

Jarak haversine antara (garis lintang persimpangan, garis bujur tinggi) ke (garis lintang, garis bujur) P sama dengan jarak d.

Inti Python di sini https://gist.github.com/alexcpn/f95ae83a7ee0293a5225

Berikut ini adalah implementasi sederhana menggunakan javascript yang didasarkan pada konversi derajat lintang ke kms dimana 1 degree latitude ~ 111.2 km.

Saya menghitung batas peta dari garis lintang dan bujur tertentu dengan lebar 10 km.

function getBoundsFromLatLng(lat, lng){
     var lat_change = 10/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;
}

Saya mengadaptasi skrip PHP yang saya temukan untuk melakukan hal ini. Anda dapat menggunakannya untuk mencari sudut-sudut kotak di sekitar suatu titik (katakanlah, jarak 20 km). Contoh spesifik saya adalah untuk Google Maps API:

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

Saya sedang mengerjakan masalah kotak pembatas sebagai masalah sampingan untuk menemukan semua titik dalam radius SrcRad dari titik LAT, PANJANG statis. Ada cukup banyak perhitungan yang digunakan

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

untuk menghitung batas bujur, tetapi saya menemukan ini tidak memberikan semua jawaban yang diperlukan. Karena yang benar-benar ingin Anda lakukan adalah

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

Saya tahu, saya tahu jawabannya harus sama, tetapi ternyata tidak. Tampaknya dengan tidak memastikan bahwa saya melakukan (SRCrad / RadEarth) Pertama dan kemudian membaginya dengan bagian Cos saya meninggalkan beberapa titik lokasi.

Setelah Anda mendapatkan semua titik kotak pembatas, jika Anda memiliki fungsi yang menghitung Jarak Titik ke Titik dengan garis lintang, panjang mudah untuk hanya mendapatkan titik-titik yang berada dalam radius jarak tertentu dari titik tetap. Inilah yang saya lakukan. Saya tahu perlu beberapa langkah ekstra, tetapi itu membantu saya

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

Sangat sederhana, cukup buka situs web panoramio dan kemudian buka Peta Dunia dari situs web panoramio. Lalu buka lokasi yang ditentukan yang diperlukan garis lintang dan bujur.

Kemudian Anda menemukan lintang dan bujur di bilah alamat misalnya di alamat ini.

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

lt = 32,739485 => lintang ln = 70,491211 => bujur

widget Panoramio JavaScript API ini membuat kotak pembatas di sekitar pasangan lintang / bujur dan kemudian mengembalikan semua foto dengan batas tersebut.

Jenis lain dari widget Panoramio JavaScript API di mana Anda juga dapat mengubah warna latar belakang dengan contoh dan kode ada di sini .

Itu tidak muncul dalam mood menulis. Itu muncul setelah penerbitan.

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

Di sini saya telah mengubah jawaban Federico A. Ramponi menjadi PHP jika ada yang tertarik:

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

Terima kasih @Fedrico A. untuk implementasi Phyton, saya telah memasukkannya ke dalam kelas kategori Objective C. Disini adalah:

#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

Saya telah mengujinya dan tampaknya bekerja dengan baik. Struct BoundsLocation harus diganti dengan kelas, saya telah menggunakannya hanya untuk membagikannya di sini.

Semua jawaban di atas hanya sebagian yang benar . Khususnya di wilayah seperti Australia, mereka selalu menyertakan tiang dan menghitung persegi panjang yang sangat besar bahkan untuk 10 km.

Khususnya algoritme oleh Jan Philip Matuschek di http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates#UsingIndex menyertakan persegi panjang yang sangat besar dari (-37, -90, -180, 180) untuk hampir setiap titik di Australia. Ini mengenai pengguna besar dalam database dan jarak harus dihitung untuk semua pengguna di hampir setengah negara.

Saya menemukan bahwa Drupal API Earth Algorithm oleh Rochester Institute of Technology bekerja lebih baik di sekitar kutub serta di tempat lain dan jauh lebih mudah untuk diterapkan.

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

Gunakan earth_latitude_rangedan earth_longitude_rangedari algoritma di atas untuk menghitung persegi panjang pembatas

Dan gunakan rumus penghitungan jarak yang didokumentasikan oleh google maps untuk menghitung jarak

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

Untuk mencari berdasarkan kilometer, bukan mil, ganti 3959 dengan 6371. Untuk (Lat, Lng) = (37, -122) dan tabel Markers dengan kolom lat dan lng , rumusnya adalah:

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;

Baca jawaban rinci saya di https://stackoverflow.com/a/45950426/5076414

Inilah jawaban Federico Ramponi di Go. Catatan: tidak ada pengecekan kesalahan :(

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