对于那些寻找基于WGS-84和GRS-80标准的Excel公式的人:

=ACOS(COS(RADIANS(90-Lat1))*COS(RADIANS(90-Lat2))+SIN(RADIANS(90-Lat1))*SIN(RADIANS(90-Lat2))*COS(RADIANS(Long1-Long2)))*6371

源

我已经创建了这个小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;
    }
};

下面是VB的实现。NET，这个实现将根据您传递的Enum值以KM或Miles为单位给您结果。

Public Enum DistanceType
    Miles
    KiloMeters
End Enum

Public Structure Position
    Public Latitude As Double
    Public Longitude As Double
End Structure

Public Class Haversine

    Public Function Distance(Pos1 As Position,
                             Pos2 As Position,
                             DistType As DistanceType) As Double

        Dim R As Double = If((DistType = DistanceType.Miles), 3960, 6371)

        Dim dLat As Double = Me.toRadian(Pos2.Latitude - Pos1.Latitude)

        Dim dLon As Double = Me.toRadian(Pos2.Longitude - Pos1.Longitude)

        Dim a As Double = Math.Sin(dLat / 2) * Math.Sin(dLat / 2) + Math.Cos(Me.toRadian(Pos1.Latitude)) * Math.Cos(Me.toRadian(Pos2.Latitude)) * Math.Sin(dLon / 2) * Math.Sin(dLon / 2)

        Dim c As Double = 2 * Math.Asin(Math.Min(1, Math.Sqrt(a)))

        Dim result As Double = R * c

        Return result

    End Function

    Private Function toRadian(val As Double) As Double

        Return (Math.PI / 180) * val

    End Function

End Class

计算距离——尤其是大距离——的主要挑战之一是解释地球的曲率。如果地球是平的，计算两点之间的距离就会像计算直线一样简单!哈弗辛公式包括一个常数(下面是R变量)，它表示地球的半径。根据你是用英里还是公里来测量，它分别等于3956英里或6367公里。 基本公式是:

Dlon = lon2 - lon1 dat = lat2 - lat1 = (sin (dlat / 2)) ^ 2 + cos (lat1) * cos (lat2) * (sin (dlon / 2)) ^ 2 C = 2 * atan2(√(a)，√(1-a)) distance = R * c(其中R为地球半径) R = 6367公里OR 3956英里

     lat1, lon1: The Latitude and Longitude of point 1 (in decimal degrees)
     lat2, lon2: The Latitude and Longitude of point 2 (in decimal degrees)
     unit: The unit of measurement in which to calculate the results where:
     'M' is statute miles (default)
     'K' is kilometers
     'N' is nautical miles

样本

function distance(lat1, lon1, lat2, lon2, unit) {
    try {
        var radlat1 = Math.PI * lat1 / 180
        var radlat2 = Math.PI * lat2 / 180
        var theta = lon1 - lon2
        var radtheta = Math.PI * theta / 180
        var dist = Math.sin(radlat1) * Math.sin(radlat2) + Math.cos(radlat1) * Math.cos(radlat2) * Math.cos(radtheta);
        dist = Math.acos(dist)
        dist = dist * 180 / Math.PI
        dist = dist * 60 * 1.1515
        if (unit == "K") {
            dist = dist * 1.609344
        }
        if (unit == "N") {
            dist = dist * 0.8684
        }
        return dist
    } catch (err) {
        console.log(err);
    }
}

这是一个简单的javascript函数，从这个链接可能是有用的。不知何故相关，但我们使用谷歌地球javascript插件而不是地图

function getApproximateDistanceUnits(point1, point2) {

    var xs = 0;
    var ys = 0;

    xs = point2.getX() - point1.getX();
    xs = xs * xs;

    ys = point2.getY() - point1.getY();
    ys = ys * ys;

    return Math.sqrt(xs + ys);
}

单位不是距离，而是相对于坐标的比率。还有其他相关的计算，你可以在这里代替getApproximateDistanceUnits函数链接

然后我使用这个函数来查看经纬度是否在半径内

function isMapPlacemarkInRadius(point1, point2, radi) {
    if (point1 && point2) {
        return getApproximateDistanceUnits(point1, point2) <= radi;
    } else {
        return 0;
    }
}

点可以定义为

 $$.getPoint = function(lati, longi) {
        var location = {
            x: 0,
            y: 0,
            getX: function() { return location.x; },
            getY: function() { return location.y; }
        };
        location.x = lati;
        location.y = longi;

        return location;
    };

然后你可以做你的事情，看看一个点是否在一个半径范围内，比如:

 //put it on the map if within the range of a specified radi assuming 100,000,000 units
        var iconpoint = Map.getPoint(pp.latitude, pp.longitude);
        var centerpoint = Map.getPoint(Settings.CenterLatitude, Settings.CenterLongitude);

        //approx ~200 units to show only half of the globe from the default center radius
        if (isMapPlacemarkInRadius(centerpoint, iconpoint, 120)) {
            addPlacemark(pp.latitude, pp.longitude, pp.name);
        }
        else {
            otherSidePlacemarks.push({
                latitude: pp.latitude,
                longitude: pp.longitude,
                name: pp.name
            });

        }

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