#python #gis
#python #гис
Вопрос:
Я пытаюсь придумать функцию, где…
Ввод: геодезическое расстояние в милях или км
Вывод: евклидово расстояние между любыми двумя точками gps, которые находятся на расстоянии ввода друг от друга
Я чувствую, что у меня есть некоторые компоненты
import numpy as np
from numpy import linalg as LA
from geopy.distance import geodesic
loc1 = np.array([40.099993, -83.166000])
loc2 = np.array([40.148652, -82.903962])
Это евклидово расстояние между этими двумя точками
LA.norm(loc1-loc2)
#0.2665175636332336
Это геодезическое расстояние в милях между этими двумя точками
geodesic(loc1,loc2).miles
#14.27909749425243
У моего мозга сейчас не хватает сока, у кого-нибудь есть идеи о том, как я могу создать такую функцию, как:
geodesic_to_euclidean(14.27909749425243)
#0.2665175636332336
Комментарии:
1. Вам нужно расстояние дуги, соединяющей две геолокации, предполагая, что поверхность земли представляет собой сферу?
2. Я смотрю только на диапазоны в пределах 50 миль, поэтому я думаю, что было бы справедливо предположить даже плоскую поверхность?
3.
LA.normв любом случае не даст вам правильного значения, потому что вы не работаете с декартовыми координатами. вам нужно вычислить угол между двумя точками относительно центра земли (в радианах), затем умножить на радиус земли.4. Итак, я использую евклидово расстояние в качестве параметра в алгоритме кластеризации DBSCAN. Единственными измерениями, вводимыми в алгоритм, являются широта и долгота. Поэтому мне не нужно евклидово расстояние для представления чего-либо реального, просто 2-нормальное расстояние между любыми двумя векторами. Это проблематично ?… Я новичок в ГИС…
5. @Jamalan вы теряете точность, чем ближе вы подходите к полюсам. По этой же причине картографическая проекция Меркатора раздувает Антарктиду намного больше, чем она есть на самом деле.
Ответ №1:
Если вас устраивает расстояние по большому кругу, как указано в комментариях, то это должно сработать. Это расстояние haversine:
def haversine(origin, destination, units='mi'):
# Radian deltas
origin_lat = radians(float(origin[0]))
origin_lon = radians(float(origin[1]))
destination_lat = radians(float(destination[0]))
destination_lon = radians(float(destination[1]))
lat_delta = destination_lat - origin_lat
lon_delta = destination_lon - origin_lon
# Radius of earth in meters
r = 6378127
# Haversine formula
a = sin(lat_delta / 2) ** 2 cos(origin_lat) *
cos(destination_lat) * sin(lon_delta / 2) ** 2
c = 2 * asin(sqrt(a))
meters_traveled = c * r
scaling_factors = {
"m:": 1,
"km": 1 / 1000,
"ft": 3.2808, # meters to feet
"mi:": 0.000621371 # meters to miles
}
return meters_traveled * scaling_factors[units]
Если у вас уже есть геодезическое (большое окружное) расстояние в метрах и вам нужна длина хорды, вы можете сделать следующее
def chord(geodesic_distance):
"""
Chord length
C = 2 * r * sin(theta/2)
Arc length; which is geodesic distance in this case
AL = R * theta
therefore
C = 2 * R * sin(AL/(2*R))
"""
r = 6378127 # Radius of earth in meters
return 2 * r * sin(geodesic_distance / (2 * r))
Комментарии:
1. для справки: Формула haversine определяет расстояние по большому кругу между двумя точками на сфере с учетом их долготы и широты.
2. @Jamalan определенно не равномерно распределенная сетка. Я знаю, что это может показаться не тем, о чем вы думаете в своей голове, но то, что вы говорите, — это спросить разницу в длине между хордой и дугой с одинаковым углом. Для малых углов (небольших расстояний по отношению к радиусу земли) эта разница будет минимальной.
3. @Jamalan
cэтот угол уже в радианах. просто умножьте на 180 / пи4. @Jamalan пожалуйста, посмотрите приведенное выше редактирование длины аккорда. Хотя Аарон прав — для близлежащих расстояний на поверхности Земли это будет минимально отличаться
5. Вы, ребята, потрясающие. Это хорошая часть stack overflow. Не циничные, просто полезные люди. Огромное вам обоим спасибо. Я преобразоваю вашу мудрость и код в изящную вспомогательную функцию.