Разница между вычислениями Винсенти и большого круга?

Пакет geopy для Python содержит две методики измерения расстояний: формулы Большого круга и Винсенти .

>>> from geopy.distance import great_circle
>>> from geopy.distance import vincenty
>>> p1 = (31.8300167,35.0662833) # (lat, lon) - https://goo.gl/maps/TQwDd
>>> p2 = (31.8300000,35.0708167) # (lat, lon) - https://goo.gl/maps/lHrrg
>>> vincenty(p1, p2).meters
429.16765838976664
>>> great_circle(p3, p4).meters
428.4088367903001

В чем разница? Какое измерение расстояния является предпочтительным?

Согласно Википедии, формула Винсенти медленнее, но точнее :

Формулы Винсенти - это два взаимосвязанных итерационных метода, используемых в геодезии для расчета расстояния между двумя точками на поверхности сфероида, разработанного Таддеем Винсенти (1975a). Они основаны на предположении, что фигура Земли представляет собой сжатый сфероид, и, следовательно, более точны, чем методы, такие как расстояние по большому кругу, которые предполагают сферическую Землю.

Разница ~0.17%в точности составляет 428 метров в Израиле. Я сделал быстрый и быстрый тест скорости:

<class 'geopy.distance.vincenty'>       : Total 0:00:04.125913, (0:00:00.000041 per calculation)
<class 'geopy.distance.great_circle'>   : Total 0:00:02.467479, (0:00:00.000024 per calculation)

Код:

import datetime
from geopy.distance import great_circle
from geopy.distance import vincenty
p1 = (31.8300167,35.0662833)
p2 = (31.83,35.0708167)

NUM_TESTS = 100000
for strategy in vincenty, great_circle:
    before = datetime.datetime.now()
    for i in range(NUM_TESTS):
        d=strategy(p1, p2).meters
    after = datetime.datetime.now()
    duration = after-before
    print "%-40s: Total %s, (%s per calculation)" % (strategy, duration, duration/NUM_TESTS)

В заключение: формула Винсенти удваивает время расчета по сравнению с большим кругом, а ее точность в проверенной точке составляет ~ 0,17%.

Поскольку время расчета незначительно, формула Винсенти предпочтительна для любых практических нужд.

Обновление : После проницательных замечаний whuber и cffk «ами и cffk » ы ответа, я согласен , что выигрыш точности должен быть сопоставлен с ошибкой, а не измерением. Следовательно, формула Винсенти на несколько порядков точнее, а не ~ 0,17%.

Если вы используете geopy, то расстояния great_circle и vincenty одинаково удобно получить. В этом случае вы почти всегда должны использовать тот, который дает вам более точный результат, т.е. Винсент. Двумя соображениями (как вы указываете) являются скорость и точность.

Винсенти в два раза медленнее. Но, вероятно, в реальном приложении увеличенное время выполнения незначительно. Даже если ваше приложение требовало расчета миллиона расстояний, мы говорим только о разнице в пару секунд.

Для точек, которые вы используете, ошибка в Винсенте составляет 6 мкм, а ошибка на расстоянии большого круга - 0,75 м. Тогда я бы сказал, что Винсенти в 120000 раз точнее (а не на 0,17% точнее). Для общих точек ошибка на расстоянии большого круга может составлять до 0,5%. Так можно ли жить с ошибкой 0,5% на расстоянии? Для повседневного использования (какое расстояние от Кейптауна до Каира?), Вероятно, вы можете. Однако многие приложения ГИС предъявляют более строгие требования к точности. (0,5% - это 5 м на 1 км. Это действительно имеет значение.)

Почти все серьезные картографические работы выполняются на эталонном эллипсоиде, и поэтому имеет смысл измерять расстояния и на эллипсоиде. Может быть, вы можете сойти с дистанции большого круга сегодня. Но для каждого нового приложения вам нужно будет проверить, является ли это приемлемым. Лучше всего использовать эллипсоидальное расстояние от начала. Вы будете лучше спать по ночам.

ADDENDUM (май 2017 г.)

В ответ на ответ, данный @ craig-hicks. У метода vincenty () в geopy есть потенциально фатальный недостаток: он выдает ошибку для почти антиподальных точек. Документация в коде предлагает увеличить количество итераций. Но это не общее решение, потому что итерационный метод, используемый vincenty (), нестабилен для таких точек (каждая итерация уводит вас от правильного решения).

Почему я характеризую проблему как «потенциально смертельную»? Потому что любое использование функции расстояния в другой библиотеке программного обеспечения должно быть в состоянии обработать исключение. Обработка этого с помощью возврата NaN или расстояния большого круга может быть неудовлетворительной, потому что результирующая функция расстояния не будет подчиняться неравенству треугольника, которое исключает его использование, например, в деревьях точек обзора.

Ситуация не совсем безрадостная. Мой пакет Python geographiclib точно вычисляет геодезическое расстояние без каких-либо сбоев. Запрос geopy тянуть # 144 изменяет расстояние функцию geopy, чтобы пакет использование geographiclib , если он доступен. К сожалению, этот запрос на получение доступа находится в подвешенном состоянии с августа 2016 года.

ДОБАВЛЕНИЕ (май 2018 г.)

geopy 1.13.0 теперь использует пакет geographiclib для вычисления расстояний. Вот пример вызова (на основе примера из исходного вопроса):

>>> from geopy.distance import great_circle
>>> from geopy.distance import geodesic
>>> p1 = (31.8300167,35.0662833) # (lat, lon) - https://goo.gl/maps/TQwDd
>>> p2 = (31.8300000,35.0708167) # (lat, lon) - https://goo.gl/maps/lHrrg
>>> geodesic(p1, p2).meters
429.1676644986777
>>> great_circle(p1, p2).meters
428.28877358686776

Приношу свои извинения за публикацию второго ответа здесь, но я воспользовался возможностью, чтобы ответить на запрос @ craig-hicks, чтобы обеспечить сравнение точности и сроков для различных алгоритмов для вычисления геодезического расстояния. Это перефразирует комментарий, который я делаю к моему запросу # 144 для геопии, который позволяет использовать одну из двух реализаций моего алгоритма для геодезических, которые будут использоваться в геопии, одна - это собственная реализация Python, геодезическая (geographiclib) , а другая использует реализация в C, геодезический (pyproj) .

Вот некоторые временные данные. Время в микросекундах за звонок

method                          dist    dest
geopy great_circle              20.4    17.1
geopy vincenty                  40.3    30.4
geopy geodesic(pyproj)          37.1    31.1
geopy geodesic(geographiclib)  302.9   124.1

Вот точность геодезических расчетов, основанных на моем наборе геодезических испытаний . Ошибки приведены в микронах (1е-6 м)

method                        distance destination
geopy vincenty                 205.629  141.945
geopy geodesic(pyproj)           0.007    0.013
geopy geodesic(geographiclib)    0.011    0.010

Я включил в запрос Hannosche № 194, который исправляет ошибку в функции назначения. Без этого исправления ошибка в расчете назначения для Винсента составляет 8,98 метра.

19,2% тестовых случаев не удалось с vincenty.distance (итераций = 20). Тем не менее, набор тестов ориентирован на случаи, которые могут привести к этому отказу.

При наличии случайных точек на эллипсоиде WGS84 алгоритм Винсенти гарантированно потерпит неудачу 16,6 из 1000000 раз (правильное решение - неустойчивая фиксированная точка метода Винсенти).

При геополитической реализации Vincenty и итераций = 20 частота отказов составляет 82,8 на 1000000. При итерациях = 200 частота отказов составляет 21,2 на 1000000.

Даже если эти показатели невелики, сбои могут быть довольно распространенными. Например, в наборе данных из 1000 случайных точек (возможно, подумайте об аэропортах мира), вычисление матрицы полного расстояния в среднем приведет к сбою в 16 раз (с итерациями = 20).

Похоже, что пакет geopy.distance предлагает функцию «distance ()», которая по умолчанию имеет значение vincenty (). Я бы порекомендовал использовать distance () в принципе, так как это пакетная рекомендация, на случай, если когда-нибудь будет отклонен от vincenty () в будущем (вряд ли это так). Продолжить чтение:

Это примечание к документации включено в исходный код указанной вами функции vincenty ():

Примечание: эта реализация расстояния Винсенти не сходится для некоторых действительных точек. В некоторых случаях результат может быть получен путем увеличения числа итераций ( iterationsключевое слово аргумента, заданное в классе __init__, по умолчанию 20). Может быть предпочтительнее использовать: class:, .great_circleкоторый немного менее точен, но всегда дает результат.

Исходный код с указанным выше комментарием / заметкой можно найти по адресу https://github.com/geopy/geopy/blob/master/geopy/distance.py. Прокрутите вниз до определения vincenty ().

Тем не менее, функцией расстояния по умолчанию, используемой этим пакетом при вызове distance (), является функция vincenty (), которая подразумевает, что сбой не является катастрофическим, и возвращается разумный ответ - самое важное, исключение не генерируется.

Обновление: как отмечено в "cffk", функция vincenty () явно генерирует исключение ValueError, когда алгоритм не сходится - хотя это не описано в описании функции. Следовательно, документация глючит.

Независимо от того, используете ли вы Винсенти, Гаверсин или сферический закон косинусов, есть мудрость в осознании любых потенциальных проблем с кодом, который вы планируете использовать, о вещах, которые нужно остерегаться и смягчать, а также о том, как бороться с проблемами Винсента и Хаверсайна против Слока. будет отличаться по мере того, как кто-то узнает о скрывающихся проблемах / крайностях каждого из них, которые могут быть или не быть общеизвестными. Опытный программист знает это. Новички не могут. Я надеюсь избавить некоторых из них от разочарования, когда фрагмент из форума делает что-то неожиданное, в некоторых случаях. Если кто-то серьезно собирается использовать какую-либо из этих версий, vincenty, haversine, sloc, то SE, SO, Reddit, Quora и т. Д., Возможно, оказали ограниченную помощь в некотором начальном кодировании решения, но это не означает, что их решение или принятый «ответ» свободны от вопросов. Если проект достаточно важен, он заслуживает соответствующего разумного количества исследований. Прочитайте руководство, прочитайте документы, и если обзор кода этого кода существует, прочитайте это. Копирование и вставка фрагмента или текста, за который проголосовали сто и более раз, не означает, что его безопасность является всеобъемлющей и гарантированной.

Интригующий ответ, опубликованный cffk, поднимает вопрос о том, как скрываться в пограничных случаях в пакетных решениях, которые могут вызвать исключения или другие трудности . Конкретные претензии, выдвинутые на этом посту, выходят за рамки моего временного бюджета, который я могу реализовать в настоящее время, но я забираю из этого, что в некоторых пакетах действительно есть скрывающиеся проблемы, в том числе, по крайней мере, одна реализация Vincenty, в отношении которых, по крайней мере, один человек предложил улучшить так или иначе, чтобы минимизировать или устранить риск возникновения этих трудностей. Я не буду вдаваться в подробности этой темы, касающейся Винсентства (будучи слишком невежественным в этом), но вместо этого обратимся к хаверсину, по крайней мере, частично по теме с ОП.

Популярная формула haversine, будь то на python или другом языке, потому что она, скорее всего, будет использовать спецификацию IEEE 754 с плавающей запятой на большинстве современных систем intel и intel-like, а также на процессорах ARM, powerPC и т. Д. также быть восприимчивым к редким, но реальным и повторяющимся ошибкам исключения очень близко или на расстоянии дуги 180 градусов, антиподальных точках, из-за приближений с плавающей точкой и округления. Некоторые новички, возможно, еще не были укушены этой ситуацией. Поскольку эта спецификация fp аппроксимирует и округляет, это не означает, что любой код, вызывающий fp64, может вызвать ошибки исключений, нет. Но какой-то код, некоторые формулы могут иметь не столь очевидные случаи краев, где аппроксимации и округления IEEE 754 fp64 могут привести к тому, что значение будет немного отклоняться от области математического метода, который, как ожидается, безошибочно оценивает такое значение. Пример ... sqrt (). Если отрицательное значение попадает в sqrt (), например sqrt (-0.00000000000000000122739), возникнет ошибка исключения. В формуле haversine, по которой она движется к решению, в atan2 () есть два метода sqrt (). a, который вычисляется и затем используется в sqrt (), может в антиподальных точках земного шара слегка отклоняться ниже 0,0 или выше 1,0, очень незначительно из-за приближений fp64 и округления, редко, но с повторением. Постоянная надежная повторяемость, в этом контексте, делает это исключительным риском, крайним случаем для защиты, чтобы смягчить, а не изолированный случайный случайный случай. Вот пример короткого фрагмента haversine для python3 без необходимой защиты:

import math as m

a = m.sin(dlat / 2)**2 + m.cos(lat1) * m.cos(lat2) * m.sin(dlon / 2)**2
c = 2 * m.atan2(m.sqrt(a), m.sqrt(1 - a))
distance = Radius * c

Очень рядом или в диаметрально противоположных точек, рассчитывается в первой строке формулы могут отклониться отрицательные, редко, но повторяемостью с теми же координатами широты Lon. Для защиты / исправить те редкие случаи, можно просто добавить, после в расчет, как показано ниже:

import math as m

note = ''

a = m.sin(dlat / 2)**2 + m.cos(lat1) * m.cos(lat2) * m.sin(dlon / 2)**2
if a < 0.0: a = 0.0 ; note = '*'
if a > 1.0: a = 1.0 ; note = '**'
c = 2 * m.atan2(m.sqrt(a), m.sqrt(1 - a))
distance = Radius * c

# note = '*'  # a went below 0.0 and was normalized back to 0.0
# note = '**' # a went above 1.0 and was normalized back to max of 1.0

Конечно, я не показывал всю функцию здесь, но короткий фрагмент, как это часто публикуется. Но этот показывает защиту для sqrt (), проверяя a и нормализуя его при необходимости, также избавляя от необходимости ставить все целиком, кроме. Вверху примечание = '' предотвращает протестирование на стадии байт-кода того, что примечание используется перед присвоением значения, если оно возвращается с результатом функции.

С этим простым изменением, добавив два теста a , функции sqrt () будут довольны, и теперь в коде есть дополнительная заметка, которая может быть возвращена вызывающему коду, чтобы предупредить, что результат был слегка нормализован, и почему. Некоторых это может волновать, а некоторых - нет, но его, предотвращая ошибку исключения, которая «может» возникнуть в противном случае. Попытка кроме блока может перехватить исключение, но не исправить его, если это явно не написано. Кажется , проще код коррекции линии (ы) сразу после расчета линии. Тщательно очищенный ввод не должен требовать попытки, кроме блока здесь.

Резюме, при использовании гаверсинуса, кодируются в явном виде , а не не с помощью пакета или библиотеки, независимо от того , ваш язык по выбору, было бы хорошей идеей для проверки и нормализации в заднюю часть в потребу диапазоне 0,0 <= а <= 1,0 в порядке защитить следующую строку с ее c расчетами. Но большинство фрагментов кода haversine не показывают это и не упоминают о риске.

Опыт работы: во время тщательного тестирования по всему миру, с шагом 0,001 градуса, я заполнил жесткий диск комбинациями латонов, которые вызвали исключение, надежное постоянное повторяемое исключение, в течение месяца также проводилось совместное тестирование надежности охлаждения процессора. фанат, и мое терпение. Да, с тех пор я удалил большинство из этих журналов, поскольку их целью было в основном доказать свою точку зрения (если каламбур разрешен). Но у меня есть несколько более коротких журналов «значений проблемных длин», которые хранятся в целях тестирования.

Точность. Потеряет ли a и весь результат haversine некоторую точность, нормализуя его обратно в область? Немного, может быть, не больше, чем уже вводили приближения и округления fp64, которые вызвали этот небольшой дрейф вне области. Если вы нашли haversine приемлемым по сравнению с vincenty - проще, быстрее, проще в настройке, устранении неполадок и обслуживании, то haversine может стать хорошим решением для вашего проекта.

Я использовал haversine для проецируемой над головой небесной сферы для измерения угловых расстояний между объектами в небе, если смотреть с позиции на земле, отображая азимут и альт в небосферу по широтно-эквивалентным координатам, не рассматривая элипсоид вообще, так как Прогнозируемая теоретическая небесная сфера является идеальной сферой, когда речь идет об измерении угловых углов взгляда расстояния между двумя объектами из положения на поверхности земли. Это соответствует моим потребностям идеально. Таким образом, haversine все еще очень полезен и очень точен в определенных приложениях (в моих целях) ... но если вы его используете, будь то на Земле для ГИС или навигации, или в наблюдениях и измерениях небесных объектов, защитите это в случае диаметрально противоположных точек или очень близко диаметрально противоположных точек, путем тестирования Aи подталкивать его обратно в нужную область, когда это необходимо.

Незащищенный хаверсин есть по всему Интернету, и я видел только одну старую статью по usenet, которая показала некоторую защиту, я думаю, от кого-то из JPL, и это, возможно, было до 1985 года, до IEEE 754 спецификации с плавающей запятой. На двух других страницах упоминались возможные проблемы рядом с антиподальными точками, но не описывались эти проблемы или как их можно смягчить. Таким образом, существует обеспокоенность для новичков (таких как я), которые не всегда могут достаточно хорошо понимать передовую практику для дальнейшего исследования и тестирования граничных сценариев некоторого кода, который они скопировали и вставили в доверительный проект. Интригующий пост cffk был освежающим в том смысле, что он был публичным с этими типами проблем, которые не часто упоминаются, редко публично кодируются для защиты во фрагментах и ​​редко обсуждаются таким образом, по сравнению с количеством незащищенных и не обсуждаемых версий, которые публикуются.

По состоянию на 20190923 г. на вики-странице формулы haversine действительно упоминается проблема, возможная в антиподальных точках из-за проблем с плавающей запятой в вычислительных устройствах ... обнадеживающих ...

https://en.wikipedia.org/wiki/Haversine_formula

(потому что на этой вики-странице в настоящее время нет html-якоря для раздела, на который я бы непосредственно ссылался, поэтому после загрузки страницы выполните поиск на этой странице браузера для «При использовании этих формул», и вы см. проблему haversine с упомянутыми пунктами antipodal, более официально.)

И у этого другого сайта также есть очень краткое упоминание об этом:

https://www.movable-type.co.uk/scripts/latlong.html

Если кто-то найдет на этой странице слово «включая защиту от ошибок округления», это будет ...

Если atan2 недоступен, c можно рассчитать из 2 ⋅ asin (мин (1, √a)) (включая защиту от ошибок округления).

Теперь есть редкий случай, когда упоминаются ошибки округления, и защита показана для версии asin (), но не упоминается и не показывается для версии atan2 (). Но по крайней мере риск ошибки округления упомянут.

imho, любому приложению 24/7/365, использующему haversine, нужна эта защита вблизи антиподальных точек как важная и простая деталь.

Я не знаю, какие пакеты haversine включают или не включают в себя эту защиту, но если вы новичок во всем этом и собираетесь использовать популярные версии «сниппетов», теперь вы знаете, что он нуждается в защите, и эта защита очень проста в реализации, то есть, если вы не используете vincenty и не используете упакованную haversine без легкого доступа для изменения кода пакета.

В любом случае, будь то Винсенти, Хаверсин или Слок, нужно знать о любых проблемах с кодом, о вещах, которые нужно остерегаться и смягчать, а также о том, как решать проблемы Винсента и Хаверсайна и Слока, будет по-разному, так как каждый станет осознавать проблемы каждого из них. скрывающиеся проблемы / крайние случаи, которые могут быть или не быть общеизвестными.