如何计算给定纬度/经度位置的边界框?
我给出了一个由经度和纬度定义的位置。 现在我想计算一个例如10公里内的边界框。
边界框应该定义为latmin,lngmin和latmax,lngmax。
我需要这些东西才能使用panoramio API 。
有人知道如何得到积分的公式吗?
编辑:我正在寻找一个公式/函数,它需要经纬度为input,并返回一个边界框作为latmin&lngmin和latmax&latmin。 MySQL,PHP,C#,JavaScript是好的,但也伪代码应该没问题。
编辑:我不是在寻找一个解决scheme,它显示了2点的距离
我build议把地球表面近似看作一个在给定的纬度上由WGS84椭球给出的半径的球体。 我怀疑latMin和latMax的精确计算需要椭圆函数,并且不会产生可观的精度增加(WGS84本身就是一个近似值)。
我的实现如下(它是用Python编写的;我没有testing过):
# 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))
编辑:下面的代码转换(度,素数,秒)度+分数的度,反之亦然(未testing):
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))
我写了一篇关于find边界坐标的文章:
http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates
文章解释了公式,并提供了一个Java实现。 (这也说明了为什么费德里科的最小/最大经度公式是不准确的。)
在这里,我已经将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)); }
我写了一个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() ]; };
你正在寻找一个椭球公式。
我发现开始编码的最好的地方是基于CPAN的Geo :: Ellipsoid库。 它为您提供了一个基线来创build您的testing,并将结果与结果进行比较。 我在之前的雇主中用它作为PHP的类似库的基础。
地球椭球::
看看location
方法。 调用它两次,你有你的Bbox。
您没有发布您正在使用的语言。 可能已经有一个地理编码库供您使用。
哦,如果你现在还没有弄清楚,谷歌地图使用WGS84椭球。
我改编了一个PHP脚本,我发现这样做。 你可以使用它来find一个点周围的一个盒子的angular落(比如20公里外)。 我的具体示例是Google Maps API:
这是一个简单的实现使用JavaScript的基础上,纬度转换为公里,其中1 degree latitude ~ 111.2 km
。
我正在计算从10km宽的给定经度和纬度的地图范围。
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; }
由于我需要一个非常粗略的估计,所以为了在弹性search查询中过滤掉一些不必要的文档,我使用了下面的公式:
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。 对于你的情况N = 10
不准确,但方便。
@Jan Philip Matuschek的一个很好的解释的插图(请把他的答案加起来,而不是这个,因为我花了一点时间来理解原来的答案)
寻找最近邻居最优化的边界框技术将需要为距离d处的点P导出最小和最大纬度,经度对。 所有落在这些之外的点绝对距离这个点大于d。 这里要注意的一点就是Jan Philip Philip Matuschek解释中强调的交叉纬度的计算。 交点的纬度不在P点的纬度,但稍微偏离点。 在确定距离d的点P的正确的最小和最大边界经度时,这是经常被忽略但是重要的部分。这在validation中也是有用的。
(P的纬度,经度高)与(纬度,经度)之间的半正向距离等于距离d。
Python的要点https://gist.github.com/alexcpn/f95ae83a7ee0293a5225
我正在研究边界框问题,作为寻找静态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部分,我离开了一些位置点。
获得所有边界框点后,如果您有一个函数可以计算给定纬度点的点到点距离,则只需要获取与固定点相距一定距离的点即可。 这是我做的。 我知道它采取了一些额外的步骤,但它帮助了我
-- 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 ;
刚刚进入panoramio网站,然后从panoramio网站打开World Map,然后转到指定的经纬度位置。
然后在地址栏中find纬度和经度,例如在这个地址。
http://www.panoramio.com/map#lt=32.739485&ln=70.491211&z=9&k=1&a=1&tab=1&pl=all
lt = 32.739485 =>纬度ln = 70.491211 =>经度
此Panoramio JavaScript API小部件创build一个经纬度/长度对的边界框,然后返回这些边界中的所有照片。
另一种types的Panoramio JavaScript API小部件,您可以在其中使用示例和代码更改背景色。
它没有performance出作曲的心情,出版后显示。
<div dir="ltr" style="text-align: center;" trbidi="on"> <script src="https://ssl.panoramio.com/wapi/wapi.js?v=1&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>
在这里,如果有人感兴趣,我已经将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.对于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
我已经testing了它,似乎工作很好。 Struct BoundsLocation应该被一个类replace,我用它来分享它。
以上所有答案都只是部分正确的 。 特别是在像澳大利亚这样的地区,它们总是包含极点,甚至在10kms的时候也计算出一个很大的矩形。
特别是Jan Philip Matuschek在http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates#UsingIndex上的algorithm包含了一个来自(-37,-90,-180,180)的非常大的矩形,用于几乎澳大利亚的每一个点。; 这对数据库中的大量用户和距离都要计算几乎一半的所有用户。
我发现罗切斯特理工大学的Drupal API Earthalgorithm在极点以及其他地方工作的更好,而且更容易实现。
https://www.rit.edu/drupal/api/drupal/sites%21all%21modules%21location%21earth.inc/7.54
使用上述algorithm中的earth_latitude_range
和earth_longitude_range
来计算边界矩形
并使用谷歌地图logging的距离计算公式来计算距离
对于(Lat,Lng)=(37,-122)和带有lat和lng列的Markers表 ,公式为:
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;
这是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), } }