• 当前位置：首页
  • IT科技
  • Python 中的 Haversine 公式（两个 GPS 点之间的方位和距离）

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

2024-12-12 08:41:00
admin
原创
197
摘要：问题描述：问题我想知道如何获取两个 GPS 点之间的距离和方位。我研究过半正弦距离。有人告诉我，我也可以使用相同的数据找到方位。一切正常，但方位尚未完全正常工作。方位输出为负值，但应在 0 - 360 度之间。设定的数据应该是水平方位96.02166666666666，并且是：Start point: 53....

问题描述：

问题

我想知道如何获取两个 GPS 点之间的距离和方位。

我研究过半正弦距离。有人告诉我，我也可以使用相同的数据找到方位。

一切正常，但方位尚未完全正常工作。方位输出为负值，但应在 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 "--------------------"

解决方案 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

解决方案 2：

这些答案大部分都是“四舍五入”地球半径。如果你用其他距离计算器（如 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：

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

解决方案 4：

您可以尝试haversine包：

示例代码：

from haversine import haversine

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

输出：

243.71209416020253

解决方案 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

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

解决方案 6：

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

解决方案 7：

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

  1. 半正矢公式

  2. 使用GeoPy测地距离

  3. 使用GeoPy大圆距离

选项 1

半正矢公式可以完成这项工作。但是，需要注意的是，这样做会将地球近似为球体，而这会产生误差（参见此答案）——因为地球不是球体。

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

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

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

  • 谷歌- 6371 公里

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

# Radius of the Earth
r = 6371.0

我们将利用math模块。

了解了半径之后，我们转到坐标，然后开始将坐标转换为弧度，以便使用数学的三角函数。为此，它导入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 kilometers`约为`6371.009 km`（对于`WGS-84`），导致误差高达约`0.5%`[来源]。

解决方案 8：

您可以通过添加 360° 来解决负方位角问题。不幸的是，这可能会导致正方位角大于 360°。这非常适合使用模数运算符，因此总的来说，您应该添加以下行

Bearing = (Bearing + 360) % 360

在您的方法结束时。

解决方案 9：

请参阅Vincenty 与大圆距离计算之间的差异。

这实际上给出了两种获取距离的方法。它们是半正弦和 Vincentys。从我的研究中，我发现 Vincentys 相对准确。还可以使用import语句来实现。

解决方案 10：

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

解决方案 11：

这里有两个函数用于计算距离和方位，它们基于之前消息中的代码以及Python 中两点之间的罗盘方位（为了清晰起见，我为这两个函数添加了纬度、经度格式的地理点元组类型）。我测试了这两个函数，它们似乎工作正常。

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

解决方案 12：

你可以在 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) + \n        math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) * \n    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中阅读更多相关信息

相关推荐
置顶 政府信创国产化的10大政策解读
政府信创国产化的10大政策解读
  政府信创国产化的10大政策解读一、信创国产化的背景与意义信创国产化，即信息技术应用创新国产化，是当前中国信息技术领域的一个重要发展方向。其核心在于通过自主研发和创新，实现信息技术应用的自主可控，减少对外部技术的依赖，并规避潜在的技术制裁和风险。随着全球信息技术竞争的加剧，以及某些国家对中国在科技领域的打压，信创国产化显...
工程项目管理   1565  
置顶 2025年20款好用的项目管理软件推荐，项目管理提效的20个工具和技巧
2025年20款好用的项目管理软件推荐，项目管理提效的20个工具和技巧
  为什么项目管理通常仍然耗时且低效？您是否还在反复更新电子表格、淹没在便利贴中并参加每周更新会议？这确实是耗费时间和精力。借助软件工具的帮助，您可以一目了然地全面了解您的项目。如今，国内外有足够多优秀的项目管理软件可以帮助您掌控每个项目。什么是项目管理软件？项目管理软件是广泛行业用于项目规划、资源分配和调度的软件。它使项...
项目管理软件   1354  
信创国产芯片如何推动自主可控生态建设
信创国产芯片如何推动自主可控生态建设
  信创国产芯片作为信息技术创新的核心领域，对于推动国家自主可控生态建设具有至关重要的意义。在全球科技竞争日益激烈的背景下，实现信息技术的自主可控，摆脱对国外技术的依赖，已成为保障国家信息安全和产业可持续发展的关键。国产芯片作为信创产业的基石，其发展水平直接影响着整个信创生态的构建与完善。通过不断提升国产芯片的技术实力、产...
国产信创系统   21  
信创生态建设中的5个关键挑战与解决方案
信创生态建设中的5个关键挑战与解决方案
  信创生态建设旨在实现信息技术领域的自主创新和安全可控，涵盖了从硬件到软件的全产业链。随着数字化转型的加速，信创生态建设的重要性日益凸显，它不仅关乎国家的信息安全，更是推动产业升级和经济高质量发展的关键力量。然而，在推进信创生态建设的过程中，面临着诸多复杂且严峻的挑战，需要深入剖析并寻找切实可行的解决方案。技术创新难题技...
信创操作系统   27  
信创国产芯片的研发难点与解决方案
信创国产芯片的研发难点与解决方案
  信创产业作为国家信息技术创新发展的重要领域，对于保障国家信息安全、推动产业升级具有关键意义。而国产芯片作为信创产业的核心基石，其研发进展备受关注。在信创国产芯片的研发征程中，面临着诸多复杂且艰巨的难点，这些难点犹如一道道关卡，阻碍着国产芯片的快速发展。然而，科研人员和相关企业并未退缩，积极探索并提出了一系列切实可行的解...
国产化替代产品目录   28  
热门文章
项目管理软件
政府信创国产化的10大政策解读
2025soft
2025年20款好用的项目管理软件推荐，项目管理提效的20个工具和技巧
img020
使用 Python 'Requests' 模块的代理
img044
10个项目管理常见问题解答：包括难点和解决方法
项目管理软件
项目延期的原因有哪些？如何有效预防和管理？
项目里程碑计划模板
项目里程碑计划模板怎么写?设定项目里程碑的5个方法
项目管理软件
项目管理指南：制定项目实施计划的6大步骤有哪些？
项目管理软件有哪些？
热门标签
更多>
IPD流程 项目管理 IT项目管理 敏捷开发 项目管理软件 Scrum 敏捷项目管理 IPD流程管理 项目风险管理 项目管理工具 项目经理 极限编程  产品管理  项目经理职责 项目路线图 项目进度管理 团队管理 测试管理工具 项目预算管理 PMP认证 瀑布式管理 项目范围管理 IT科技 项目管理平台 工程项目管理 项目资源管理 pmp证书考试
云禅道AD
禅道项目管理软件

云端的项目管理软件

尊享禅道项目软件收费版功能

无需维护，随时随地协同办公

内置subversion和git源码管理

每天备份，随时转为私有部署

免费试用