Python:计算两个lat/long之间的方位
Python: Calculate bearing between two lat/long
我正在尝试计算两个 lat/long 之间的方位。
我没有关于 function/formula 本身的问题,
提供:
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
问题是结果不是预期的。
函数的预期用途returns(很长的)列表中两个 lat/long 对之间的方位,即
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 之间形成的线在图中完全 "vertical"(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
很高兴提供更多信息,在此先感谢您的 any/all 帮助。
您是否考虑过使用 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()
的正确坐标类型。但是使用经过良好测试的现有地理空间库可能会为您省去很多麻烦。
结束更改功能:
from geographiclib.geodesic import Geodesic
...
def get_bearing(lat1, lat2, long1, long2):
brng = Geodesic.WGS84.Inverse(lat1, long1, lat2, long2)['azi1']
return brng
您可以试试这个名为 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)
找方位的代码没问题。但是你只需要在查找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
我正在尝试计算两个 lat/long 之间的方位。
我没有关于 function/formula 本身的问题,
提供:
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
问题是结果不是预期的。
函数的预期用途returns(很长的)列表中两个 lat/long 对之间的方位,即
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 之间形成的线在图中完全 "vertical"(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
很高兴提供更多信息,在此先感谢您的 any/all 帮助。
您是否考虑过使用 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()
的正确坐标类型。但是使用经过良好测试的现有地理空间库可能会为您省去很多麻烦。
结束更改功能:
from geographiclib.geodesic import Geodesic
...
def get_bearing(lat1, lat2, long1, long2):
brng = Geodesic.WGS84.Inverse(lat1, long1, lat2, long2)['azi1']
return brng
您可以试试这个名为 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)
找方位的代码没问题。但是你只需要在查找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