Формула Haversine в Python (азимут и расстояние между двумя точками GPS)

Question

проблема

Я хотел бы знать, как определить расстояние между двумя точками GPS. Я исследовал формулу have rsine. Кто-то сказал мне, что я могу также найти направление, используя те же данные.

редактировать

Все работает нормально, но подшипник пока не совсем работает. Выход подшипника отрицательный, но должен быть в пределах 0 - 360 градусов. Заданные данные должны составлять горизонтальный пеленг 96.02166666666666и является:

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

Вот мой новый код:

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

Answer 1

Вот версия 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

Answer 2

Большинство из этих ответов "округляют" радиус Земли. Если вы сравните их с другими калькуляторами расстояний (такими как геопсия), эти функции будут отключены.

Это хорошо работает:

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


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

print(haversine(lat1, lon1, lat2, lon2))

Answer 3

Существует также векторизованная реализация, которая позволяет использовать 4 числовых массива вместо скалярных значений для координат:

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

Answer 4

Вы можете попробовать следующее:

from haversine import haversine
haversine((45.7597, 4.8422),(48.8567, 2.3508),miles = True)
243.71209416020253

Answer 5

Расчет подшипника неверен, вам нужно поменять местами входы в 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

Это даст вам правильное отношение.

Answer 6

Вот небольшая векторизованная реализация формулы Haversine, которую дал @Michael Dunn, дает улучшение в 10-50 раз по сравнению с большими векторами.

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

Answer 7

Учитывая, что ваша цель — измерить расстояние между двумя точками (представленными географическими координатами), ниже оставим три варианта:

1. Формула гаверсина

2. Использование геодезического расстояния

3. Использование расстояния большого круга GeoPyGeoPy

---

Опция 1

Формула Хаверсина выполнит эту работу, однако важно отметить, что при ее выполнении Земля приближается к сфере, и в этом есть ошибка ( см. Этот ответ ) - поскольку Земля не является сферой.

Чтобы использовать формулу Хаверсина, прежде всего, нужно определить радиус Земли. Это само по себе может привести к некоторым спорам. Учитывая следующие три источника

• Центр космических полетов имени Годдарда НАСА : 6371 км

• Википедия : 6371 км (3958,8 миль)

• Гугл - 6371 км

я буду использовать значение 6371км как ссылка на радиус Земли.

      # Radius of the Earth
r = 6371.0

Мы будем использовать mathмодуль.

После радиуса переходят к координатам и начинают с преобразования координат в радианы, чтобы использовать математические тригонометрические функции . Для этого импортируется math.radians(x)и используйте их следующим образом

      #Import radians from math module
from math import radians

# Latitude and Longitude for the First Point (let's consider 40.000º and 21.000º)
lat1 = radians(40.000)
lon1 = radians(21.000)

# Latitude and Longitude for the Second Point (let's consider 30.000º and 25.000º)
lat2 = radians(30.000)
lon2 = radians(25.000)

Теперь можно применить Формулу Хаверсина. Сначала долгота точки 1 вычитается из долготы точки 2.

      dlon = lon2 - lon1
dlat = lat2 - lat1

Затем и здесь есть пара тригонометрических функций, которые можно использовать, а именно: math.sin(), math.cos(), а также math.atan2(). Мы также будем использовать math.sqrt()

      # Import sin, cos, atan2, and sqrt from math module
from math import sin, cos, atan2, sqrt

a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
c = 2 * atan2(sqrt(a), sqrt(1 - a))
d = r * c

Затем можно получить расстояние, напечатав d.

Поскольку это может помочь, давайте соберем все в функцию (вдохновленный ответом @Michael Dunn )

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

def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great-circle distance (in km) between two points 
    using their longitude and latitude (in degrees).
    """
    # Radius of the Earth
    r = 6371.0

    # Convert degrees to radians 
    # First point
    lat1 = radians(lat1)
    lon1 = radians(lon1)

    # Second Point
    lat2 = radians(lat2)
    lon2 = radians(lon2)

    # Haversine formula 
    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)) 
    return r * c

---

Вариант 2

Один из них будет использовать , а точнее, geodesic.

Мы можем получить результаты как в км, так и в милях ( )

      # Import Geopy's distance
from geopy import distance

wellington = (-41.32, 174.81)
salamanca = (40.96, -5.50)
print(distance.distance(wellington, salamanca).km) # If one wants in miles, change `km` to `miles`

[Out]: 19959.6792674

---

Вариант 3

Один из них будет использовать , а точнее, great-circle.

Мы можем получить результаты как в км, так и в милях ( расстояние GeoPyИсточникрасстояние GeoPyИсточник )

      # Import Geopy's distance
from geopy import distance

newport_ri = (41.49008, -71.312796)
cleveland_oh = (41.499498, -81.695391)

print(distance.great_circle(newport_ri, cleveland_oh).miles) # If one wants in km, change `miles` to `km`

[Out]: 536.997990696

Answer 8

Вы можете решить проблему отрицательных подшипников, добавив 360°. К сожалению, это может привести к тому, что подшипники будут больше 360 ° для положительных подшипников. Это хороший кандидат на оператор по модулю, так что в целом вы должны добавить строку

Bearing = (Bearing + 360) % 360

в конце вашего метода.

Answer 9

Перейдите по этой ссылке: https://gis.stackexchange.com/questions/84885/whats-the-difference-between-vincenty-and-great-circle-distance-calculations

это на самом деле дает два способа получения расстояния. Это Хаверсин и Винсент. Из моего исследования я узнал, что Винсент относительно точен. Также используйте оператор import для реализации.

Answer 10

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

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

Answer 11

Вы можете использовать приведенную ниже реализацию в Python

      import math

def haversine_distance(lat1, lon1, lat2, lon2, unit='K'):
r = 6371  # radius of the earth in kilometers
if unit == 'M':
    r = 3960  # radius of the earth in miles
dLat = math.radians(lat2 - lat1)
dLon = math.radians(lon2 - lon1)
a = math.sin(dLat / 2) * math.sin(dLat / 2) + \
    math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) * \
    math.sin(dLon / 2) * math.sin(dLon / 2)
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
distance = r * c
return distance

Подробнее об этом можно прочитать на сайте Haversine Formula.

Answer 12

Вот две функции для вычисления расстояния и пеленга, которые основаны на коде в предыдущих сообщениях и https://gist.github.com/jeromer/2005586 (для ясности добавлен тип кортежа для географических точек в лат, формате lon для обеих функций). Я проверил обе функции, и они, кажется, работают правильно.

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