Perbedaan antara Vincenty dan perhitungan jarak lingkaran-besar?

Paket geopy Python menampilkan dua teknik pengukuran jarak: Great Circle dan formula Vincenty .

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

Apa bedanya? Pengukuran jarak mana yang lebih disukai?

Menurut Wikipedia, rumus Vincenty lebih lambat tetapi lebih akurat :

Rumus Vincenty adalah dua metode iteratif terkait yang digunakan dalam geodesi untuk menghitung jarak antara dua titik pada permukaan spheroid, yang dikembangkan oleh Thaddeus Vincenty (1975a). Rumus ini didasarkan pada asumsi bahwa sosok Bumi adalah spheroid oblate, dan karenanya lebih akurat daripada metode seperti jarak lingkaran besar yang mengasumsikan Bumi bulat.

Perbedaan akurasi ~0.17%dalam jarak 428 meter di Israel. Saya telah membuat tes kecepatan cepat dan kotor:

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

Kode:

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)

Untuk menyimpulkan: rumus Vincenty adalah dua kali lipat waktu perhitungan dibandingkan dengan lingkaran besar, dan perolehan akurasinya pada titik yang diuji adalah ~ 0,17%.

Karena waktu perhitungan dapat diabaikan, rumus Vincenty lebih disukai untuk setiap kebutuhan praktis.

Pembaruan : Mengikuti komentar mendalam dari whuber dan cffk dan cffk , saya setuju bahwa perolehan akurasi harus dibandingkan dengan kesalahan, bukan pengukuran. Oleh karena itu, rumus Vincenty adalah beberapa pesanan yang lebih besar, tidak ~ 0,17%.

Jika Anda menggunakan geopy, maka jarak great_circle dan vincenty sama-sama nyaman untuk didapatkan. Dalam hal ini, Anda hampir selalu harus menggunakan yang memberi Anda hasil yang lebih akurat, yaitu, vincenty. Dua pertimbangan (seperti yang Anda tunjukkan) adalah kecepatan dan ketepatan.

Vincenty dua kali lebih lambat. Tetapi mungkin dalam aplikasi nyata peningkatan waktu berjalan diabaikan. Sekalipun aplikasi Anda meminta perhitungan sejuta jarak, kami hanya membicarakan perbedaan waktu beberapa detik.

Untuk poin yang Anda gunakan, kesalahan dalam vincenty adalah 6 m dan kesalahan dalam jarak lingkaran besar adalah 0,75 m. Saya kemudian akan mengatakan bahwa vincenty adalah 120000 kali lebih akurat (daripada 0,17% lebih akurat). Untuk poin umum, kesalahan dalam jarak lingkaran besar bisa sebanyak 0,5%. Jadi bisakah Anda hidup dengan kesalahan 0,5% dalam jarak? Untuk penggunaan biasa (berapa jarak dari Cape Town ke Kairo?), Mungkin Anda bisa. Namun, banyak aplikasi GIS memiliki persyaratan akurasi yang lebih ketat. (0,5% 5m lebih dari 1 km. Itu benar-benar membuat perbedaan.)

Hampir semua pekerjaan pemetaan serius dilakukan pada ellipsoid referensi dan karenanya masuk akal bahwa jarak juga harus diukur pada ellipsoid. Mungkin Anda bisa pergi dengan jarak lingkaran besar hari ini. Tetapi untuk setiap aplikasi baru, Anda harus memeriksa apakah ini masih dapat diterima. Lebih baik hanya menggunakan jarak ellipsoidal dari awal. Anda akan tidur lebih nyenyak di malam hari.

ADDENDUM (Mei 2017)

Sebagai balasan untuk jawaban yang diberikan oleh @ craig-hicks. Metode vincenty () dalam geopy memang memiliki kesalahan fatal: ia melempar kesalahan untuk hampir poin antipodal. Dokumentasi dalam kode menunjukkan peningkatan jumlah iterasi. Tetapi ini bukan solusi umum karena metode berulang yang digunakan oleh vincenty () tidak stabil untuk poin-poin tersebut (setiap iterasi membawa Anda lebih jauh dari solusi yang benar).

Mengapa saya mencirikan masalah ini sebagai "berpotensi fatal"? Karena setiap penggunaan fungsi jarak dalam pustaka perangkat lunak lain harus dapat menangani pengecualian. Menanganinya dengan mengembalikan NaN atau jarak lingkaran besar mungkin tidak memuaskan, karena fungsi jarak yang dihasilkan tidak akan mematuhi ketimpangan segitiga yang menghalangi penggunaannya, misalnya, dalam pohon titik-pandang.

Situasinya tidak sepenuhnya suram. Paket python saya geographiclib menghitung jarak geodesik secara akurat tanpa kegagalan. The geopy tarik permintaan # 144 berubah fungsi jarak geopy untuk paket penggunaan geographiclib jika tersedia. Sayangnya permintaan penarikan ini telah di limbo sejak Augest 2016.

ADDENDUM (Mei 2018)

geopy 1.13.0 sekarang menggunakan paket geographiclib untuk menghitung jarak. Berikut ini contoh panggilan (berdasarkan contoh pada pertanyaan awal):

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

Permintaan maaf saya untuk mengirim jawaban kedua di sini, tetapi saya mengambil kesempatan untuk menanggapi permintaan dengan @raig-hicks untuk memberikan perbandingan akurasi dan waktu untuk berbagai algoritma untuk menghitung jarak geodesi. Ini memparafrasekan komentar yang saya buat pada permintaan tarikan saya # 144 untuk geopy yang memungkinkan penggunaan salah satu dari dua implementasi algoritma saya untuk geodesics untuk digunakan dalam geopy, satu adalah implementasi python asli, geodesik (geographiclib) , dan kegunaan lainnya sebuah implementasi di C, geodesic (pyproj) .

Berikut adalah beberapa data waktu. Waktu dalam mikron per panggilan

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

Berikut adalah akurasi perhitungan geodesik berdasarkan Set Tes Geodesik saya . Kesalahan diberikan dalam satuan mikron (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.010

Saya telah memasukkan permintaan tarik hannosche # 194 yang memperbaiki bug buruk di fungsi tujuan. Tanpa perbaikan ini, kesalahan dalam perhitungan tujuan untuk vincenty adalah 8,98 meter.

19,2% dari kasus tes gagal dengan vincenty.distance (iterasi = 20). Namun set tes condong ke kasus yang akan menyebabkan kegagalan ini.

Dengan titik acak pada ellipsoid WGS84, algoritma Vincenty dijamin gagal 16,6 dari 1000000 kali (solusi yang benar adalah titik tetap tidak stabil dari metode Vincenty).

Dengan implementasi geopy Vincenty dan iterasi = 20, tingkat kegagalan adalah 82,8 per 10.00000. Dengan iterasi = 200, tingkat kegagalan adalah 21,2 per 10.00000.

Meskipun angka ini kecil, kegagalan bisa sangat umum. Misalnya dalam dataset 1000 titik acak (mungkin bandara dunia, mungkin), menghitung matriks jarak penuh akan gagal rata-rata 16 kali (dengan iterasi = 20).

Tampaknya paket geopy.distance menawarkan fungsi "distance ()" yang default ke vincenty (). Saya akan merekomendasikan menggunakan jarak () pada prinsipnya, karena itu adalah rekomendasi paket, dalam kasus ini pernah menyimpang dari vincenty () di masa depan (tidak mungkin seperti itu). Lanjut membaca:

Catatan dokumentasi ini termasuk dalam kode sumber untuk fungsi vincenty () yang Anda tentukan:

Catatan: Implementasi jarak Vincenty ini gagal menyatu untuk beberapa titik yang valid. Dalam beberapa kasus, hasilnya dapat diperoleh dengan meningkatkan jumlah iterasi ( iterationsargumen kata kunci, yang diberikan di kelas __init__, dengan default 20). Mungkin lebih baik menggunakan: class:, .great_circleyang sedikit kurang akurat, tetapi selalu membuahkan hasil.

Kode sumber dengan komentar / catatan di atas dapat ditemukan di https://github.com/geopy/geopy/blob/master/geopy/distance.py Gulir ke bawah ke definisi untuk vincenty ()

Namun demikian, fungsi jarak default yang digunakan oleh paket itu ketika caliing distance () adalah fungsi vincenty (), yang menyiratkan bahwa kegagalan untuk konvergen bukan merupakan bencana, dan jawaban yang masuk akal dikembalikan - yang paling penting pengecualian tidak dihasilkan.

Pembaruan: Seperti dicatat oleh "cffk", fungsi vincenty () tidak secara eksplisit melempar pengecualian ValueError ketika algoritma tidak konvergen - meskipun tidak didokumentasikan dalam deskripsi fungsi. Oleh karena itu, dokumentasinya buggy.

Apakah menggunakan vincenty atau haversine atau hukum kosinus bulat, ada kebijaksanaan untuk mengetahui setiap masalah potensial dengan kode yang Anda rencanakan untuk digunakan, hal-hal yang harus diperhatikan dan dikurangi, dan bagaimana seseorang berurusan dengan vincenty vs haversine vs masalah sloc akan berbeda ketika seseorang menyadari masalah / edgecases masing-masing, yang mungkin atau mungkin tidak dikenal. Pemrogram berpengalaman tahu ini. Pemula mungkin tidak. Saya berharap dapat menghindarkan sebagian dari mereka frustrasi ketika cuplikan dari forum melakukan sesuatu yang tidak terduga, dalam kasus tertentu. Jika seseorang benar-benar akan menggunakan beberapa versi dari semua ini, vincenty, haversine, sloc, maka SE, SO, Reddit, Quora, dll, mungkin telah memberikan bantuan terbatas dalam beberapa pengkodean awal dari suatu solusi, tetapi itu tidak berarti bahwa solusi mereka atau 'jawaban' yang diterima bebas dari masalah. Jika suatu proyek cukup penting, maka layak untuk dilakukan sejumlah penelitian yang wajar. Baca manual, baca dokumen, dan jika ada tinjauan kode dari kode itu, baca itu. Menyalin dan menempelkan potongan atau intisari yang telah di-upgrade seratus kali atau lebih tidak berarti keamanannya komprehensif dan terjamin.

Jawaban yang menarik yang diposting oleh cffk meningkatkan kesadaran akan edgecases, dalam solusi paket, yang dapat menghasilkan pengecualian atau kesulitan lain . Klaim khusus yang dibuat dalam pos itu berada di luar anggaran waktu saya untuk mengejar saat ini, tetapi saya mengambil darinya bahwa memang ada masalah yang mengintai dalam paket tertentu, termasuk setidaknya satu implementasi vincenty, mengenai yang setidaknya satu orang telah mengusulkan untuk meningkatkan dengan satu atau lain cara, untuk meminimalkan atau menghilangkan risiko menghadapi kesulitan-kesulitan itu. Saya tidak akan menambahkan lebih jauh ke topik itu mengenai vincenty (terlalu bodoh tentang hal itu), tetapi sebaliknya akan beralih ke haversine, setidaknya sebagian pada topik dengan OP.

Rumus haversine yang diterbitkan secara populer, baik dengan python atau bahasa lain, karena itu kemungkinan besar akan menggunakan spesifikasi floating point IEEE 754 pada sebagian besar semua sistem seperti intel dan intel saat ini, dan prosesor ARM, powerPC, dll. juga rentan terhadap kesalahan pengecualian yang jarang tetapi nyata dan berulang yang sangat dekat atau pada jarak busur 180 derajat, titik-titik antipodal, karena pendekatan dan pembulatan titik apung. Beberapa pemula mungkin belum tergigit oleh situasi ini. Karena spesifikasi fp ini mendekati dan bulat, ini tidak berarti bahwa kode apa pun yang memanggil fp64 dapat menyebabkan kesalahan pengecualian, tidak. Tetapi beberapa kode, beberapa rumus mungkin tidak memiliki edgecases yang sangat jelas di mana perkiraan dan pembulatan IEEE 754 fp64 dapat menyebabkan nilai menyimpang sedikit dari domain metode matematika yang diharapkan dengan sempurna mengevaluasi nilai tersebut. Contoh ... sqrt (). Jika nilai negatif menemukan jalannya menjadi sqrt (), seperti sqrt (-0.00000000000000000000122739), akan ada kesalahan pengecualian. Dalam rumus haversine, cara pengembangannya menuju solusi, ada dua metode sqrt () di atan2 (). Itua yang dihitung dan kemudian digunakan dalam sqrt (), dapat, pada titik-titik antipodal di dunia, sedikit tersesat di bawah 0,0 atau di atas 1,0, sangat sedikit karena pendekatan fp64 dan pembulatan, jarang, tetapi berulang. Pengulangan yang konsisten dan dapat diandalkan, dalam konteks ini, menjadikan ini sebagai risiko pengecualian, suatu dasar untuk melindungi, memitigasi, dan bukan kebetulan acak yang terisolasi. Berikut adalah contoh potongan python3 pendek haversine, tanpa perlindungan yang diperlukan:

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

Sangat dekat atau pada titik-titik antipodal, sebuah dihitung di baris pertama dari rumus mungkin nyasar negatif, jarang, tapi repeatably dengan koordinat tersebut lon lat yang sama. Untuk melindungi / memperbaiki mereka kejadian langka, satu hanya dapat menambahkan, setelah sebuah perhitungan, seperti yang terlihat di bawah ini:

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

Tentu saja saya tidak menunjukkan seluruh fungsi di sini, tetapi potongan pendek seperti yang sering diposting. Tapi ini menunjukkan perlindungan untuk sqrt (), dengan menguji a , dan menormalkannya jika perlu, juga menghemat kebutuhan untuk mencoba semuanya kecuali. Note = '' di bagian atas adalah untuk mencegah tahap bytecode dari memprotes note yang sedang digunakan sebelum diberi nilai, jika dikembalikan dengan hasil fungsi.

Dengan perubahan sederhana ini, menambahkan dua sebuah tes, sqrt () fungsi akan senang, dan kode sekarang memiliki tambahan catatan yang dapat dikembalikan ke memanggil kode, untuk peringatan bahwa hasilnya telah sedikit dinormalisasi, dan mengapa. Beberapa mungkin peduli, beberapa mungkin tidak, tetapi ada di sana, mencegah kesalahan pengecualian, yang 'dapat' terjadi. Coba kecuali blok dapat menangkap pengecualian, tetapi tidak memperbaikinya, kecuali secara tertulis ditulis untuk melakukannya. Tampaknya lebih mudah untuk kode garis koreksi (s) segera setelah sebuah garis perhitungan. Masukan yang telah digosok harus tidak perlu dicoba kecuali memblokir di sini sama sekali.

Singkatnya, jika menggunakan haversine, kode eksplisit daripada menggunakan paket atau perpustakaan, tidak peduli bahasa pilihan Anda, itu akan menjadi ide yang baik untuk menguji dan untuk menormalkan sebuah kembali ke kisaran yg diperlukan dari 0,0 <= a <= 1,0 dalam rangka untuk melindungi baris berikutnya dengan perhitungan c . Tetapi sebagian besar cuplikan kode haversine tidak menunjukkannya, dan tidak menyebutkan risikonya.

Pengalaman: selama pengujian menyeluruh di seluruh dunia, dengan kenaikan 0,001 derajat, saya telah mengisi hard drive dengan kombinasi lat yang menyebabkan pengecualian, pengecualian berulang yang dapat diandalkan dan konsisten, selama sebulan juga menguji secara reliabilitas keandalan pendinginan CPU penggemar, dan kesabaran saya. Ya, saya sejak itu menghapus sebagian besar log itu, karena tujuan mereka sebagian besar adalah untuk membuktikan maksudnya (jika pun kata diizinkan). Tapi saya punya beberapa log pendek 'nilai masalah lat lon', disimpan untuk tujuan pengujian.

Akurasi: Akankah suatu dan seluruh hasil haversine kehilangan beberapa akurasi dengan menormalkan bahwa sedikit kembali kecil ke domain? Tidak banyak, mungkin tidak lebih dari perkiraan fp64 dan pembulatan sudah diperkenalkan, yang menyebabkan sedikit keluar dari domain. Jika Anda telah menemukan haversine dapat diterima lebih dari vincenty - lebih sederhana, lebih cepat, lebih mudah untuk menyesuaikan, memecahkan masalah dan memelihara, maka haversine mungkin merupakan solusi yang baik untuk proyek Anda.

Saya telah menggunakan haversine pada skysphere yang diproyeksikan overhead untuk mengukur jarak sudut antara objek di langit, seperti yang dilihat dari posisi di bumi, pemetaan azimuth dan alt ke skysphere lat lon seperti koordinat, tidak ada elipsoid untuk dipertimbangkan sama sekali, karena proyeksi langit skysphere adalah bola sempurna, ketika datang untuk mengukur jarak sudut melihat sudut antara dua objek dari posisi di permukaan bumi. Ini sesuai dengan kebutuhan saya dengan sempurna. Jadi, haversine masih sangat berguna, dan sangat akurat, dalam aplikasi tertentu (baik untuk keperluan saya) ... tetapi jika Anda menggunakannya, apakah di bumi untuk GIS atau navigasi, atau dalam pengamatan dan pengukuran objek langit, lakukan perlindungan dalam hal titik antipodal atau titik antipodal yang sangat dekat, dengan menguji adan mendorongnya kembali ke domain yang diperlukan saat dibutuhkan.

Haversine yang tidak terlindungi ada di internet, dan saya hanya melihat satu pos usenet lama yang menunjukkan perlindungan, saya pikir dari seseorang di JPL, dan itu mungkin pra-1985, pra-IEEE 754 floating point spec. Dua halaman lain menyebutkan kemungkinan masalah di dekat poin antipodal, tetapi tidak menjelaskan masalah itu, atau bagaimana orang dapat mengatasinya. Jadi ada kekhawatiran bagi para pemula (seperti saya) yang mungkin tidak selalu memahami praktik yang baik dengan cukup baik untuk penelitian lebih lanjut, dan menguji edgecases, dari beberapa kode yang telah mereka salin dan tempelkan ke dalam proyek dengan kepercayaan. Posting menarik cffk menyegarkan karena bersifat publik dengan jenis masalah ini, yang tidak sering disebutkan, jarang dikodekan secara publik untuk perlindungan dalam cuplikan, dan jarang dibahas dengan cara ini, dibandingkan dengan jumlah versi yang tidak dilindungi dan tidak didiskusikan yang diposting.

Pada 20190923, halaman wiki untuk formula haversine memang menyebutkan masalah yang mungkin terjadi di titik-titik antipodal, karena masalah titik mengambang di perangkat komputasi ... mendorong ...

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

(karena halaman wiki itu, pada saat ini, tidak memiliki jangkar html untuk bagian yang akan saya tautkan secara langsung, oleh karena itu, setelah halaman dimuat, lakukan pencarian pada halaman browser itu untuk 'Ketika menggunakan formula ini' dan Anda akan lihat masalah haversine dengan poin antipodal yang disebutkan, lebih resmi.)

Dan situs lain ini juga memiliki penyebutan yang sangat singkat:

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

Jika seseorang menemukan di halaman itu untuk 'termasuk perlindungan terhadap kesalahan pembulatan', ada ini ...

Jika atan2 tidak tersedia, c dapat dihitung dari 2 ⋅ asin (min (1, √a)) (termasuk perlindungan terhadap kesalahan pembulatan).

Sekarang ada contoh langka di mana kesalahan pembulatan disebutkan, dan perlindungan ditampilkan untuk versi asin (), namun tidak disebutkan atau ditampilkan untuk versi atan2 (). Tetapi setidaknya risiko kesalahan pembulatan disebutkan.

imho, aplikasi 24/7/365 apa pun yang menggunakan haversine, memerlukan perlindungan ini di dekat titik antipodal sebagai detail penting dan sederhana.

Saya tidak tahu paket haversine mana yang termasuk atau tidak termasuk perlindungan ini, tetapi jika Anda masih baru dalam hal ini, dan Anda akan menggunakan versi 'snippet' yang dipublikasikan secara populer, sekarang Anda tahu itu membutuhkan perlindungan, dan perlindungan itu sangat sederhana untuk diterapkan, yaitu, jika Anda tidak menggunakan vincenty, dan tidak menggunakan haversine paket tanpa akses mudah untuk mengubah kode paket.

TKI, apakah menggunakan vincenty atau haversine atau sloc, orang harus menyadari masalah apa pun dengan kode, hal-hal yang harus diperhatikan dan dimitigasi, dan bagaimana seseorang berurusan dengan vincenty vs haversine vs sloc akan berbeda ketika seseorang mengetahui masing-masing masalah tersembunyi / edgecases, yang mungkin atau mungkin tidak dikenal secara luas.