如何计算由经纬度指定的两点之间的距离?
为了澄清,我想用千米来表示距离;这些点使用WGS84系统,我想了解可用方法的相对准确性。
如何计算由经纬度指定的两点之间的距离?
为了澄清,我想用千米来表示距离;这些点使用WGS84系统,我想了解可用方法的相对准确性。
当前回答
在提供的代码中有一些错误,我在下面修复了它。
以上所有答案都假定地球是一个球体。然而,更精确的近似是扁球体。
a= 6378.137#equitorial radius in km
b= 6356.752#polar radius in km
def Distance(lat1, lons1, lat2, lons2):
lat1=math.radians(lat1)
lons1=math.radians(lons1)
R1=(((((a**2)*math.cos(lat1))**2)+(((b**2)*math.sin(lat1))**2))/((a*math.cos(lat1))**2+(b*math.sin(lat1))**2))**0.5 #radius of earth at lat1
x1=R1*math.cos(lat1)*math.cos(lons1)
y1=R1*math.cos(lat1)*math.sin(lons1)
z1=R1*math.sin(lat1)
lat2=math.radians(lat2)
lons2=math.radians(lons2)
R2=(((((a**2)*math.cos(lat2))**2)+(((b**2)*math.sin(lat2))**2))/((a*math.cos(lat2))**2+(b*math.sin(lat2))**2))**0.5 #radius of earth at lat2
x2=R2*math.cos(lat2)*math.cos(lons2)
y2=R2*math.cos(lat2)*math.sin(lons2)
z2=R2*math.sin(lat2)
return ((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)**0.5
其他回答
在其他答案中,r中的实现是缺失的。
用地质圈包中的distm函数计算两点之间的距离非常简单:
distm(p1, p2, fun = distHaversine)
地点:
p1 = longitude/latitude for point(s)
p2 = longitude/latitude for point(s)
# type of distance calculation
fun = distCosine / distHaversine / distVincentySphere / distVincentyEllipsoid
由于地球不是完美的球形,所以椭球体的文森提公式可能是计算距离的最佳方法。因此,在地质圈包中,您可以使用:
distm(p1, p2, fun = distVincentyEllipsoid)
当然,你不一定要使用geosphere包,你也可以用一个函数来计算以R为基底的距离:
hav.dist <- function(long1, lat1, long2, lat2) {
R <- 6371
diff.long <- (long2 - long1)
diff.lat <- (lat2 - lat1)
a <- sin(diff.lat/2)^2 + cos(lat1) * cos(lat2) * sin(diff.long/2)^2
b <- 2 * asin(pmin(1, sqrt(a)))
d = R * b
return(d)
}
我在这里发布了我的工作示例。
在MySQL中列出表中指定点(我们使用一个随机点- lat:45.20327, long:23.7806)之间距离小于50 KM的所有点(表中字段为coord_lat和coord_long):
列出所有距离<50,单位:公里(地球半径6371公里):
SELECT denumire, (6371 * acos( cos( radians(45.20327) ) * cos( radians( coord_lat ) ) * cos( radians( 23.7806 ) - radians(coord_long) ) + sin( radians(45.20327) ) * sin( radians(coord_lat) ) )) AS distanta
FROM obiective
WHERE coord_lat<>''
AND coord_long<>''
HAVING distanta<50
ORDER BY distanta desc
上面的例子是在MySQL 5.0.95和5.5.16 (Linux)中测试的。
我不喜欢添加另一个答案,但谷歌地图API v.3具有球形几何(以及更多)。在将你的WGS84转换为十进制度后,你可以这样做:
<script src="http://maps.google.com/maps/api/js?sensor=false&libraries=geometry" type="text/javascript"></script>
distance = google.maps.geometry.spherical.computeDistanceBetween(
new google.maps.LatLng(fromLat, fromLng),
new google.maps.LatLng(toLat, toLng));
关于谷歌的计算有多精确,甚至使用了什么模型都没有任何消息(尽管它说的是“球面”而不是“大地水准面”。顺便说一下,“直线”距离显然不同于一个人在地球表面旅行的距离,而这似乎是每个人都在假设的。
FSharp版本,使用里程:
let radialDistanceHaversine location1 location2 : float =
let degreeToRadian degrees = degrees * System.Math.PI / 180.0
let earthRadius = 3959.0
let deltaLat = location2.Latitude - location1.Latitude |> degreeToRadian
let deltaLong = location2.Longitude - location1.Longitude |> degreeToRadian
let a =
(deltaLat / 2.0 |> sin) ** 2.0
+ (location1.Latitude |> degreeToRadian |> cos)
* (location2.Latitude |> degreeToRadian |> cos)
* (deltaLong / 2.0 |> sin) ** 2.0
atan2 (a |> sqrt) (1.0 - a |> sqrt)
* 2.0
* earthRadius
可能有一个更简单、更正确的解决方案:地球的周长在赤道上是40000公里,在格林威治(或任何经度)周期上约为37000公里。因此:
pythagoras = function (lat1, lon1, lat2, lon2) {
function sqr(x) {return x * x;}
function cosDeg(x) {return Math.cos(x * Math.PI / 180.0);}
var earthCyclePerimeter = 40000000.0 * cosDeg((lat1 + lat2) / 2.0);
var dx = (lon1 - lon2) * earthCyclePerimeter / 360.0;
var dy = 37000000.0 * (lat1 - lat2) / 360.0;
return Math.sqrt(sqr(dx) + sqr(dy));
};
我同意它应该被微调,我自己说过它是一个椭球,所以半径乘以余弦值是不同的。但它更准确一点。与谷歌map相比,误差明显减小。