在其他答案中，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)
}

我已经创建了这个小Javascript LatLng对象，可能对某人有用。

var latLng1 = new LatLng(5, 3);
var latLng2 = new LatLng(6, 7);
var distance = latLng1.distanceTo(latLng2); 

代码:

/**
 * latLng point
 * @param {Number} lat
 * @param {Number} lng
 * @returns {LatLng}
 * @constructor
 */
function LatLng(lat,lng) {
    this.lat = parseFloat(lat);
    this.lng = parseFloat(lng);

    this.__cache = {};
}

LatLng.prototype = {
    toString: function() {
        return [this.lat, this.lng].join(",");
    },

    /**
     * calculate distance in km to another latLng, with caching
     * @param {LatLng} latLng
     * @returns {Number} distance in km
     */
    distanceTo: function(latLng) {
        var cacheKey = latLng.toString();
        if(cacheKey in this.__cache) {
            return this.__cache[cacheKey];
        }

        // the fastest way to calculate the distance, according to this jsperf test;
        // http://jsperf.com/haversine-salvador/8
        // http://stackoverflow.com/questions/27928
        var deg2rad = 0.017453292519943295; // === Math.PI / 180
        var lat1 = this.lat * deg2rad;
        var lng1 = this.lng * deg2rad;
        var lat2 = latLng.lat * deg2rad;
        var lng2 = latLng.lng * deg2rad;
        var a = (
            (1 - Math.cos(lat2 - lat1)) +
            (1 - Math.cos(lng2 - lng1)) * Math.cos(lat1) * Math.cos(lat2)
            ) / 2;
        var distance = 12742 * Math.asin(Math.sqrt(a)); // Diameter of the earth in km (2 * 6371)

        // cache the distance
        this.__cache[cacheKey] = distance;

        return distance;
    }
};

数学有问题，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

干杯!

精确计算中长点之间距离所需的函数是复杂的，陷阱也很多。我不推荐哈弗辛或其他球形的解决方案，因为有很大的不准确性(地球不是一个完美的球体)。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椭球，是大多数用途的最佳选择。

在我的项目中，我需要计算很多点之间的距离，所以我继续尝试优化我在这里找到的代码。平均而言，在不同的浏览器中，我的新实现的运行速度比获得最多好评的答案快2倍。

function distance(lat1, lon1, lat2, lon2) {
  var p = 0.017453292519943295;    // Math.PI / 180
  var c = Math.cos;
  var a = 0.5 - c((lat2 - lat1) * p)/2 + 
          c(lat1 * p) * c(lat2 * p) * 
          (1 - c((lon2 - lon1) * p))/2;

  return 12742 * Math.asin(Math.sqrt(a)); // 2 * R; R = 6371 km
}

您可以在这里使用我的jsPerf并查看结果。

最近我需要在python中做同样的事情，所以这里是一个python实现:

from math import cos, asin, sqrt, pi

def distance(lat1, lon1, lat2, lon2):
    p = pi/180
    a = 0.5 - cos((lat2-lat1)*p)/2 + cos(lat1*p) * cos(lat2*p) * (1-cos((lon2-lon1)*p))/2
    return 12742 * asin(sqrt(a)) #2*R*asin...

为了完整起见:维基百科上的Haversine。

由于这是关于这个话题最受欢迎的讨论，我将在这里补充我从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)