为啥这个用于 lat/long 的 Lua Haversine 代码不起作用？

Posted 2023-03-05

【中文标题】为啥这个用于 lat/long 的 Lua Haversine 代码不起作用？

【英文标题】：How come this Lua Haversine code for lat/long won't work?为什么这个用于 lat/long 的 Lua Haversine 代码不起作用？

【发布时间】：2020-11-20 05:13:48

【问题描述】：

我正在尝试创建一个半正弦“距离计算器”，我可以在其中输入两个纬度/经度坐标，并通过使用半正弦公式给出它们之间的距离。它在 Lua 中，代码如下：

local R = 6371000 -- metres
local lat1 =  la1 * math.pi/180
local lat2 = la2 * math.pi/180
local dlat = (lat2 - lat1) * math.pi/180
local dlong = (long2 - long1) * math.pi/180

local a = math.sin(dlat/2) * math.sin(dlat/2) + math.cos(lat1) * math.cos(lat2) * math.sin(dlong/2) * math.sin(dlong/2)
local c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))

local d = R * c -- in metres
d = d / 1000 -- in kilometers

print(d)

我已经针对使用此公式的多个在线距离检查器对其进行了测试，结果始终相差很大。而且，我从这里获取了代码，然后将其更改为 Lua：https://www.movable-type.co.uk/scripts/latlong.html

知道为什么这不起作用吗？谢谢。

【问题讨论】：

顺便说一句，您可以使用math.rad 将度数转换为弧度。

【参考方案1】：

假设您只是省略了设置坐标的行，您的问题是您将dlat 乘以math.pi/180，尽管lat1 和lat2 已经是弧度，而不是度数。 local dlat 应该是 lat2 - lat1。

以下是代码的改进变体：

local function haversine (lat1, lat2, long1, long2) -- in radians.
    local cos, sin = math.cos, math.sin
    local dlat, dlong = lat2 - lat1, long2 - long1
    return sin (dlat / 2) ^ 2 + cos (lat1) * cos (lat2) * sin (dlong / 2) ^ 2
end

local function distance (p1, p2) -- in degrees.
    local pi, arcsin, sqrt = math.pi, math.asin, math.sqrt
    local d2r = pi / 180

    local R = 6371000 -- in metres

    local lat1, lat2 =  p1[1] * d2r, p2[1] * d2r
    local long1, long2 = p1[2] * d2r, p2[2] * d2r
    
    local a = haversine (lat1, lat2, long1, long2)

    return 2 * R * arcsin (sqrt (a)) / 1000 -- in km
end


local Moscow, Novokuznetsk = 55.7558, 37.6173, 53.7596, 87.1216
print (distance (Moscow, Novokuznetsk)) -- 3126 km as reported by https://www.distancefromto.net/distance-from-moscow-to-novokuznetsk-ru.

distance 和 haversine 已移至函数中， math 中的函数已本地化，以实现更整洁的外观和性能， 角度到弧度的转换率是针对性能预先计算的， ^是用来供电的， 反正切替换为反正弦，以简化公式。

【讨论】：

谢谢！我还有另一个问题；给定地球上两点的纬度/经度，您会碰巧知道另一种非常准确的计算距离的方法吗？我对文森提公式很熟悉，但也许还有比这更准确的东西？

我可以推荐以下内容：1) 将sin (x / 2) ^ 2 替换为(1 - cos (x)) / 2，将arcsin (sqrt (a)) 替换为arccos (1 - 2 * 2) / 2。这不会提高精度，但会为您节省两个平方和一个sqrt, 2) 通过计算 (0,0) 和 (0,180) 或 (90,0) 和 (-90,0) 之间的距离来校准公式) 并将结果与​​地球圆周的一半进行比较。目前的方法是参考值（40 075 公里）的 0.1%。可以通过将R 替换为地球赤道半径（6378136.6 m，en.wikipedia.org/wiki/Equator#Exact_length）来改进它。

这将给出赤道半径 40075.014172304。任何进一步的改进只能基于考虑地球的真实形状，它不是球体 (en.wikipedia.org/wiki/Geoid)。我不知道有任何简单的公式可以实现这一点。我不熟悉的文森特公式将考虑地球的扁率；任何进一步的改进都必须使用描述地球真实形状的数据集。