У меня есть набор данных, состоящий из долготы и широты, и мне нужно преобразовать эти значения в значения расстояния и азимута. По сути, моя цель - как рассчитать расстояние и азимут между двумя точками широты и долготы.
Я знаю, что есть много подобных вопросов, но почему-то ни один из них не дал мне того, что мне точно нужно. Есть много ресурсов о том, как рассчитать расстояние, но сложно вычислить азимут. Я реализовал то, что нашел в интернете, и также скопировал некоторые решения людей, которые ответили на подобные вопросы.
Вот мой окончательный код:
from math import *
from geopy.distance import geodesic, distance
from geographiclib.geodesic import Geodesic
from haversine import haversine
from pygc import great_distance
import pyproj
lat1, lon1, lat2, lon2 = 48.863846, 8.863657, 48.846757, 8.927386 # Stuttgart: 48.1351, 11.5820
print("distance using the haversine lib => ", haversine(point1=(lat1, lon1), point2=(lat2, lon2)))
geod = pyproj.Geod(ellps='WGS84')
azimuth1, azimuth2, dist = geod.inv(lon1, lat1, lon2, lat2)
print('pyproj distance: ', dist/1000) # pyproj distance: 5.048040324129042
print('pyproj azimuth: ', azimuth1, azimuth2) # pyproj azimuth: 112.09094775492177 -67.86106109849165
d = great_distance(start_latitude=lat1, start_longitude=lon1, end_latitude=lat2, end_longitude=lon2)
print("distance using pygc lib => ", d) # distance using pygc lib => 'distance': 50.48 'azimuth': 112.09094775472433, 'reverse_azimuth': 292.1389389013081}
dis = geodesic((lat1, lon1), (lat2, lon2)).kilometers
print("distance using geopy = ", dis) # distance using geopy = 5.048040324129042
dis1 = distance((lat1, lon1), (lat2, lon2))
print("distance using distance geopy = ", dis1) # distance using distance geopy = 5.048040324129042 km
res = Geodesic.WGS84.Inverse(lat1, lon1, lat2, lon2)
print("distance and azi using geographiclib => ", res['s12']/1000, res['azi1']) # distance and azi using geographiclib => 5.048040324129042 112.09094775492177
print("bearing using geographiclib => ", res['azi1']) # bearing => 112.09094775492177
def get_bearing(lat1, lon1, lat2, lon2, method='fab'):
rad_lat1 = lat1 * pi / 180
rad_lon1 = lon1 * pi / 180
rad_lat2 = lat1 * pi / 180
rad_lon2 = lat1 * pi / 180
if method.lower() == 'fab':
rad_theta = rad_lon2 - rad_lon1
y = sin(rad_theta) * cos(rad_lat2)
x = cos(rad_lat1) * sin(rad_lat2) - sin(rad_lat1) * cos(rad_lat2) * cos(rad_theta)
brng = atan2(y, x) * 180 / pi
return (brng + 360) % 360
elif method.lower() == 'python':
import numpy as np
dLon = rad_lon2 - rad_lon1
y = sin(dLon) * cos(rad_lat2)
x = cos(rad_lat1) * sin(rad_lat2) - sin(rad_lat1) * cos(rad_lat2) * cos(dLon)
brng = np.rad2deg(atan2(y, x))
if brng < 0:
brng += 360
return brng
elif method.lower() == 'so':
bearing = atan2(sin(rad_lon2 - rad_lon1) * cos(rad_lat2),
cos(rad_lat1) * sin(rad_lat2) - sin(rad_lat1) * cos(rad_lat2) * cos(rad_lon2 - rad_lon1))
bearing = degrees(bearing)
bearing = (bearing + 360) % 360
return bearing
elif method.lower() == 'initial':
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)
initial_bearing = degrees(initial_bearing)
compass_bearing = (initial_bearing + 360) % 360
return compass_bearing
return initial_bearing(pointA=(lat1, lon1), pointB=(lat2, lon2))
else:
dLon = rad_lon2 - rad_lon1
y = sin(dLon) * cos(rad_lat2)
x = cos(rad_lat1)*sin(rad_lat2) - sin(rad_lat1)*cos(rad_lat2)*cos(dLon)
brng = atan2(y, x) * 180 / pi
return brng
def get_distance_in_km(lat1, lon1, lat2, lon2, method='fab'):
earth_diameter = 6371 * 2
if method.lower() == 'fab':
fi1 = lat1 * pi / 180
fi2 = lat2 * pi / 180
deltafi = (lat2 - lat1) * pi / 180
deltasigma = (lon2 - lon1) * pi / 180
a = sin(deltafi/2) * sin(deltafi/2) + cos(fi1) * cos(fi2) * sin(deltasigma/2) * sin(deltasigma/2)
return earth_diameter * atan2(sqrt(a), sqrt(1-a))
elif method.lower() == 'haversian':
# 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
return c * r
else:
def haversin(x):
return sin(x / 2) ** 2
return 2 * asin(
sqrt(haversin(lat2 - lat1) + cos(lat1) * cos(lat2) * haversin(lon2 - lon1))
) * (earth_diameter/2)
res = get_distance_in_km(lat1, lon1, lat2, lon2, method='fab')
res1 = get_distance_in_km(lat1, lon1, lat2, lon2, method='haversian')
print("distance using fab method => ", res) # distance using fab method => 5.034895107102602
print("distance using haversian => ", res1) # distance using haversian => 5.034895107102733
br = get_bearing(lat1, lat2, lon1, lon2)
br1 = get_bearing(lat1, lat2, lon1, lon2, method='other')
br2 = get_bearing(lat1, lat2, lon1, lon2, method='python')
br3 = get_bearing(lat1, lat2, lon1, lon2, method='so')
br4 = get_bearing(lat1, lat2, lon1, lon2, method='initial')
print("bearing using the fab method : ", br) # bearing using the fab method : 89.99356472317174
print("bearing using other: ", br1) # bearing using other: 89.99356472317177
print("bearing using python: ", br2) # bearing using python: 89.99356472317177
print("bearing using method found on stackoverflow: ", br3)# bearing using method found on stackoverflow: 89.99356472317174
print("bearing using initial method: ", br4) # bearing using initial method: 233.48863432180315
Как вы можете видеть вкод, я использую разные методы для расчета расстояния и пеленга. Я нашел эти методы в Интернете на разных форумах и в научных статьях. Проблема в том, что они не дают одинаковый выход.
Разница между выходами расстояния очень мала и может быть проигнорирована, но выход азимута сбивает с толку. 3 метода выводят одно и то же значение для азимута, а другие методы выводят другие значения.
Я предоставил выходные данные в комментариях вместе с кодом, но вы можете запустить код и посмотреть, что я имею в виду. Проект, над которым я работаю, требует высокой точности, и поэтому я пробую разные методы, чтобы увидеть, какой из них точнее других, а также библиотеки, которые предоставляют более точные значения, чем другие, как здесь, я пытаюсь geopy, geographiclib, pygc и pyproj.
Вы можете видеть, что эти библиотеки предоставляют разные значения для обратного азимута или обратного азимута, а исходный метод также предоставляет совершенно другое значение. Я использовал этот калькулятор расстояний веб-сайт, чтобы проверить свои результаты. Согласно веб-сайту, значения библиотек более точны, чем реализация.
Это сбивает с толку, так как я реализовал эти функции в соответствии с научными статьями, а некоторые другие детали реализации были найдены на stackoverflow и других форумах. Я очень запутался в этих результатах, и, как я уже сказал, мой Проект требует высокой точности выбора значений, поэтому я должен быть осторожен.
Я надеюсь, что кто-то может подсказать мне об этом и сказать, какой метод или библиотекая должен использовать для этого проекта или какой метод является более точным для расчета значений расстояния и азимута.
Я использую Python 3.7. Вот ссылки на библиотеки, которые я использовал: geopy geographiclib pygc pyproj