如何计算给定的lat/lng位置的包围盒?

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.

我已经给出了一个由经度和纬度定义的位置。现在我想计算一个包围盒在这个点的10千米内。

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

边界框应该被定义为latmin, lngmin和latmax, lngmax。

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

我需要这些东西来使用全景API。

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

有人知道如何得到thos点的公式吗?

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.

编辑:我正在寻找一个公式/函数，它以lat & lng作为输入，并返回一个边界框作为latmin & lngmin和latmax & latmin。Mysql, php, c#， javascript很好，但是伪代码也应该没问题。

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

编辑:我不是在寻找一个能显示2点距离的解。

15 个解决方案

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

我建议将地球表面近似为一个球体，在给定纬度上由WGS84椭球体给出。我猜想，对latMin和latMax的精确计算需要椭圆函数，而且精度不会有明显提高(WGS84本身就是一个近似)。

My implementation follows (It's written in Python; I have not tested it):

我的实现如下(它是用Python编写的;我没有测试它):

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

#2

I wrote an article about finding the bounding coordinates:

我写了一篇关于寻找边界坐标的文章:

http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates

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

本文解释了公式，并提供了Java实现。(这也说明了为什么费德里科的最小/最大经度公式不准确。)

#3

Here I have converted Federico A. Ramponi's answer to C# for anybody interested:

在这里，我已经将费德里科A. Ramponi的答案转换为对任何感兴趣的人的c#:

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

#4

I wrote a JavaScript function that returns the four coordinates of a square bounding box, given a distance and a pair of coordinates:

我编写了一个JavaScript函数，它返回一个正方形边界框的四个坐标，给定一个距离和一对坐标:

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

#5

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::来自CPAN的椭圆体库。它为您提供了创建测试的基线，并将结果与结果进行比较。我在以前的雇主那里用它作为一个类似的PHP库的基础。

Geo::Ellipsoid

地理:椭球

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

看一下位置方法。叫它两次，你就得到了你的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.

哦，如果你还没有算出来，谷歌地图使用的是WGS84椭球。

#6

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:

我修改了一个PHP脚本。你可以用它来找到一个盒子的角(比如说，20公里外)。我的具体例子是谷歌地图API:

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

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

#7

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

N = kms所需的位置。对你的案子N = 10

Not accurate but handy.

不准确而方便。

#8

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.

这是一个使用javascript的简单实现，它是基于对kms的纬度的转换，其中1度纬度~ 111.2公里。

I am calculating bounds of the map from a given latitude and longitude with 10km width.

我正在计算地图的边界，从给定的经纬度和10公里宽。

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

#9

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)

@Jan Philip Matuschek出色的解释说明。(请举手表决他的答案，而不是这个;我加了这个，因为我花了一点时间来理解原始答案

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.

寻找最近邻的边界箱技术需要得到最小和最大纬度，经度对，在距离d处的点P。所有落在这些点之外的点肯定是大于d的点。这里需要注意的一点是，在Jan Philip Matuschek的解释中强调了交叉的纬度的计算。交点的纬度不在点P的纬度，而是稍微偏离点。这是在确定距离d点P的正确最小和最大跳跃经度时经常忽略但很重要的部分。这在验证中也很有用。

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

在(交点的纬度，经度高)到(纬度，经度)之间的距离等于距离d。

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

Python要点这里https://gist.github.com/alexcpn/f95ae83a7ee0293a5225

#10

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

我正在研究边界框的问题，作为一个次要的问题来找到一个静态LAT的SrcRad半径内的所有点。有相当多的计算使用。

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.

我知道，我知道答案应该是一样的，但我发现事实并非如此。似乎没有确保我先做(SRCrad/RadEarth)然后除以Cos的部分，我省略了一些位置点。

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
;

#11

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

http://www.panoramio.com/map lt = 32.739485 ln = 70.491211 z = 9 = 1和大部分= 1选项卡= 1 pl =

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

lt = 32.739485 = >纬度ln = 70.491211 = >经度

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

这个全景JavaScript API小部件在一个lat/长对上创建一个包围盒，然后在这些边界中返回所有的照片。

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

另一种类型的全景JavaScript API小部件，您还可以通过示例和代码更改背景颜色。

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>

#12

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

这里，我已经将费德里科A. Ramponi的答案转换为PHP，如果有人感兴趣的话:

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

#13

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

感谢@Fedrico . for the Phyton implementation，我将它移植到Objective C类中。这里是:

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

我已经测试过了，而且看起来效果不错。Struct BoundsLocation应该被一个类替换，我用它来在这里共享它。

#14

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.

以上所有答案都只是部分正确。特别是在澳大利亚这样的地区，它们总是包括极点，即使是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.

这一算法由Jan Philip Matuschek在http://janmatuschek.de/latitudelongitudeboundingindex中提出，该算法包含了一个非常大的矩形(-37，-90，-180,180)，几乎适用于澳大利亚的每一个点。这对数据库中的大型用户产生了影响，在全国几乎一半的用户中都需要计算距离。

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.

我发现罗彻斯特理工学院的Drupal API地球算法在杆子和其他地方更好，而且更容易实现。

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

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

从上面的算法中使用earth_latitude_range和earth_de_range计算边界矩形。

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

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:

用公里代替英里来搜索，用6371代替3959。For (Lat, Lng) =(37， -122)和一个有列Lat和Lng的标记表，公式为:

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://*.com/a/45950426/5076414

请阅读我在https://*.com/a/45950426/5076414的详细答案。

#15

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