Pakiet geopolityczny Pythona zawiera dwie techniki pomiaru odległości: Wielki Krąg i wzory Vincenta .
>>> 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.4088367903001Jaka jest różnica? Który pomiar odległości jest preferowany?
Odpowiedzi:
Według Wikipedii formuła Vincenta jest wolniejsza, ale dokładniejsza :
Wzory Vincenta to dwie powiązane metody iteracyjne stosowane w geodezji do obliczania odległości między dwoma punktami na powierzchni sferoidy, opracowane przez Thaddeusa Vincentego (1975a). Opierają się one na założeniu, że Ziemia jest spłaszczoną sferoidą, a zatem są dokładniejsze niż metody takie jak odległość od wielkiego koła, która zakłada sferyczną Ziemię.
Różnica dokładności ~0.17%w Izraelu wynosi 428 metrów. Zrobiłem szybki i brudny test prędkości:
<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)Kod:
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)Podsumowując: Formuła Vincenty'ego podwaja czas obliczeń w porównaniu do wielkiego koła, a jego przyrost dokładności w testowanym punkcie wynosi ~ 0,17%.
Ponieważ czas obliczeń jest znikomy, formuła Vincenta jest preferowana dla każdej praktycznej potrzeby.
Aktualizacja : Po wnikliwych komentarzach Whucera i cffk i cffk zgadzam się, że przyrost dokładności należy porównać z błędem, a nie z pomiarem. Dlatego wzór Vincenta jest kilka rzędów wielkości bardziej dokładny, a nie ~ 0,17%.
Jeśli korzystasz z geopoli, odległości great_circle i vincenty są równie wygodne do uzyskania. W takim przypadku prawie zawsze powinieneś użyć tego, który daje dokładniejszy wynik, tj. Vincenty. Dwie kwestie (jak wskazałeś) to szybkość i dokładność.
Vincenty jest dwa razy wolniejszy. Ale prawdopodobnie w prawdziwej aplikacji wydłużony czas działania jest znikomy. Nawet jeśli twoja aplikacja wymagała miliona obliczeń odległości, mówimy tylko o różnicy w czasie kilku sekund.
Dla użytych punktów błąd w vincenty wynosi 6 μm, a błąd w odległości dużego koła wynosi 0,75 m. Powiedziałbym wtedy, że vincenty jest 120000 razy bardziej dokładny (zamiast 0,17% bardziej dokładny). W przypadku punktów ogólnych błąd w odległości dużego koła może wynosić nawet 0,5%. Czy możesz żyć z 0,5% błędem w odległości? Do codziennego użytku (jaka jest odległość od Kapsztadu do Kairu?), Prawdopodobnie możesz. Jednak wiele aplikacji GIS ma znacznie bardziej rygorystyczne wymagania dotyczące dokładności. (0,5% to 5 m na 1 km. To naprawdę robi różnicę).
Prawie wszystkie poważne prace związane z mapowaniem przeprowadzane są na elipsoidzie odniesienia, dlatego sensowne jest, aby odległości również mierzyć na elipsoidzie. Może dziś uda ci się uciec z dużymi odległościami. Ale dla każdej nowej aplikacji będziesz musiał sprawdzić, czy jest to nadal dopuszczalne. Lepiej jest użyć elipsoidalnej odległości od samego początku. Będziesz lepiej spać w nocy.
DODATEK (maj 2017 r.)
W odpowiedzi na odpowiedź udzieloną przez @ craig-hicks. Metoda vincenty () w geopy ma potencjalnie śmiertelną wadę: rzuca błąd na punkty prawie antypodalne. Dokumentacja w kodzie sugeruje zwiększenie liczby iteracji. Ale to nie jest ogólne rozwiązanie, ponieważ metoda iteracyjna używana przez vincenty () jest niestabilna dla takich punktów (każda iteracja prowadzi cię dalej od poprawnego rozwiązania).
Dlaczego opisuję problem jako „potencjalnie śmiertelny”? Ponieważ każde użycie funkcji odległości w innej bibliotece oprogramowania musi być w stanie obsłużyć wyjątek. Obsługa go przez zwracanie NaN lub odległości do wielkiego koła może nie być zadowalająca, ponieważ wynikowa funkcja odległości nie będzie przestrzegać nierówności trójkąta, co wyklucza jej użycie, np. W drzewach punktów obserwacyjnych.
Sytuacja nie jest całkowicie ponura. Mój pakiet geograficzny geographiclib dokładnie oblicza odległość geodezyjną bez żadnych awarii. Prośba geopy ciągnąć nr 144 zmienia funkcję oddali geopy do pakietu użycie geographiclib jeśli jest ona dostępna. Niestety to żądanie ściągnięcia jest w zawieszeniu od sierpnia 2016 r.
DODATEK (maj 2018 r.)
geopy 1.13.0 używa teraz pakietu geograficznego do obliczania odległości. Oto przykładowe wywołanie (na podstawie przykładu z pierwotnego pytania):
>>> 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.28877358686776Przepraszam za opublikowanie tutaj drugiej odpowiedzi, ale korzystam z okazji, by odpowiedzieć na prośbę @ craig-hicks, aby zapewnić porównanie dokładności i czasu dla różnych algorytmów obliczania odległości geodezyjnej. Parafrazuje komentarz, który kieruję do mojego żądania ściągnięcia # 144 dotyczącego geopolityki, co pozwala na użycie jednej z dwóch implementacji mojego algorytmu dla geodetyki w geologii, jedna jest natywną implementacją Pythona, geodezyjną (geographiclib) , a inne zastosowania wdrożenie w C geodezyjnej (pyproj) .
Oto niektóre dane dotyczące czasu. Czasy w mikrosekundach na połączenie
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.1Oto dokładność obliczeń geodezyjnych na podstawie mojego zestawu testów geodezyjnych . Błędy podano w jednostkach mikronów (1e-6 m)
method distance destination
geopy vincenty 205.629 141.945
geopy geodesic(pyproj) 0.007 0.013
geopy geodesic(geographiclib) 0.011 0.010Uwzględniłem żądanie ściągnięcia nr 194 od hannosche, które naprawia zły błąd w funkcji docelowej. Bez tej poprawki błąd w obliczeniach miejsca docelowego dla vincenty wynosi 8,98 metra.
19,2% przypadków testowych zakończyło się niepowodzeniem z vincenty.distance (iteracje = 20). Jednak zestaw testowy jest przechylony w kierunku przypadków, które spowodowałyby tę awarię.
W przypadku losowych punktów na elipsoidzie WGS84, algorytm Vincenty'ego z pewnością zawiedzie 16,6 na 1000000 razy (poprawne rozwiązanie jest niestabilnym stałym punktem metody Vincenty'ego).
Przy implementacji geopinii Vincenta i iteracjach = 20 wskaźnik awaryjności wynosi 82,8 na 1000000. Przy iteracjach = 200 wskaźnik awaryjności wynosi 21,2 na 1000000.
Mimo że stawki te są małe, awarie mogą być dość powszechne. Na przykład w zestawie danych 1000 losowych punktów (być może na światowych lotniskach) obliczenie macierzy pełnej odległości zawiodłoby średnio 16 razy (przy iteracjach = 20).
Wygląda na to, że pakiet geopy.distance oferuje funkcję „odległość ()”, która domyślnie jest ustawiona na vincenty (). Zasadniczo poleciłbym używanie funkcji distance (), ponieważ jest to zalecenie dotyczące pakietu, na wypadek, gdyby w przyszłości kiedykolwiek odbiegało od vincenty () (co jest mało prawdopodobne). Kontynuuj czytanie:
Ta uwaga w dokumentacji jest zawarta w kodzie źródłowym określonej przez ciebie funkcji vincenty ():
Uwaga: Ta implementacja odległości Vincenta nie jest zbieżna dla niektórych ważnych punktów. W niektórych przypadkach wynik można uzyskać poprzez zwiększenie liczby iteracji (
iterationsargument słowa kluczowego, podany w klasie__init__z domyślną wartością 20). Preferowane może być użycie: class:,.great_circlektóra jest nieznacznie mniej dokładna, ale zawsze daje wynik.
Kod źródłowy z powyższym komentarzem / notatką można znaleźć na stronie https://github.com/geopy/geopy/blob/master/geopy/distance.py Przewiń w dół do definicji vincenty ()
Niemniej jednak domyślną funkcją odległości używaną przez ten pakiet podczas skalowania odległości () jest funkcja vincenty (), co oznacza, że niepowodzenie zbieżności nie jest katastrofalne, i zwracana jest rozsądna odpowiedź - co najważniejsze, wyjątek nie jest generowany.
Aktualizacja: Jak zauważa „cffk”, funkcja vincenty () jawnie zgłasza wyjątek ValueError, gdy algorytm nie jest zbieżny - chociaż nie jest to udokumentowane w opisie funkcji. Dlatego dokumentacja jest błędna.
Niezależnie od tego, czy używasz vincenty, haverine, czy sferycznego prawa cosinusów, mądrość polega na uświadomieniu sobie wszelkich potencjalnych problemów z kodem, którego planujesz użyć, rzeczy, na które należy uważać i które można złagodzić, oraz tego, jak radzić sobie z problemami vincenty vs. haversine vs. sloc będą się różnić, gdy uświadomimy sobie czyjeś czyhające problemy / obramowania, które mogą, ale nie muszą być powszechnie znane. Doświadczony programista wie o tym. Początkujący mogą nie. Mam nadzieję zaoszczędzić niektórym z nich frustracji, gdy fragment kodu z forum robi coś nieoczekiwanego, w niektórych przypadkach. Jeśli ktoś poważnie zamierza zastosować jakąś wersję któregokolwiek z nich, Vincenty, haverine, sloc, to SE, SO, Reddit, Quora itp., Może zapewnić ograniczoną pomoc w początkowym kodowaniu rozwiązania, ale to nie znaczy, że ich rozwiązanie lub zaakceptowana „odpowiedź” jest wolna od problemów. Jeśli projekt jest wystarczająco ważny, zasługuje na odpowiednią rozsądną ilość badań. Przeczytaj instrukcję, przeczytaj dokumentację, a jeśli istnieje przegląd kodu tego kodu, przeczytaj go. Skopiowanie i wklejenie fragmentu lub treści, która została oceniona sto lub więcej razy, nie oznacza, że jego bezpieczeństwo jest wszechstronne i pewne.
Intrygująca odpowiedź wysłana przez cffk podnosi świadomość, że czają się przypadkowe przypadki w opakowanych rozwiązaniach, które mogą powodować wyjątki lub inne trudności . Konkretne twierdzenia zawarte w tym poście przekraczają obecnie mój budżet czasowy, ale zabieram z tego, że w niektórych pakietach, w tym co najmniej jednej implementacji vincenty, rzeczywiście istnieją problemy, w których co najmniej jedna osoba zaproponowała poprawę w taki czy inny sposób, aby zminimalizować lub wyeliminować ryzyko napotkania tych trudności. Nie będę dodawał więcej do tego tematu dotyczącego Vincenty (będąc zbyt zbyt nieświadomym), ale zamiast tego przejdę do haversine, przynajmniej częściowo na temat z PO.
Popularnie opublikowana formuła haverine, czy to w Pythonie, czy w innym języku, ponieważ najprawdopodobniej będzie korzystała ze specyfikacji zmiennoprzecinkowej IEEE 754 w większości wszystkich systemów Intel i podobnych do Intel, a także procesorów ARM, powerPC itp. być również podatnym na rzadkie, ale rzeczywiste i powtarzalne błędy wyjątków bardzo bliskie lub w odległości 180 stopni łuku, punkty antypodalne, z powodu przybliżeń zmiennoprzecinkowych i zaokrąglania. Niektórzy nowicjusze mogą jeszcze nie zostać ukąszeni przez tę sytuację. Ponieważ ta specyfikacja fp jest przybliżona i zaokrągla, nie oznacza to, że każdy kod wywołujący fp64 może powodować błędy wyjątków, nie. Ale trochę kodu niektóre formuły mogą nie mieć tak oczywistych przypadków, w których aproksymacje i zaokrąglenia IEEE 754 fp64 mogą powodować, że wartość nieznacznie wykracza poza dziedzinę metody matematycznej, która, jak się oczekuje, bezbłędnie oceni taką wartość. Przykład ... sqrt (). Jeśli wartość ujemna znajdzie się w sqrt (), takiej jak sqrt (-0.00000000000000000122739), wystąpi błąd wyjątku. We wzorze haverine, w którym postępuje on w kierunku rozwiązania, istnieją dwie metody sqrt () w atan2 (). Thea, który jest obliczany, a następnie stosowany w sqrt (), może w punktach antypodalnych na kuli ziemskiej nieco zbłądzić poniżej 0,0 lub powyżej 1,0, bardzo nieznacznie z powodu przybliżeń fp64 i zaokrąglania, rzadko, ale powtarzalnie. Konsekwentna niezawodna powtarzalność w tym kontekście sprawia, że jest to wyjątkowe ryzyko, przypadek do ochrony, łagodzenia, a nie izolowany losowy przypadek. Oto przykład krótkiego fragmentu haversine w Pythonie bez niezbędnej ochrony:
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 * cBardzo blisko lub na antypodyczne punktów oblicza się w pierwszym wierszu wzoru może zboczyć ujemny, rzadko, ale z tymi samymi repeatably łac współrzędnych Lon. Aby zabezpieczyć / skorygowania tych rzadkich przypadków, można po prostu dodać, po w obliczeniach, jak widać poniżej:
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.0Oczywiście nie pokazałem tutaj całej funkcji, ale krótki fragment, który jest tak często publikowany. Ale ten pokazuje ochronę sqrt (), testując a i normalizując go, jeśli to konieczne, oszczędzając również potrzeby wypróbowania całości. Uwaga = '' góra do góry ma zapobiegać protestowaniu przez kod bajtowy, że nuta jest używana przed przypisaniem wartości, jeśli zostanie zwrócona z wynikiem funkcji.
Po tej prostej zmianie, dodając dwa testy a , funkcje sqrt () będą zadowolone, a kod ma teraz dodatkową notatkę, którą można zwrócić do kodu wywołującego, aby ostrzec, że wynik został nieco znormalizowany i dlaczego. Niektórzy mogą się tym przejmować, inni nie, ale to tam, zapobiegając błędowi wyjątku, który w przeciwnym razie może wystąpić. Blok try try może przechwycić wyjątek, ale go nie naprawić, chyba że jest to wyraźnie napisane. Wydaje się łatwiej kodzie Correction (s) natychmiast po linią obliczeniowej. Dokładnie wyszorowane wejście nie powinno wtedy wymagać próby, z wyjątkiem bloku tutaj.
Podsumowanie, jeśli używasz haversine, kodowanego jawnie zamiast pakietu lub biblioteki, bez względu na wybrany język, dobrym pomysłem byłoby przetestowanie i znormalizowanie z powrotem do wymaganego zakresu 0,0 <= a <= 1,0 w kolejności aby chronić następny wiersz z jego obliczeniami c . Ale większość fragmentów kodu haverine nie pokazuje tego i nie wspomina o ryzyku.
Doświadczenie: podczas dokładnych testów na całym świecie, w przyrostach co 0,001 stopnia, napełniłem dysk twardy kombinacjami lat lon, które spowodowały wyjątek, niezawodny spójny powtarzalny wyjątek, w ciągu miesiąca równolegle testującego niezawodność chłodzenia procesora fanem i moją cierpliwością. Tak, od tego czasu usunąłem większość tych dzienników, ponieważ ich celem było głównie udowodnienie sensu (jeśli kalambur jest dozwolony). Mam jednak krótsze dzienniki „problematycznych wartości lontu”, przechowywane do celów testowych.
Dokładność: czy a i cały wynik haverine straci pewną dokładność poprzez normalizację tego małego kawałka z powrotem do dziedziny? Niewiele, może nie więcej niż przybliżenia fp64 i zaokrąglenia już były wprowadzane, co spowodowało to nieznaczne przesunięcie poza domenę. Jeśli uznasz, że haversine jest już akceptowalny w stosunku do vincenty - prostsze, szybsze, łatwiejsze do dostosowania, rozwiązywania problemów i konserwacji, haversine może być dobrym rozwiązaniem dla twojego projektu.
Użyłem haversine na rzutowanej sferze nieba do pomiaru odległości kątowych między obiektami na niebie, widzianej z pozycji na ziemi, mapowania azymutu i alt na sferę lat współrzędnych podobnych do współrzędnych, żadnych elipsoid w ogóle do rozważenia, ponieważ rzutowana teoretyczna kula nieba jest idealną kulą, jeśli chodzi o pomiar odległości kątowej kątów patrzenia między dwoma obiektami z pozycji na powierzchni ziemi. Idealnie odpowiada moim potrzebom. Tak więc, haversine jest nadal bardzo przydatny i bardzo dokładny w niektórych aplikacjach (dobrze w moich celach) ... ale jeśli go użyjesz, czy to na ziemi do GIS lub nawigacji, lub w obserwacjach i pomiarach obiektów nieba, chroń to w przypadku antypodyczne punktów lub bardzo blisko antypodyczne punktów, testując Ai w razie potrzeby umieszczając go z powrotem w niezbędnej domenie.
Niechroniony haverine jest dostępny w całym Internecie i widziałem tylko jeden stary post usenet, który pokazał pewną ochronę, myślę, że od kogoś z JPL, i może to być wcześniejsza wersja zmiennoprzecinkowa sprzed 1985 r., IEEE 754. Dwie inne strony wspomniały o możliwych problemach w pobliżu punktów antypodalnych, ale nie opisały tych problemów ani sposobów ich złagodzenia. Pojawiają się więc obawy początkujących (takich jak ja), którzy nie zawsze rozumieją dobre praktyki na tyle dobrze, aby dalej badać i testować edgecases niektórych kodów, które skopiowali i wkleili w zaufany projekt. Intrygujący post cffk był odświeżający, ponieważ był publiczny z tego rodzaju problemami, które nie są często wspominane, rzadko publicznie kodowane dla ochrony we fragmentach i rzadko omawiane w ten sposób, w porównaniu do liczby opublikowanych niechronionych i nieopisanych wersji.
W wersji 20190923 strona wiki dotycząca formuły haverine faktycznie wspomina o problemie możliwym w punktach antypodalnych z powodu problemów z zmiennoprzecinkowymi urządzeniami komputerowymi ... zachęcając ...
https://en.wikipedia.org/wiki/Haversine_formula
(ponieważ ta strona wiki nie ma w tej chwili kotwicy HTML dla sekcji, do której prowadziłbym bezpośredni link, dlatego po załadowaniu strony wyszukaj na tej stronie przeglądarki hasło „Podczas korzystania z tych formuł”, a zobacz problem haversine ze wspomnianymi punktami antypodalnymi, bardziej oficjalnie).
I ta druga strona również bardzo krótko o tym wspomina:
https://www.movable-type.co.uk/scripts/latlong.html
Jeśli ktoś znajdzie na tej stronie wyrażenie „w tym ochrona przed błędami zaokrąglania”, istnieje ...
Jeśli atan2 nie jest dostępny, c można obliczyć z 2 ⋅ asin (min (1, √a)) (w tym ochrona przed błędami zaokrąglania).
Teraz jest rzadki przypadek, w którym wspomniane są błędy zaokrąglania, a ochrona jest pokazywana dla wersji asin (), ale nie jest wymieniona ani pokazana dla wersji atan2 (). Ale wspomniane jest przynajmniej ryzyko błędów zaokrąglania.
Imho, każda aplikacja 24/7/365 korzystająca z haverine, potrzebuje tej ochrony w pobliżu punktów antypodalnych jako ważny i prosty szczegół.
Nie wiem, które pakiety haversine zawierają lub nie obejmują tej ochrony, ale jeśli jesteś nowy w tym wszystkim i zamierzasz używać popularnie opublikowanych wersji „snippet”, teraz wiesz, że potrzebuje ochrony, i ta ochrona jest bardzo łatwa do wdrożenia, to znaczy, jeśli nie używasz vincenty i nie używasz spakowanego haverine bez łatwego dostępu do modyfikacji kodu pakietu.
IOW, niezależnie od tego, czy używasz vincenty, haversine czy sloc, należy zdawać sobie sprawę z wszelkich problemów z kodem, rzeczy, na które należy uważać i je łagodzić, oraz to, jak radzimy sobie z problemami vincenty vs. haversine vs. czające się problemy / edgecases, które mogą, ale nie muszą być powszechnie znane.