我通过简化公式来简化计算。

下面是Ruby版本:

include Math
earth_radius_mi = 3959
radians = lambda { |deg| deg * PI / 180 }
coord_radians = lambda { |c| { :lat => radians[c[:lat]], :lng => radians[c[:lng]] } }

# from/to = { :lat => (latitude_in_degrees), :lng => (longitude_in_degrees) }
def haversine_distance(from, to)
  from, to = coord_radians[from], coord_radians[to]
  cosines_product = cos(to[:lat]) * cos(from[:lat]) * cos(from[:lng] - to[:lng])
  sines_product = sin(to[:lat]) * sin(from[:lat])
  return earth_radius_mi * acos(cosines_product + sines_product)
end

这是我的java实现计算距离经过一些搜索。我用的是世界平均半径(来自维基百科)，单位是千米。İf你想要的结果英里，然后使用世界半径英里。

public static double distanceLatLong2(double lat1, double lng1, double lat2, double lng2) 
{
  double earthRadius = 6371.0d; // KM: use mile here if you want mile result

  double dLat = toRadian(lat2 - lat1);
  double dLng = toRadian(lng2 - lng1);

  double a = Math.pow(Math.sin(dLat/2), 2)  + 
          Math.cos(toRadian(lat1)) * Math.cos(toRadian(lat2)) * 
          Math.pow(Math.sin(dLng/2), 2);

  double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));

  return earthRadius * c; // returns result kilometers
}

public static double toRadian(double degrees) 
{
  return (degrees * Math.PI) / 180.0d;
}

下面是移植到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);
}

下面是Erlang实现

lat_lng({Lat1, Lon1}=_Point1, {Lat2, Lon2}=_Point2) ->
  P = math:pi() / 180,
  R = 6371, % Radius of Earth in KM
  A = 0.5 - math:cos((Lat2 - Lat1) * P) / 2 +
    math:cos(Lat1 * P) * math:cos(Lat2 * P) * (1 - math:cos((Lon2 - Lon1) * P))/2,
  R * 2 * math:asin(math:sqrt(A)).

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

下面是一个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;
    }

}