如何计算由经纬度指定的两点之间的距离?

为了澄清,我想用千米来表示距离;这些点使用WGS84系统,我想了解可用方法的相对准确性。


当前回答

你可以用Haversine公式计算它,它是:

a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)
c = 2 ⋅ atan2( √a, √(1−a) )
d = R ⋅ c

下面给出了一个计算两点之间距离的例子

假设我要计算从新德里到伦敦的距离,那么我该如何使用这个公式:

New delhi co-ordinates= 28.7041° N, 77.1025° E
London co-ordinates= 51.5074° N, 0.1278° W

var R = 6371e3; // metres
var φ1 = 28.7041.toRadians();
var φ2 = 51.5074.toRadians();
var Δφ = (51.5074-28.7041).toRadians();
var Δλ = (0.1278-77.1025).toRadians();

var a = Math.sin(Δφ/2) * Math.sin(Δφ/2) +
        Math.cos(φ1) * Math.cos(φ2) *
        Math.sin(Δλ/2) * Math.sin(Δλ/2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));

var d = R * c; // metres
d = d/1000; // km

其他回答

下面是移植到Java的已接受的答案实现,以备任何人需要。

package com.project529.garage.util;


/**
 * Mean radius.
 */
private static double EARTH_RADIUS = 6371;

/**
 * Returns the distance between two sets of latitudes and longitudes in meters.
 * <p/>
 * Based from the following JavaScript SO answer:
 * http://stackoverflow.com/questions/27928/calculate-distance-between-two-latitude-longitude-points-haversine-formula,
 * which is based on https://en.wikipedia.org/wiki/Haversine_formula (error rate: ~0.55%).
 */
public double getDistanceBetween(double lat1, double lon1, double lat2, double lon2) {
    double dLat = toRadians(lat2 - lat1);
    double dLon = toRadians(lon2 - lon1);

    double a = Math.sin(dLat / 2) * Math.sin(dLat / 2) +
            Math.cos(toRadians(lat1)) * Math.cos(toRadians(lat2)) *
                    Math.sin(dLon / 2) * Math.sin(dLon / 2);
    double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
    double d = EARTH_RADIUS * c;

    return d;
}

public double toRadians(double degrees) {
    return degrees * (Math.PI / 180);
}

我不喜欢添加另一个答案,但谷歌地图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));

关于谷歌的计算有多精确,甚至使用了什么模型都没有任何消息(尽管它说的是“球面”而不是“大地水准面”。顺便说一下,“直线”距离显然不同于一个人在地球表面旅行的距离,而这似乎是每个人都在假设的。

这是一个简单的PHP函数,它将给出一个非常合理的近似值(误差小于+/-1%)。

<?php
function distance($lat1, $lon1, $lat2, $lon2) {

    $pi80 = M_PI / 180;
    $lat1 *= $pi80;
    $lon1 *= $pi80;
    $lat2 *= $pi80;
    $lon2 *= $pi80;

    $r = 6372.797; // mean radius of Earth in km
    $dlat = $lat2 - $lat1;
    $dlon = $lon2 - $lon1;
    $a = sin($dlat / 2) * sin($dlat / 2) + cos($lat1) * cos($lat2) * sin($dlon / 2) * sin($dlon / 2);
    $c = 2 * atan2(sqrt($a), sqrt(1 - $a));
    $km = $r * $c;

    //echo '<br/>'.$km;
    return $km;
}
?>

如前所述;地球不是一个球体。它就像马克·麦奎尔决定用来练习的一个很旧很旧的棒球——到处都是凹痕和凸起。简单的计算(像这样)把它当作一个球体。

不同的方法或多或少的精确取决于你在这个不规则的卵形上的位置以及你的点之间的距离(它们越近,绝对误差范围就越小)。你的期望越精确,计算就越复杂。

更多信息:维基百科地理距离

数学有问题,LUA的学位…如果有人知道修复,请清理这段代码!

与此同时,这里有一个Haversine在LUA中的实现(与Redis一起使用!)

function calcDist(lat1, lon1, lat2, lon2)
    lat1= lat1*0.0174532925
    lat2= lat2*0.0174532925
    lon1= lon1*0.0174532925
    lon2= lon2*0.0174532925

    dlon = lon2-lon1
    dlat = lat2-lat1

    a = math.pow(math.sin(dlat/2),2) + math.cos(lat1) * math.cos(lat2) * math.pow(math.sin(dlon/2),2)
    c = 2 * math.asin(math.sqrt(a))
    dist = 6371 * c      -- multiply by 0.621371 to convert to miles
    return dist
end

干杯!

由于这是关于这个话题最受欢迎的讨论,我将在这里补充我从2019年底到2020年初的经验。为了补充现有的答案-我的重点是找到一个准确和快速(即向量化)的解决方案。

让我们从这里最常用的答案——哈弗辛方法开始。向量化是很简单的,参见下面python中的例子:

def haversine(lat1, lon1, lat2, lon2):
    """
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)

    All args must be of equal length.
    Distances are in meters.
    
    Ref:
    https://stackoverflow.com/questions/29545704/fast-haversine-approximation-python-pandas
    https://ipython.readthedocs.io/en/stable/interactive/magics.html
    """
    Radius = 6.371e6
    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])

    dlon = lon2 - lon1
    dlat = lat2 - lat1

    a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2

    c = 2 * np.arcsin(np.sqrt(a))
    s12 = Radius * c
    
    # initial azimuth in degrees
    y = np.sin(lon2-lon1) * np.cos(lat2)
    x = np.cos(lat1)*np.sin(lat2) - np.sin(lat1)*np.cos(lat2)*np.cos(dlon)
    azi1 = np.arctan2(y, x)*180./math.pi

    return {'s12':s12, 'azi1': azi1}

就精确度而言,它是最不准确的。维基百科在没有任何来源的情况下表示相对偏差平均为0.5%。我的实验显示偏差较小。以下是10万个随机点与我的库的比较,应该精确到毫米级:

np.random.seed(42)
lats1 = np.random.uniform(-90,90,100000)
lons1 = np.random.uniform(-180,180,100000)
lats2 = np.random.uniform(-90,90,100000)
lons2 = np.random.uniform(-180,180,100000)
r1 = inverse(lats1, lons1, lats2, lons2)
r2 = haversine(lats1, lons1, lats2, lons2)
print("Max absolute error: {:4.2f}m".format(np.max(r1['s12']-r2['s12'])))
print("Mean absolute error: {:4.2f}m".format(np.mean(r1['s12']-r2['s12'])))
print("Max relative error: {:4.2f}%".format(np.max((r2['s12']/r1['s12']-1)*100)))
print("Mean relative error: {:4.2f}%".format(np.mean((r2['s12']/r1['s12']-1)*100)))

输出:

Max absolute error: 26671.47m
Mean absolute error: -2499.84m
Max relative error: 0.55%
Mean relative error: -0.02%

因此,在10万对随机坐标上,平均偏差为2.5km,这可能对大多数情况都是好的。

下一个选择是Vincenty公式,精确到毫米,这取决于收敛标准,也可以向量化。它确实有在对跖点附近收敛的问题。你可以通过放宽收敛标准使其收敛于这些点,但准确度会下降到0.25%甚至更多。在对映点之外,Vincenty将提供与地理库相近的结果,相对误差小于1。平均是E-6。

这里提到的Geographiclib实际上是当前的黄金标准。它有几个实现,而且相当快,特别是如果你使用的是c++版本。

Now, if you are planning to use Python for anything above 10k points I'd suggest to consider my vectorized implementation. I created a geovectorslib library with vectorized Vincenty routine for my own needs, which uses Geographiclib as fallback for near antipodal points. Below is the comparison vs Geographiclib for 100k points. As you can see it provides up to 20x improvement for inverse and 100x for direct methods for 100k points and the gap will grow with number of points. Accuracy-wise it will be within 1.e-5 rtol of Georgraphiclib.

Direct method for 100,000 points
94.9 ms ± 25 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
9.79 s ± 1.4 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

Inverse method for 100,000 points
1.5 s ± 504 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
24.2 s ± 3.91 s per loop (mean ± std. dev. of 7 runs, 1 loop each)