我建议在给定纬度下，将地球表面近似为一个半径由WGS84椭球体给出的球体。我怀疑精确计算纬度最小值和纬度最大值需要椭圆函数，并且不会产生可感知的精度增加（WGS84本身就是一个近似）。

我的实现如下（用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))

编辑：以下代码将（度数、分、秒）转换为度数+度的分数，以及相反的转换（未经测试）：

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

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

将此翻译为中文：http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates

这篇文章解释了公式，并提供了Java实现。(它还展示了为什么Federico的最小/最大经度公式是不准确的。)

在这里，我已经将Federico 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));
}

由于我需要一个非常粗略的估算，以便在elasticsearch查询中过滤掉一些不必要的文件，因此我采用了以下公式：

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 = 从给定位置需要的公里数。 对于您的情况，N = 10

不够准确但很实用。

我编写了一个 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()
  ];
};

这是一个简单的javascript实现，它基于将纬度度数转换为公里，其中 1度纬度 ~ 111.2千米。

我正在计算地图的范围，根据给定的纬度、经度和半径（单位为公里）。

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

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

从（纬度交点，高经度）到P的（纬度，经度）的aversine距离等于距离d。

Python代码片段在这里 https://gist.github.com/alexcpn/f95ae83a7ee0293a5225

请把这个翻译成中文： 无法翻译图片，图片显示为英文“enter image description here”。

您正在寻找椭球体公式。

我发现开始编码的最佳地点是基于CPAN的Geo::Ellipsoid库。它为您提供创建测试和将您的结果与其结果进行比较的基准线。我以此为基础在我的前雇主为PHP编写了类似的库。

Geo::Ellipsoid 的中文翻译为“地理椭球体”。

看一下位置方法。调用它两次，你就得到了你的bbox。

你没有发布你正在使用哪种语言。可能已经有一个适用于你的地理编码库。

哦，如果你到現在還沒有搞清楚，Google地圖使用WGS84橢球體。

我适应了一个我发现的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

这是基于经纬度和距离获取边界框坐标的JavaScript代码。已测试并正常工作。

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

使用以下的 getBoundingBox 函数来画出边界框。

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

我正在处理边缘框问题，这是寻找所有位于静态LAT，LONG点的SrcRad半径内的点的副问题。有相当多的计算使用。

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

计算经度范围，但我发现这并不能提供所有需要的答案。因为实际上，你想做的是

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

我知道，我知道答案应该是一样的，但我发现它并不是。我发现，如果我不确定我首先执行（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
;

如果有人感兴趣，这里是我把Federico 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));
}
?>

谢谢@Fedrico A.提供的Python实现。我已经将其移植到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.

这很简单，只需访问Panoramio网站，然后从Panoramio网站中打开世界地图。然后前往所需的指定位置，该位置需要提供纬度和经度。

然后您可以在地址栏中找到纬度和经度，例如在此地址中。

将此翻译成中文：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&k=1&a=1&tab=1&pl=all

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

这个Panoramio JavaScript API小部件会创建一个边界框，围绕着一个经纬度点对，并返回该边界框内的所有照片。

另一种Panoramio JavaScript API小部件，在此小部件中，您还可以使用示例和代码更改背景颜色。

它不在撰写模式中显示，它在发布后显示。

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

所有以上答案只是部分正确的。特别是在像澳大利亚这样的地区，他们总是包括电线杆并计算一个非常大的矩形，即使是10公里。

特别是Jan Philip Matuschek在http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates#UsingIndex上的算法，为澳大利亚的几乎每一个点都包括一个非常大的矩形（-37，-90，-180，180）。这会影响到数据库中的大量用户，并且需要为几乎占据了该国一半的所有用户计算距离。

我发现罗切斯特理工学院的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

使用上述算法中的 earth_latitude_range 和 earth_longitude_range 来计算边界矩形。

并使用谷歌地图记录的距离计算公式来计算距离。

将此翻译为中文：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:

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;

请阅读我的详细回答：https://stackoverflow.com/a/45950426/5076414。

这是Federico Ramponi在Go中的回答。注意：没有错误检查:(

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

这里是 C++ 实现。

#include <array>
#include <cmath>

double degrees_to_radian(double deg)
{
    return deg * M_PI / 180.0;
}

double radian_to_degrees(double rad)
{
    return rad * 180.0 / M_PI ;
}

double WGS84EarthRadius(double lat_radians){
    static const double WGS84_a = 6378137.0;  // Major semiaxis [m]
    static const double WGS84_b = 6356752.3;  // Minor semiaxis [m]
    double An = WGS84_a*WGS84_a * cos(lat_radians);
    double Bn = WGS84_b*WGS84_b * sin(lat_radians);
    double Ad = WGS84_a * cos(lat_radians);
    double Bd = WGS84_b * sin(lat_radians);

    return std::sqrt( (An*An + Bn*Bn)/(Ad*Ad + Bd*Bd) );
}

std::array<double, 4> BoundingBox(double latitude, double longitude, uint32_t range_meters)
{
    double lat = degrees_to_radian(latitude);
    double lon = degrees_to_radian(longitude);
    
    uint32_t halfSide = range_meters;

    double radius =  WGS84EarthRadius(lat);
    double pradius = radius*cos(lat);

    double latMin = lat - halfSide/radius;
    double latMax = lat + halfSide/radius;
    double lonMin = lon - halfSide/pradius;
    double lonMax = lon + halfSide/pradius;

    return {radian_to_degrees(latMin), radian_to_degrees(lonMin), radian_to_degrees(latMax) , radian_to_degrees(lonMax)};
}