如何计算由经纬度指定的两点之间的距离?
为了澄清,我想用千米来表示距离;这些点使用WGS84系统,我想了解可用方法的相对准确性。
如何计算由经纬度指定的两点之间的距离?
为了澄清,我想用千米来表示距离;这些点使用WGS84系统,我想了解可用方法的相对准确性。
当前回答
下面是SQL实现,以km为单位计算距离,
SELECT UserId, ( 3959 * acos( cos( radians( your latitude here ) ) * cos( radians(latitude) ) *
cos( radians(longitude) - radians( your longitude here ) ) + sin( radians( your latitude here ) ) *
sin( radians(latitude) ) ) ) AS distance FROM user HAVING
distance < 5 ORDER BY distance LIMIT 0 , 5;
要获得通过编程语言实现的更多细节,您可以浏览这里给出的php脚本
其他回答
下面是一个c#实现:
static class DistanceAlgorithm
{
const double PIx = 3.141592653589793;
const double RADIUS = 6378.16;
/// <summary>
/// Convert degrees to Radians
/// </summary>
/// <param name="x">Degrees</param>
/// <returns>The equivalent in radians</returns>
public static double Radians(double x)
{
return x * PIx / 180;
}
/// <summary>
/// Calculate the distance between two places.
/// </summary>
/// <param name="lon1"></param>
/// <param name="lat1"></param>
/// <param name="lon2"></param>
/// <param name="lat2"></param>
/// <returns></returns>
public static double DistanceBetweenPlaces(
double lon1,
double lat1,
double lon2,
double lat2)
{
double dlon = Radians(lon2 - lon1);
double dlat = Radians(lat2 - lat1);
double a = (Math.Sin(dlat / 2) * Math.Sin(dlat / 2)) + Math.Cos(Radians(lat1)) * Math.Cos(Radians(lat2)) * (Math.Sin(dlon / 2) * Math.Sin(dlon / 2));
double angle = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
return angle * RADIUS;
}
}
这个链接可能对你有帮助,因为它详细介绍了使用哈弗辛公式来计算距离。
摘录:
这个脚本计算两点之间的大圆距离 也就是说,在地球表面上的最短距离-使用 “半正矢”公式。
function getDistanceFromLatLonInKm(lat1,lon1,lat2,lon2) {
var R = 6371; // Radius of the earth in km
var dLat = deg2rad(lat2-lat1); // deg2rad below
var dLon = deg2rad(lon2-lon1);
var a =
Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.cos(deg2rad(lat1)) * Math.cos(deg2rad(lat2)) *
Math.sin(dLon/2) * Math.sin(dLon/2)
;
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c; // Distance in km
return d;
}
function deg2rad(deg) {
return deg * (Math.PI/180)
}
我在R中做了一个自定义函数,使用R基本包中可用的函数来计算两个空间点之间的距离(km)。
custom_hav_dist <- function(lat1, lon1, lat2, lon2) {
R <- 6371
Radian_factor <- 0.0174533
lat_1 <- (90-lat1)*Radian_factor
lat_2 <- (90-lat2)*Radian_factor
diff_long <-(lon1-lon2)*Radian_factor
distance_in_km <- 6371*acos((cos(lat_1)*cos(lat_2))+
(sin(lat_1)*sin(lat_2)*cos(diff_long)))
rm(lat1, lon1, lat2, lon2)
return(distance_in_km)
}
样例输出
custom_hav_dist(50.31,19.08,54.14,19.39)
[1] 426.3987
PS:要计算以英里为单位的距离,请将函数R(6371)替换为3958.756(海里使用3440.065)。
在提供的代码中有一些错误,我在下面修复了它。
以上所有答案都假定地球是一个球体。然而,更精确的近似是扁球体。
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
精确计算中长点之间距离所需的函数是复杂的,陷阱也很多。我不推荐哈弗辛或其他球形的解决方案,因为有很大的不准确性(地球不是一个完美的球体)。vincenty公式更好,但在某些情况下会抛出错误,即使编码正确。
与其自己编写函数,我建议使用geopy,它已经实现了非常精确的地理库来进行距离计算(论文来自作者)。
#pip install geopy
from geopy.distance import geodesic
NY = [40.71278,-74.00594]
Beijing = [39.90421,116.40739]
print("WGS84: ",geodesic(NY, Beijing).km) #WGS84 is Standard
print("Intl24: ",geodesic(NY, Beijing, ellipsoid='Intl 1924').km) #geopy includes different ellipsoids
print("Custom ellipsoid: ",geodesic(NY, Beijing, ellipsoid=(6377., 6356., 1 / 297.)).km) #custom ellipsoid
#supported ellipsoids:
#model major (km) minor (km) flattening
#'WGS-84': (6378.137, 6356.7523142, 1 / 298.257223563)
#'GRS-80': (6378.137, 6356.7523141, 1 / 298.257222101)
#'Airy (1830)': (6377.563396, 6356.256909, 1 / 299.3249646)
#'Intl 1924': (6378.388, 6356.911946, 1 / 297.0)
#'Clarke (1880)': (6378.249145, 6356.51486955, 1 / 293.465)
#'GRS-67': (6378.1600, 6356.774719, 1 / 298.25)
这个库的唯一缺点是它不支持向量化计算。 对于向量化计算,您可以使用新的gevectorslib。
#pip install geovectorslib
from geovectorslib import inverse
print(inverse(lats1,lons1,lats2,lons2)['s12'])
lat和lon是numpy数组。Geovectorslib是非常准确和非常快!我还没有找到改变椭球的方法。标准采用WGS84椭球,是大多数用途的最佳选择。