Почему при расчете расстояния между точками на земле мои расчеты Хаверсин и Геодезический расходятся? - PullRequest
Новая
       34

Почему при расчете расстояния между точками на земле мои расчеты Хаверсин и Геодезический расходятся?

3 голосов
/ 15 октября 2019

Я получаю дико расходящиеся расстояния, используя два приближения для вычисления расстояния между точками на поверхности Земли. Я использую приближение Haversine ( векторизация ) и более точное (предположительно) geopy.distance.geodesic .

enter image description hereenter image description here

Как вы видите, я отклонился на пять процентов, так как расстояния между точками становятся большими. Это расхождение из-за ошибки округления в Haversine? Действительно ли я доверяю геодезической? Вот код: 

import numpy as np
lat = np.linspace(35,45,100)
lon = np.linspace(-120,-110,100)

data = pd.DataFrame({'Latitude':lat,'Longitude':lon})




def Haversine(v):
    """
    distance between two lat,lon coordinates 
    using the Haversine formula. Assumes one
    radius. r = 3,950 to 3,963 mi 
    """
    from timeit import default_timer as timer
    start = timer()
    R = 3958 # radius at 40 deg 750 m elev
    v = np.radians(v)

    dlat = v[:, 0, np.newaxis] - v[:, 0]
    dlon = v[:, 1, np.newaxis] - v[:, 1]
    c = np.cos(v[:,0,None])

    a = np.sin(dlat / 2.0) ** 2 + c * c.T * np.sin(dlon / 2.0) ** 2

    c = 2 * np.arcsin(np.sqrt(a))
    result = R * c
    print(round((timer() - start),3))
    return result



def slowdistancematrix(data):

    from geopy.distance import geodesic
    distance = np.zeros((data.shape[0],data.shape[0]))
    for i in range(data.shape[0]):

        lat_lon_i = data.Latitude.iloc[i],data.Longitude.iloc[i]

        for j in range(i):

            lat_lon_j = data.Latitude.iloc[j],data.Longitude.iloc[j]

            distance[i,j] = geodesic(lat_lon_i, lat_lon_j).miles
            distance[j,i] = distance[i,j] # make use of symmetry

    return distance

distanceG = slowdistancematrix(data)
distanceH = Haversine(data.values)



plt.scatter(distanceH.ravel(),distanceG.ravel()/distanceH.ravel(),s=.5)
plt.ylabel('Geodesic/Haversine')
plt.xlabel('Haversine distance (miles)')
plt.title('all points in distance matrix')

Я бы лучше использовал векторизованную версию, потому что она быстрая. Тем не менее, 5% слишком велики для меня, чтобы чувствовать себя комфортно с ним. Предполагается, что Haversine только на 0,5%.

ОБНОВЛЕНИЕ:

Обнаружена ошибка. при реализации векторизованной версии я рассчитывал не все расстояния между точками, а только между некоторыми. Я обновил код, чтобы отразить это. Вот какая разница между Haversine и Geodesic для моего домена (25-55 * по -125-1-1):

enter image description here

чертовски хорошо!

Ответы [ 2 ]

1 голос
/ 16 октября 2019

В формуле Хаверсина была ошибка матричной алгебры. Я обновил код в вопросе. Теперь я получаю гораздо лучшее согласие между Haversine и геодезической:

enter image description here

По моему фактическому набору данных:

enter image description here enter image description here

1 голос
/ 15 октября 2019

Формула Хаверсайна вычисляет расстояния между точками на сфере (расстояние большого круга), как и geopy.distance.great_circle.

С другой стороны, geopy.distance.geodesic вычисляет расстояния междууказывает на эллипсоидальную модель Земли, которую вы можете рассматривать как «сплюснутую» сферу.

Разница не в том, чтобы округлять, а в том, что они используют разные формулы, а геодезическая формула более точно моделируетистинная форма земли.

...