Почему этот код Lua Haversine для lat/long не работает?

Я пытаюсь создать калькулятор расстояния гаверсинуса, где я могу ввести две координаты широты/долготы, и он даст расстояние между ними с помощью формулы гаверсинуса. Это на 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

Любая идея, почему это не работает? Спасибо.


debugging latitude-longitude math lua haversine
person Matthew T.    schedule 20.11.2020    source источник
comment
Кстати, вы можете использовать math.rad для преобразования градусов в радианы.   -  person lhf    schedule 20.11.2020

Ответы (1)


arrow_upward
1
arrow_downward

Ваша проблема, если предположить, что вы просто пропустили строки, в которых вы задали координаты, заключается в том, что вы умножаете 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 локализованы для более аккуратного вида и производительности,
  • коэффициент преобразования градусов в радианы предварительно рассчитывается для производительности,
  • ^ используется для питания,
  • арктангенс заменен арксинусом, чтобы упростить формулу.
person Alexander Mashin    schedule 20.11.2020
comment
Благодарю вас! У меня есть еще один вопрос; Вы случайно не знаете другой очень точный способ расчета расстояний с учетом широты и долготы двух точек на Земле? Я знаком с формулой Винсенти, но, может быть, есть что-то еще более точное? - person Matthew T.; 21.11.2020
comment
Могу порекомендовать следующее: 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) и сравнивая результат с половиной окружности Земли. Текущий метод отличается от эталонного значения на 0,1% (40 075 км). Его можно улучшить, заменив R экваториальным радиусом Земли (6378136,6 м, en.wikipedia.org/ wiki/Equator#Exact_length). - person Alexander Mashin; 21.11.2020
comment
Что даст экваториальный радиус 40075,014172304. Любые дальнейшие улучшения могут быть основаны только на учете реальной формы Земли, которая не является сферой (en.wikipedia .org/wiki/Geoid). Я не знаю никаких простых формул для достижения этого. Формулы Винсенти, с которыми я не знаком, будут учитывать сплюснутость Земли; и любое дальнейшее улучшение должно осуществляться с использованием набора данных, описывающего реальную форму Земли. - person Alexander Mashin; 21.11.2020