logo标签列表

Python中的Haversine公式（计算两个GPS点之间的方位角和距离）

pythongpsdistancehaversinebearing
157

问题

我想知道如何计算两个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 Haversine实现可以在http://www.codecodex.com/wiki/Calculate_Distance_Between_Two_Points_on_a_Globe#Python找到。然而，对于短距离计算，存在非常简单的方法。那么，您预期的最大距离是多少？您能否在某些本地笛卡尔坐标系中获取您的坐标？ - eat
一些Python实现：
  • http://code.activestate.com/recipes/576779-calculating-distance-between-two-geographic-points/
  • http://www.platoscave.net/blog/2009/oct/5/calculate-distance-latitude-longitude-python/
 - Fábio Diniz
1
@James Dyson：像15公里这样的距离，创建圆圈不起作用。我的建议是：首先找出欧几里得距离的解决方案！这将为您提供一个可行的解决方案，然后如果您的距离要长得多得多，则调整应用程序。谢谢。 - eat
1
@James Dyson：如果你上面的评论是针对我（以及我的早期建议），那么答案肯定是（而且相当“平凡”）。我可以给出一些示例代码，但它不会利用三角学，而是几何学（所以我不确定它是否会对你有所帮助。你是否熟悉向量的概念？在你的情况下，位置和方向可以用向量最直接的方式处理）。 - eat
1
atan2(sqrt(a), sqrt(1-a))与asin(sqrt(a))相同。 - user102008
显示剩余2条评论
12个回答
322

这是 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
5
你可以这样写：import math，但是你必须要指定使用 math.pi、math.sin 等函数。如果使用 from math import *，则可以直接访问模块中的所有内容。请参阅 Python 教程中的 "命名空间" 部分(例如 http://docs.python.org/tutorial/modules.html)。 - Michael Dunn
2
你为什么使用atan2(sqrt(a), sqrt(1-a))而不是只用asin(sqrt(a))？在这种情况下，atan2更准确吗？ - Eyal
2
应该使用浮点数除法来处理dlat|dlon为整数的极端情况： a = sin(dlat/2.0)**2 + cos(lat1) * cos(lat2) * sin(dlon/2.0)**2 - Dmitriy
4
如果平均地球半径被定义为6371公里，则相当于3959英里，而不是3956英里。请参阅“全球平均半径”了解计算这些值的各种方法。 - ekhumoro
3
这返回的是什么？方位角还是距离？ - AesculusMaximus
显示剩余7条评论
20

大多数答案都是“舍入”地球半径。如果将这些答案与其他距离计算器（如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))
3
这个比上面的例子更准确！ - Alex van Es
5
这并未涉及到R的变化，即从极地的6356.752公里到赤道的6378.137公里。 - ldmtwo
3
那个错误对你的应用程序真的很重要吗？ - Tejas Kale
1
@AlexvanEs 在精度不准确的情况下增加精度没有任何好处。事实上，相反它会给人一种虚假的高精度感。https://dev59.com/KWsy5IYBdhLWcg3w3Rvc#8270846 - Louis Cottereau
关于“这些函数将被关闭”的问题：您能量化一下吗？ 10％？ 1％？ 0.000001％？ - Peter Mortensen
16

还有一种向量化实现方式，它允许使用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))
11

你可以尝试使用haversine包：

示例代码：

from haversine import haversine

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

输出：

243.71209416020253
这个怎么在 Django 的 ORM 查询中使用？ - Gocht
Haversine库是否有计算方位角的函数？ - Anmol Deep
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

这将给您正确的方向。
我实际上正在努力理解这些方程式是如何推导出来的，因为我正在阅读一篇论文。你给了我一个指针：haversine formula，这是我第一次听说，谢谢。 - arilwan
3

这里是一个NumPy向量化实现的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 the 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
2
这段代码不起作用；你从NumPy导入了所有的函数，然后在整个代码中都使用np。 - Dj Sushi
2
你可以通过加上360°来解决负向轴承问题。不幸的是，这可能会导致正向轴承大于360°。这是模数运算符的一个好例子，因此总体而言，你应该添加以下行:
Bearing = (Bearing + 360) % 360

在你的方法结束时。

2
我认为只需要这样写：Bearing = Bearing % 360 - Holger Bille
这是使用增强赋值的好时机：Bearing %= 360。 - Michael Dunn
1
默认情况下，atan2中的Y是第一个参数。这里是文档。您需要切换输入以获得正确的方位角。
bearing = atan2(sin(lon2-lon1)*cos(lat2), cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(lon2-lon1))
bearing = degrees(bearing)
bearing = (bearing + 360) % 360
这看起来像是一个虚假的答案。它是吗？"in"应该做什么？ - Peter Mortensen
@PeterMortensen 是的，没错...看起来过去7年中有一个错别字...不过我记得上次看的时候它没错...谢谢你指出来。我已经编辑了它。 - gisdude
1
请参考Vincenty与大圆距离计算的区别。
实际上，这提供了两种获取距离的方式。它们是haversine和Vincentys。经过我的研究，我得知Vincentys相对精确。还需要使用import语句进行实现。
0

考虑到您的目标是测量两点之间的距离（表示为地理坐标），将提供以下三个选项：

  1. Haversine公式

  2. 使用GeoPy测地线距离

  3. 使用GeoPy大圆距离

选项1

哈弗辛公式可以完成这项工作。然而，需要注意的是这种方法是将地球近似为一个球体，因此存在误差 (请参考这个答案) - 因为地球并不是一个球体。

要使用哈弗辛公式，首先需要定义地球的半径。这本身可能会引起一些争议。以下是三个来源：

我将使用值为6371公里作为地球半径的参考。

# Radius of the Earth
r = 6371.0

我们将利用math模块。

在半径之后，我们转向坐标，并开始将坐标转换为弧度，以便使用数学三角函数。为此，我们导入{{link3：math.radians(x)}}并按以下方式使用：

# Import radians from the '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 the '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

一个方法是使用GeoPy的距离，更具体地说，是使用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 it in miles, change `km` to `miles`

[Out]: 19959.6792674

选项3

一个选择是使用GeoPy的距离，更具体地说，是使用great-circle。

我们可以得到以公里或英里为单位的结果（来源）。

# 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 it in km, change `miles` to `km`

[Out]: 536.997990696

注意事项：

  • 由于大圆距离通常使用Haversine公式进行计算(如Willem Hendriks所指出的)，选项1和3是相似的，但使用不同的半径。

    • GeoPy的大圆距离使用地球的球形模型，使用国际大地测量和地球物理学联合会定义的平均地球半径6371.0087714150598公里，约为6371.009公里(对于WGS-84)，结果误差最高可达0.5%[来源]。
1
大圆通常使用Haversine算法实现。因此，它们可以具有不同的半径但是相同。 - Willem Hendriks
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接