Geopy：计算 GPS 航向/方位答案

【问题标题】：Geopy: calculating GPS heading / bearingGeopy：计算 GPS 航向/方位
【发布时间】：2013-07-12 21:29:59
【问题描述】：

第一次在这里发布海报。

我正在对为桥梁检查 ROV 八旋翼收集的 GPS 数据进行一些数据分析。我们使用 3D 扫描激光雷达、立体视觉、INS 和其他一些简洁技术在 ROS 上运行八旋翼。我目前在与Doug Weibel's 设置类似的设置中使用ublox LEA-6T 来收集原始GPS 数据，例如载波相位、多普勒频移和卫星星历。然后我使用一个开源项目RTKLIB对本地NOAA CORS站点进行一些DGPS后处理，以获得cm精度，以便在重建桥梁的3D点云时更好地进行姿态估计。

无论如何，我使用大部分scipy 来统计验证我的测试结果。
不过，特别是对于这部分，我只是在使用：

我一直在使用 geopy 的方便距离函数来研究我的位置协方差与我测量的地面实况的偏移量。几乎没有按摩参数，我可以找到矩阵中每个标准差元素所描绘的每个方向的距离；北、东、上和之间的三个方向。

但是，这些距离是绝对的，并不描述方向。
说：正、负分别与向北或向南相关。

我可以简单地使用纬度和经度来检测方向的极性，
但我希望能够找到所描述距离的精确点对点方位，
因为我相信全局航向的值可能对我目前的应用之外的其他应用有用。

我发现其他人提出了类似的问题
但它似乎假设了一个很大的圆近似
我希望至少使用 WGS-84 椭球模型，或者任何可以在 geopy 中使用的相同模型：
Jump to Calculating distances

任何建议表示赞赏，
-ruffsl

感兴趣的来源：

python-3.3：http:// www.python.org/download/releases/3.3.0/
numpy:http:// www.numpy.org/
地理：https:// code.google.com/p/geopy/
scipy: http:// www.scipy.org/
ublox LEA-6T：http:// www.u-blox.com/en/gps-modules/u-blox-6-timing-module/lea-6t.html
Doug Weibel 的：http:// diydrones.com/profiles/blogs/proof-of-concept-test-extremely-accurate-3d-velocity-measurement
RTKLIB：http:// www.rtklib.com/
NOAA CORS：http:// geodesy.noaa.gov/CORS/
ROS: http:// www.ros.org/wiki/

【问题讨论】：

你说的很多，但问的很少。你到底需要什么？计算两个纬度/经度 WG84 坐标之间方位角的函数？
@AlexWien，是的，简而言之，这就是我的目标。当人们自己坐飞机时，我发现它在从过去的帖子中寻找答案时更有帮助，它可以为其他人提供更相关和可搜索的关键字。
但这对我来说很有帮助，尤其是 RTKlib

标签： python python-3.x gps geo geopy

【解决方案1】：

使用 Python 的 geographiclib 包。这可以计算椭圆体上的距离和方位等等。（可以内插路径、测量区域等）例如，在

pip install geographiclib

你可以的

>>> from geographiclib.geodesic import Geodesic
>>> Geodesic.WGS84.Inverse(-41.32, 174.81, 40.96, -5.50)
{'lat1': -41.32, 'a12': 179.6197069334283, 's12': 19959679.26735382, 'lat2': 40.96, 'azi2': 18.825195123248392, 'azi1': 161.06766998615882, 'lon1': 174.81, 'lon2': -5.5}

这会计算从新西兰惠灵顿 (41.32S 174.81E) 到西班牙萨拉曼卡 (40.96N 5.50W) 的测地线。距离由 s12（19959679 米）给出，初始方位角（方位角）由 azi1（从北顺时针方向 161.067...度）给出。

【讨论】：

如果您提供有关Geodesic.WGS84.Inverse 的结果集的更多详细信息并实际计算问题所要求的 GPS 航向/方位，将会很有帮助。
此文档链接解释了返回的字典。 geographiclib.sourceforge.io/1.52/python/interface.html

【解决方案2】：

两个纬度/经度坐标之间的方位角：(lat1, lon1), (lat2, lon2)

在下面的代码中，lat1,lon1,lat2,lon2 假定为弧度。
将之前的度数转换为弧度。

dLon = lon2 - lon1;
y = Math.sin(dLon) * Math.cos(lat2);
x = Math.cos(lat1)*Math.sin(lat2) -
        Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
brng = Math.atan2(y, x).toDeg();

方位现在在 -180/180 范围内。

归一化为罗盘方位 (0-360)

if brng < 0: brng+= 360

【讨论】：

我确实喜欢优雅的数学，但我相信这是使用球坐标，因此仍然是大圆近似。对于导航和一般方向应用，包括我自己检查极性的目的，我想这与 WG84 相比就足够了。但对于测量和现场检查，这可能过于粗略，或者不适用于局部事件估计。不过很棒。
是的，这些公式假定地球是一个椭球体。更准确的是，您必须使用基于椭球的公式/实现。
我很困惑...为什么 geopy 没有根据计算距离的 vincenty 函数来计算这个函数？

【解决方案3】：

@AlexWien 在 Python 中的回答

import math, numpy as np

def get_bearing(lat1,lon1,lat2,lon2):
    dLon = lon2 - lon1;
    y = math.sin(dLon) * math.cos(lat2);
    x = math.cos(lat1)*math.sin(lat2) - math.sin(lat1)*math.cos(lat2)*math.cos(dLon);
    brng = np.rad2deg(math.atan2(y, x));
    if brng < 0: brng+= 360
    return brng

【讨论】：