可能有一个更简单、更正确的解决方案:地球的周长在赤道上是40000公里，在格林威治(或任何经度)周期上约为37000公里。因此:

pythagoras = function (lat1, lon1, lat2, lon2) {
   function sqr(x) {return x * x;}
   function cosDeg(x) {return Math.cos(x * Math.PI / 180.0);}

   var earthCyclePerimeter = 40000000.0 * cosDeg((lat1 + lat2) / 2.0);
   var dx = (lon1 - lon2) * earthCyclePerimeter / 360.0;
   var dy = 37000000.0 * (lat1 - lat2) / 360.0;

   return Math.sqrt(sqr(dx) + sqr(dy));
};

我同意它应该被微调，我自己说过它是一个椭球，所以半径乘以余弦值是不同的。但它更准确一点。与谷歌map相比，误差明显减小。

我在这里发布了我的工作示例。

在MySQL中列出表中指定点(我们使用一个随机点- lat:45.20327, long:23.7806)之间距离小于50 KM的所有点(表中字段为coord_lat和coord_long):

列出所有距离<50，单位:公里(地球半径6371公里):

SELECT denumire, (6371 * acos( cos( radians(45.20327) ) * cos( radians( coord_lat ) ) * cos( radians( 23.7806 ) - radians(coord_long) ) + sin( radians(45.20327) ) * sin( radians(coord_lat) ) )) AS distanta 
FROM obiective 
WHERE coord_lat<>'' 
    AND coord_long<>'' 
HAVING distanta<50 
ORDER BY distanta desc

上面的例子是在MySQL 5.0.95和5.5.16 (Linux)中测试的。

function distance($lat1, $lon1, $lat2, $lon2) { 
    $pi80 = M_PI / 180; 
    $lat1 *= $pi80; $lon1 *= $pi80; $lat2 *= $pi80; $lon2 *= $pi80; 
    $dlat = $lat2 - $lat1; 
    $dlon = $lon2 - $lon1; 
    $a = sin($dlat / 2) * sin($dlat / 2) + cos($lat1) * cos($lat2) * sin($dlon / 2) * sin($dlon / 2);  
    $km = 6372.797 * 2 * atan2(sqrt($a), sqrt(1 - $a)); 
    return $km; 
}

这个链接可能对你有帮助，因为它详细介绍了使用哈弗辛公式来计算距离。

摘录:

这个脚本计算两点之间的大圆距离 也就是说，在地球表面上的最短距离-使用 “半正矢”公式。

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)
}

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

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

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