Python：计算两个纬度/经度之间的方位

Posted 2023-02-16

【中文标题】Python：计算两个纬度/经度之间的方位

【英文标题】：Python: Calculate bearing between two lat/long

【发布时间】：2019-07-19 07:25:35

【问题描述】：

我正在尝试计算两个纬度/经度之间的方位。

我没有关于函数/公式本身的问题，

提供：

def get_bearing(lat1, long1, lat2, long2):
    dLon = (long2 - long1)

    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)

    brng = np.rad2deg(brng)

    return brng

问题是结果不是预期的。

该函数的预期用途返回（非常长）列表中两个纬度/经度对之间的方位，即

    lat1 = path[int(len(path) * location / 1000)][0]
    lat2 = path[int(len(path) * location / 1000) + 1][0]
    lng1 = path[int(len(path) * location / 1000)][1]
    lng2 = path[int(len(path) * location / 1000) + 1][1]

然后方位角结果会改变绘图的视图方向，其中方位角可以采用 [-180, 180] 范围内的值。理想情况下，结果应该是 lat1、lng1 和 lat2、lng2 之间形成的线在图中完全“垂直”（lat/lon 注释在图中切换），见下文

我希望有人能够从函数返回的方位以及预期的方位中推断出问题所在。以下几个例子：

Current Location: 30.07134 -97.23076
Next in path: 30.0709 -97.22907
Calculated Bearing: 88.39967863143139
Expected Bearing: ~-70.67

Current Location: 29.91581 -96.85068
Next in path: 29.91556 -96.85021
Calculated Bearing: 118.9170342272798
Expected Bearing: ~122.67

Current Location: 29.69419 -96.53487
Next in path: 29.69432 -96.53466
Calculated Bearing 141.0271357781952
Expected Bearing: ~56

Current Location: 29.77357 -96.07924
Next in path: 29.77349 -96.07876
Calculated Bearing 165.24612555483893
Expected Bearing: ~104

很高兴提供更多信息，在此先感谢您的任何/所有帮助。

【参考方案1】：

您是否考虑过使用pyproj 进行计算而不是自己滚动？：

import pyproj
geodesic = pyproj.Geod(ellps='WGS84')
fwd_azimuth,back_azimuth,distance = geodesic.inv(long1, lat1, long2, lat2)

在此示例中，fwd_azimuth 是您所追求的方位，back_azimuth 是反向方位（方向相反）。

我在这里使用 WGS84，因此您需要替换为正确的坐标系，并且需要重写以确保 lat/long 是 geodesic.inv() 的正确坐标类型。但是使用经过充分测试的现有地理空间库可能会为您省去很多麻烦。

【讨论】：

我收到了一个nan 为fwd_azimuth 返回的float(lat)、float(lon)...等 - 我肯定想使用地理空间库，但是我不想处理特殊问题，nan 上有什么想法吗？

您输入的纬度/经度对使用什么坐标系？查看API docs for pyproj 以了解它们如何表示纬度/经度。就像我上面说的，你必须修改你的代码以这种格式传递 lat/long。

tbh，idk - 我正在使用 Plot.ly 的 ScatterMapBox 对象，但我在他们的参考资料中找不到任何内容：plot.ly/python/reference/#scattermapbox。我猜它和 GoogleMaps 一样，因为他们的 API 可以很好地与ScatterMapBox 配合使用。 Google Maps and Microsoft Virtual Earth use a Mercator projection based on the World Geodetic System (WGS) 1984 geographic coordinate system (datum).

我找到了一个我喜欢并且有效的库（geographiclib）！我会赞成你的回答，因为你引导我朝那个方向发展

你的 lats 和 longs 是错误的，例如见pyproj4.github.io/pyproj/stable/api/geod.html： >>> from pyproj import Geod >>> g = Geod(ellps='clrk66') # Use Clarke 1866 椭球体。 >>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.) >>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.) >>> # 计算前后方位角，加上距离 >>> # 在波士顿和波特兰之间。 >>> az12,az21,dist = g.inv(boston_lon,boston_lat,portland_lon,portland_lat) >>> "%7.3f %6.3f %12.3f" % (az12,az21,dist) '-66.531 75.654 4164192.708'

【参考方案2】：

最终改变了功能：

from geographiclib.geodesic import Geodesic
...
def get_bearing(lat1, lat2, long1, long2):
    brng = Geodesic.WGS84.Inverse(lat1, long1, lat2, long2)['azi1']
    return brng

【讨论】：

此文档链接解释了从 Geodesic.WGS84.Inverse geographiclib.sourceforge.io/1.52/python/interface.html 返回的字典

【参考方案3】：

查找轴承的代码很好。但是你只需要在找到 X 和 Y 时添加 math.radians 即可。

import numpy
import math

def get_bearing(lat1, long1, lat2, long2):
    dLon = (long2 - long1)
    x = math.cos(math.radians(lat2)) * math.sin(math.radians(dLon))
    y = math.cos(math.radians(lat1)) * math.sin(math.radians(lat2)) - math.sin(math.radians(lat1)) * math.cos(math.radians(lat2)) * math.cos(math.radians(dLon))
    brng = numpy.arctan2(x,y)
    brng = numpy.degrees(brng)

    return brng

【讨论】：

代码 sn-p 缺少导入 - 从哪里导入 arctan2 和 degrees，应该是 numpy 吗？

是的，它是从 numpy 导入的。对不起，我没有把它放在那里。 @Wolkenarchitekt

【参考方案4】：

你可以试试这个名为 Turfpy 的新包，它是 turf.js 在 python 中的移植。

from turfpy import measurement
from geojson import Point, Feature
start = Feature(geometry=Point((-75.343, 39.984)))
end = Feature(geometry=Point((-75.534, 39.123)))
measurement.bearing(start,end)

https://pypi.org/project/turfpy/