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

这是一个 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

---

这些答案中的大多数都是“四舍五入”地球的半径。如果您将这些与其他距离计算器（例如 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))

---

还有一个矢量化实现，它允许使用 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))

---

您可以尝试使用 hasrsine 包：https://pypi.org/project/haversine/

示例代码：

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

---

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

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

---

这是@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

---

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

Bearing = (Bearing + 360) % 360

在你的方法结束时。

---

默认情况下，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

---

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

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

---

这里有两个计算距离和方位的函数，它们基于之前消息中的代码和 https://gist.github.com/jeromer/2005586（为清楚起见，为两个函数添加了 lat、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)

---

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

使用 GeoPy 测地线距离的 Haversine 公式 使用 GeoPy 大圆距离

选项1

Haversine 公式将完成这项工作，但重要的是要注意，这样做会将地球近似为球体，并且存在错误（see this answer） - 因为地球不是球体。

为了使用Haversine公式，首先需要定义地球的半径。这本身可能会引起一些争议。考虑以下三个来源

美国宇航局戈达德太空飞行中心：6371公里

维基百科：6371 公里（3958.8 英里）

谷歌 - 6371 公里

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

# Radius of the Earth
r = 6371.0

我们将利用 math 模块。

在半径之后，一个移动到坐标，一个开始将坐标转换为弧度，以便使用 math's trigonometric functions。为此导入 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)

现在准备应用Haversine Formula。第一个将点 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's answer 启发）

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's distance，更具体地说，是 geodesic。

我们可以在公里或英里 (Source) 上获得结果

# 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

我们将使用 GeoPy's distance，更具体地说，是 great-circle。

我们可以在公里或英里 (Source) 上获得结果

# 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

---