Python中的Haversine公式（两个GPS点之间的方位和距离）答案

【问题标题】：Haversine Formula in Python (Bearing and Distance between two GPS points)Python中的Haversine公式（两个GPS点之间的方位和距离）
【发布时间】：2011-06-22 06:01:28
【问题描述】：

问题

我想知道如何获得 2 个 GPS 点之间的距离和方位。我研究了haversine公式。有人告诉我，我也可以使用相同的数据找到方位。

编辑

一切正常，但轴承还不能完全正常工作。轴承输出负值，但应在 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 "--------------------"

【问题讨论】：

Python hasrsine 实现可以在codecodex.com/wiki/… 找到。然而，对于短距离计算，存在非常简单的方法。现在，您期望的最大距离是多少？你能在一些当地的笛卡尔坐标系中得到你的坐标吗？
python中的一些实现：-code.activestate.com/recipes/…-platoscave.net/blog/2009/oct/5/…
@James Dyson：距离像 15 公里，创造圈不算什么。我的建议：首先找出欧几里得距离的解决方案！这将为您提供一个可行的解决方案，然后如果您的距离会更长，然后调整您的应用程序。谢谢
@James Dyson：如果您的上述评论是针对我的（以及针对我之前的建议），那么答案肯定是（而且非常“微不足道”）。我也许可以给出一些示例代码，但它不会使用三角函数，而是几何（所以我不确定它是否会对你有所帮助。你是否熟悉向量的概念？在你的情况下，位置和方向可以用向量最直接的方式处理）。
atan2(sqrt(a), sqrt(1-a)) 与asin(sqrt(a)) 相同

标签： python gps distance haversine bearing

【解决方案1】：

这是一个 Python 版本：

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

def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance in kilometers 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. Determines return value units.
    return c * r

【讨论】：

可以使用 math.radians() 函数而不是乘以 pi/180 - 效果相同，但更多的是自我记录。
你可以，但如果你说import math，那么你必须指定math.pi、math.sin等。使用from math import *，你可以直接访问所有模块内容。在 python 教程中查看“命名空间”（例如docs.python.org/tutorial/modules.html）
你为什么使用 atan2(sqrt(a), sqrt(1-a)) 而不是 asin(sqrt(a))？在这种情况下 atan2 是否更准确？
如果平均地球半径定义为 6371 公里，那么这相当于 3959 英里，而不是 3956 英里。有关计算这些值的各种方法，请参阅 Global average radii。
返回的是什么？方位还是距离？

【解决方案2】：

您可以通过添加 360° 来解决负方位问题。不幸的是，这可能导致正向轴承的轴承大于 360°。这是模运算符的一个很好的候选，所以总而言之，你应该添加行

Bearing = (Bearing + 360) % 360

在你的方法结束时。

【讨论】：

我认为只是：轴承 = 轴承 % 360
这也是使用增强赋值的好时机：Bearing %= 360

【解决方案3】：

方位计算不正确，需要将输入交换到 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

这将为您提供正确的方位。

【讨论】：

在阅读论文时，我实际上很难理解这些方程是如何推导出来的。你给了我指点：haversine formula我第一次听到这个，谢谢。

【解决方案4】：

你可以试试 hasrsine 包： https://pypi.org/project/haversine/

示例代码：

from haversine import haversine
haversine((45.7597, 4.8422),(48.8567, 2.3508), unit='mi')
243.71209416020253

【讨论】：

如何在 Django 的 ORM 查询中使用它？

【解决方案5】：

默认情况下，atan2 中的 Y 是第一个参数。这是documentation。您需要切换输入以获得正确的方位角。

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

【讨论】：

【解决方案6】：

参考此链接：https://gis.stackexchange.com/questions/84885/whats-the-difference-between-vincenty-and-great-circle-distance-calculations

这实际上提供了两种获取距离的方法。他们是Haversine 和Vincentys。根据我的研究，我知道 Vincentys 是相对准确的。也可以使用import语句来实现。

【讨论】：

【解决方案7】：

这里有两个计算距离和方位的函数，它们基于之前消息中的代码和https://gist.github.com/jeromer/2005586（为清楚起见，为这两个函数添加了纬度、经度格式的地理点的元组类型）。我测试了这两个功能，它们似乎都可以正常工作。

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

【讨论】：

此方法提供的结果与上述所有其他方法不同！

【解决方案8】：

这些答案中的大多数都是“四舍五入”地球的半径。如果您将这些与其他距离计算器（例如 geopy）进行对比，这些功能将被关闭。

这很好用：

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

【讨论】：

这个比上面的例子准确多了！
这并没有解决两极的 R. 6356.752 公里到赤道的 6378.137 公里的变化
这个错误对您的应用程序真的很重要吗？ cs.nyu.edu/visual/home/proj/tiger/gisfaq.html

【解决方案9】：

这是 @Michael Dunn 给出的 Haversine 公式的 numpy 矢量化实现，比大型矢量提高了 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

【讨论】：

【解决方案10】：

还有一个矢量化实现，它允许使用 4 个 numpy 数组而不是坐标的标量值：

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

【讨论】：