如何计算给定 lat/lng 位置的边界框？

你要找的是椭球公式。

我发现开始编码的最佳位置是基于 CPAN 的 Geo: : 椭圆体库。它为您提供了一个基线来创建关闭的测试并将结果与其结果进行比较。我以前的雇主将它作为类似 PHP 库的基础。

椭圆体

看一下 location方法，调用两次就得到了 bbox。

You didn't post what language you were using. There may already be a geocoding library available for you.

哦，如果你现在还没弄明白，谷歌地图使用的是 WGS84椭圆体。

我建议将地球表面局部近似为一个半径由 WGS84椭球体在给定纬度给出的球体。我怀疑 latMin 和 latMax 的精确计算将需要椭圆函数，并且不会产生可观的精度提高(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

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

我改编了一个 PHP 脚本来完成这个任务。你可以用它来找到一个盒子的角落周围的一个点(例如，20公里外)。我的具体例子是 Google Maps API:

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

在这里，我把费德里科 · A · 兰波尼的答案转换成了 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));
}

这是非常简单的只是去 panoramio 网站，然后从 panoramio 网站打开世界地图。然后前往指定的经纬度。

然后你在地址栏里发现了经纬度，比如这个地址。

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

纬度32.739485 ln=70.491211 =>longitude

这个 Panoramio JavaScript API 小部件围绕一个 lat/long 对创建一个边界框，然后返回在这个边界内的所有照片。

另一种类型的 PanoramioJavaScriptAPI 小部件，也可以使用 例子和代码在这里更改背景颜色。

它不表现在作曲的情绪上，它表现在出版之后。

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

我正在研究的边界盒问题作为一个侧面的问题，以找到所有的点在 SrcRad 半径的静态 LAT，长点。有相当多的计算使用

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. 我知道这需要一些额外的步骤，但它帮助了我

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

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

@ Jan Philip Matuschek 的插图很好的解释。(请对他的回答投赞成票，而不是这个; 我添加这个，因为我花了一点时间来理解原来的答案)

优化寻找最近邻的边界盒技术需要求出距离 d 上点 P 的最小和最大纬度、经度对。所有落在这些点之外的点，距离都大于 d。 这里需要注意的一点是交点纬度的计算，正如 Jan Philip Matuschek 解释中所强调的那样。交点的纬度不在点 P 的纬度上，而是与点 P 的纬度略有偏移。在确定距离 d 的点 P 的正确最小和最大边界经度时，这是一个经常被忽略的重要部分。这在验证中也很有用。

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

巨蟒的要点在这里 https://gist.github.com/alexcpn/f95ae83a7ee0293a5225

这里我把 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));
}
?>

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

不准确，但很方便。

下面是一个使用 javascript 的简单实现，它基于纬度到公里的转换，其中 1 degree latitude ~ 111.2 km。

我正在根据给定的纬度、经度和半径(以公里为单位)计算地图的边界。

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

感谢@Fedrico A 提供的 Phyton 实现，我已经把它移植到 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. 结构的边界位置应该被一个类替换，我使用它只是为了在这里共享它。

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

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;