Formula Haversine dengan Python (Bearing dan Jarak antara dua titik GPS)

120

Masalah

Saya ingin tahu bagaimana mendapatkan jarak dan arah antara 2 titik GPS . Saya telah meneliti formula haversine. Seseorang memberi tahu saya bahwa saya juga dapat menemukan arah menggunakan data yang sama.

Edit

Semuanya berfungsi dengan baik tetapi bantalannya belum berfungsi dengan baik. Output bearing negatif tetapi harus antara 0 - 360 derajat. Data yang ditetapkan harus membuat bantalan horizontal 96.02166666666666 dan:

Start point: 53.32055555555556 , -1.7297222222222221   
Bearing:  96.02166666666666  
Distance: 2 km  
Destination point: 53.31861111111111, -1.6997222222222223  
Final bearing: 96.04555555555555

Ini kode baru saya:

from math import *

Aaltitude = 2000
Oppsite  = 20000

lat1 = 53.32055555555556
lat2 = 53.31861111111111
lon1 = -1.7297222222222221
lon2 = -1.6997222222222223

lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])

dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * atan2(sqrt(a), sqrt(1-a))
Base = 6371 * c


Bearing =atan2(cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(lon2-lon1), sin(lon2-lon1)*cos(lat2)) 

Bearing = degrees(Bearing)
print ""
print ""
print "--------------------"
print "Horizontal Distance:"
print Base
print "--------------------"
print "Bearing:"
print Bearing
print "--------------------"


Base2 = Base * 1000
distance = Base * 2 + Oppsite * 2 / 2
Caltitude = Oppsite - Aaltitude

a = Oppsite/Base
b = atan(a)
c = degrees(b)

distance = distance / 1000

print "The degree of vertical angle is:"
print c
print "--------------------"
print "The distance between the Balloon GPS and the Antenna GPS is:"
print distance
print "--------------------"
avitex
sumber
Implementasi Python haversine dapat ditemukan codecodex.com/wiki/… . Namun untuk perhitungan jarak pendek ada cara yang sangat sederhana. Sekarang, berapa jarak maksimum yang Anda harapkan? Apakah Anda bisa mendapatkan koordinat Anda di beberapa sistem koordinat kartesian lokal?
makan
Beberapa implementasi dalam python: - code.activestate.com/recipes/… - platoscave.net/blog/2009/oct/5/…
Fábio Diniz
1
@ James Dyson: dengan jarak seperti 15km, lingkaran pencipta tidak menghitung apa-apa. Saran saya: cari tahu dulu solusinya dengan jarak euclidean! Itu akan memberi Anda solusi yang berfungsi dan kemudian jika jarak Anda akan jauh lebih panjang, maka sesuaikan aplikasi Anda. Terima kasih
makan
1
@ James Dyson: Jika komentar Anda di atas ditujukan untuk saya (dan untuk saran saya sebelumnya), jawabannya pasti (dan cukup 'sepele' juga). Saya mungkin dapat memberikan beberapa kode contoh, tetapi itu tidak akan menggunakan trigonometri, melainkan geometri (jadi saya tidak yakin apakah itu akan membantu Anda sama sekali. Apakah Anda paham sama sekali dengan konsep vektor? Dalam posisi kasus dan arah dapat ditangani dengan cara paling mudah dengan vektor).
makan
1
atan2(sqrt(a), sqrt(1-a))sama denganasin(sqrt(a))
user102008

Jawaban:

241

Ini versi Python:

from math import radians, cos, sin, asin, sqrt

def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance between two points 
    on the earth (specified in decimal degrees)
    """
    # convert decimal degrees to radians 
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])

    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a)) 
    r = 6371 # Radius of earth in kilometers. Use 3956 for miles
    return c * r
Michael Dunn
sumber
10
Bisa menggunakan fungsi math.radians () daripada mengalikan dengan pi / 180 - efek yang sama, tetapi sedikit lebih mendokumentasikan diri.
Hugh Bothwell
4
Anda bisa, tetapi jika Anda berkata import mathmaka Anda harus menentukan math.pi, math.sindll. Dengan from math import *Anda mendapatkan akses langsung ke semua konten modul. Lihat "ruang nama" dalam tutorial python (seperti docs.python.org/tutorial/modules.html )
Michael Dunn
2
Kenapa Anda menggunakan atan2 (akar (a), akar (1-a)), bukan hanya asin (akar (a))? Apakah atan2 lebih akurat dalam kasus ini?
Eyal
4
Jika radius Bumi rata-rata didefinisikan sebagai 6371 km, maka itu setara dengan 3959 mil, bukan 3956 mil. Lihat Jari - jari rata-rata global untuk berbagai cara menghitung nilai-nilai ini.
ekhumoro
3
apa ini kembali? Bantalan atau jarak?
AesculusMaximus
11

Sebagian besar jawaban ini "membulatkan" jari-jari bumi. Jika Anda memeriksanya dengan kalkulator jarak lain (seperti geopy), fungsi ini akan dimatikan.

Ini bekerja dengan baik:

from math import radians, cos, sin, asin, sqrt

def haversine(lat1, lon1, lat2, lon2):

      R = 3959.87433 # this is in miles.  For Earth radius in kilometers use 6372.8 km

      dLat = radians(lat2 - lat1)
      dLon = radians(lon2 - lon1)
      lat1 = radians(lat1)
      lat2 = radians(lat2)

      a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2
      c = 2*asin(sqrt(a))

      return R * c

# Usage
lon1 = -103.548851
lat1 = 32.0004311
lon2 = -103.6041946
lat2 = 33.374939

print(haversine(lat1, lon1, lat2, lon2))
Tanah liat
sumber
2
Yang ini jauh lebih akurat daripada contoh di atas!
Alex van Es
Ini tidak mengatasi variasi R. 6356.752 km di kutub hingga 6378.137 km di ekuator
ldmtwo
3
Apakah kesalahan itu penting untuk aplikasi Anda? cs.nyu.edu/visual/home/proj/tiger/gisfaq.html
Tejas Kale
8

Ada juga implementasi vektor , yang memungkinkan untuk menggunakan 4 larik numpy alih-alih nilai skalar untuk koordinat:

def distance(s_lat, s_lng, e_lat, e_lng):

   # approximate radius of earth in km
   R = 6373.0

   s_lat = s_lat*np.pi/180.0                      
   s_lng = np.deg2rad(s_lng)     
   e_lat = np.deg2rad(e_lat)                       
   e_lng = np.deg2rad(e_lng)  

   d = np.sin((e_lat - s_lat)/2)**2 + np.cos(s_lat)*np.cos(e_lat) * np.sin((e_lng - s_lng)/2)**2

   return 2 * R * np.arcsin(np.sqrt(d))
Sergey Malyutin
sumber
5

Perhitungan bantalan salah, Anda perlu menukar input ke atan2.

    bearing = atan2(sin(long2-long1)*cos(lat2), cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(long2-long1))
    bearing = degrees(bearing)
    bearing = (bearing + 360) % 360

Ini akan memberi Anda arah yang benar.

Jon Anderson
sumber
Saya sebenarnya berjuang untuk memahami bagaimana persamaan ini diturunkan saat saya membaca makalah. Anda telah memberi saya petunjuk: haversine formulapertama kali saya mendengar ini, terima kasih.
arilwan
ini benar, dan mod 360 adalah sentuhan yang bagus
Marc Compere
4

Anda dapat mencoba yang berikut ini:

from haversine import haversine
haversine((45.7597, 4.8422),(48.8567, 2.3508), unit='mi')
243.71209416020253
Vamshi G
sumber
Bagaimana ini bisa digunakan dalam kueri ORM Django?
Gocht
3

Berikut ini adalah implementasi vektor numpy dari Haversine Formula yang diberikan oleh @Michael Dunn, memberikan peningkatan 10-50 kali lipat dibandingkan vektor besar.

from numpy import radians, cos, sin, arcsin, sqrt

def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance between two points 
    on the earth (specified in decimal degrees)
    """

    #Convert decimal degrees to Radians:
    lon1 = np.radians(lon1.values)
    lat1 = np.radians(lat1.values)
    lon2 = np.radians(lon2.values)
    lat2 = np.radians(lat2.values)

    #Implementing Haversine Formula: 
    dlon = np.subtract(lon2, lon1)
    dlat = np.subtract(lat2, lat1)

    a = np.add(np.power(np.sin(np.divide(dlat, 2)), 2),  
                          np.multiply(np.cos(lat1), 
                                      np.multiply(np.cos(lat2), 
                                                  np.power(np.sin(np.divide(dlon, 2)), 2))))
    c = np.multiply(2, np.arcsin(np.sqrt(a)))
    r = 6371

    return c*r
Shubham Singh Yadav
sumber
2

Anda dapat menyelesaikan soal baring negatif dengan menambahkan 360 °. Sayangnya, hal ini mungkin menghasilkan bantalan yang lebih besar dari 360 ° untuk bantalan positif. Ini adalah kandidat yang baik untuk operator modulo, jadi secara keseluruhan Anda harus menambahkan baris

Bearing = (Bearing + 360) % 360

di akhir metode Anda.

OBu
sumber
1
Saya pikir itu hanya: Bearing = Bearing% 360
Holger Bille
1

Y di atan2, secara default, adalah parameter pertama. Berikut dokumentasinya . Anda perlu mengganti input Anda untuk mendapatkan sudut bantalan yang benar.

bearing = atan2(sin(lon2-lon1)*cos(lat2), cos(lat1)*sin(lat2)in(lat1)*cos(lat2)*cos(lon2-lon1))
bearing = degrees(bearing)
bearing = (bearing + 360) % 360
gisdude
sumber
0

Berikut adalah dua fungsi untuk menghitung jarak dan arah, yang didasarkan pada kode di pesan sebelumnya dan https://gist.github.com/jeromer/2005586 (menambahkan tipe tupel untuk titik geografis dalam format lat, lon untuk kedua fungsi untuk kejelasan ). Saya menguji kedua fungsi tersebut dan tampaknya berfungsi dengan benar.

#coding:UTF-8
from math import radians, cos, sin, asin, sqrt, atan2, degrees

def haversine(pointA, pointB):

    if (type(pointA) != tuple) or (type(pointB) != tuple):
        raise TypeError("Only tuples are supported as arguments")

    lat1 = pointA[0]
    lon1 = pointA[1]

    lat2 = pointB[0]
    lon2 = pointB[1]

    # convert decimal degrees to radians 
    lat1, lon1, lat2, lon2 = map(radians, [lat1, lon1, lat2, lon2]) 

    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a)) 
    r = 6371 # Radius of earth in kilometers. Use 3956 for miles
    return c * r


def initial_bearing(pointA, pointB):

    if (type(pointA) != tuple) or (type(pointB) != tuple):
        raise TypeError("Only tuples are supported as arguments")

    lat1 = radians(pointA[0])
    lat2 = radians(pointB[0])

    diffLong = radians(pointB[1] - pointA[1])

    x = sin(diffLong) * cos(lat2)
    y = cos(lat1) * sin(lat2) - (sin(lat1)
            * cos(lat2) * cos(diffLong))

    initial_bearing = atan2(x, y)

    # Now we have the initial bearing but math.atan2 return values
    # from -180° to + 180° which is not what we want for a compass bearing
    # The solution is to normalize the initial bearing as shown below
    initial_bearing = degrees(initial_bearing)
    compass_bearing = (initial_bearing + 360) % 360

    return compass_bearing

pA = (46.2038,6.1530)
pB = (46.449, 30.690)

print haversine(pA, pB)

print initial_bearing(pA, pB)
Oleksiy Muzalyev
sumber
metode ini memberikan hasil selain semua metode lain di atas!
basilisk