Skip to content

Instantly share code, notes, and snippets.

@tomowarkar

tomowarkar/distance.md

Last active July 24, 2022 01:27
Show Gist options
  • Save tomowarkar/dd483a6aaaa6e32a8c8ccef2a908ebb5 to your computer and use it in GitHub Desktop.
Save tomowarkar/dd483a6aaaa6e32a8c8ccef2a908ebb5 to your computer and use it in GitHub Desktop.
Download ZIP
緯度軽度でNearestNeighborsを使う
Raw
distance.md

2地点間距離を求めたい

2地点間距離

import math
from sklearn.neighbors import DistanceMetric


def sklearn_haversine_distance(from_, to_):
    (y1, x1), (y2, x2) = map(math.radians, from_), map(math.radians, to_)
    distance = DistanceMetric.get_metric('haversine').pairwise    
    return distance([[y1, x1], [y2, x2]])[1, 0]*6378.137


def haversine_distance(from_, to_):
    (y1, x1), (y2, x2) = map(math.radians, from_), map(math.radians, to_)
    dx, dy = x2 - x1, y2 - y1
    s, c = math.sin, math.cos
    
    dist = math.asin(math.sqrt(
        s(dy/2)**2+c(y1)*c(y2)*s(dx/2)**2
    ))
    r = 6378.137
    return dist*2*r

def lambert_andoyer_distance(from_, to_):
    (y1, x1), (y2, x2) = map(math.radians, from_), map(math.radians, to_)
    
    a = 6378.140
    b = 6356.755
    f = 1 - b / a
    s, c = math.sin, math.cos

    pa = math.atan(b / a * math.tan(y1))
    pb = math.atan(b / a * math.tan(y2))
    xx = math.acos(s(pa) * s(pb) + c(pa) * c(pb) * c(x1 - x2))
    c1 = (s(xx) - xx) * (s(pa) + s(pb)) ** 2 / c(xx / 2) ** 2
    c2 = (s(xx) + xx) * (s(pa) - s(pb)) ** 2 / s(xx / 2) ** 2
    dr = f / 8 * (c1 - c2)
    return a * (xx + dr)

性能評価

import pandas as pd


d = {
    '札幌': (43.064301, 141.346869),
    '東京': (35.689608, 139.692080),
    '福岡': (33.606316, 130.418108),
    'シドニー': (-33.856960, 151.215109),
    'ワシントン': (38.897668, -77.036680),
    'ロンドン': (51.501157, -0.142491),
}

df = pd.DataFrame.from_dict(d, orient = "index", columns=['lat', 'lng'])
lat lng
札幌 43.0643 141.347
東京 35.6896 139.692
福岡 33.6063 130.418
シドニー -33.857 151.215
ワシントン 38.8977 -77.0367
ロンドン 51.5012 -0.142491 
from itertools import combinations, product


answer_df = pd.DataFrame(combinations(df.index, 2), columns=['from', 'to'])
answer_df['distance'] = [831, 1417, 8577, 10140, 8889, 881, 7792, 10928, 9583, 7778, 11494, 9417, 15709, 16990, 5913]


for func in [sklearn_haversine_distance, haversine_distance, lambert_andoyer_distance]:
    name = func.__name__.replace('_distance', '')
    
    def apply_func(x):
        ret = func(df.loc[x['from']], df.loc[x['to']])
        return ret, (ret-x['distance'])/x['distance']*100
    
    answer_df[[name, f'error_{name}']] = answer_df.apply(apply_func, axis=1, result_type='expand')


answer_df.round(2)
from to distance sklearn_haversine error_sklearn_haversine haversine error_haversine lambert_andoyer error_lambert_andoyer
0 札幌 東京 831 833.15 0.26 833.15 0.26 831.03 0
1 札幌 福岡 1417 1418.5 0.11 1418.5 0.11 1417.1 0.01
2 札幌 シドニー 8577 8621.55 0.52 8621.55 0.52 8576.73 -0
3 札幌 ワシントン 10140 10126.8 -0.13 10126.8 -0.13 10139.9 -0
4 札幌 ロンドン 8889 8874 -0.17 8874 -0.17 8888.94 -0
5 東京 福岡 881 879.99 -0.11 879.99 -0.11 880.66 -0.04
6 東京 シドニー 7792 7834.18 0.54 7834.18 0.54 7791.77 -0
7 東京 ワシントン 10928 10916.5 -0.11 10916.5 -0.11 10927.9 -0
8 東京 ロンドン 9583 9570.55 -0.13 9570.55 -0.13 9583.44 0
9 福岡 シドニー 7778 7818.2 0.52 7818.2 0.52 7777.89 -0
10 福岡 ワシントン 11494 11483.5 -0.09 11483.5 -0.09 11494.2 0
11 福岡 ロンドン 9417 9405.24 -0.12 9405.24 -0.12 9417.03 0
12 シドニー ワシントン 15709 15726.9 0.11 15726.9 0.11 15708.5 -0
13 シドニー ロンドン 16990 17013.5 0.14 17013.5 0.14 16989.9 -0
14 ワシントン ロンドン 5913 5904.16 -0.15 5904.16 -0.15 5912.87 -0

garbage collection

# 動く
def great_circular_distance(lat1, lng1, lat2, lng2):
    y1, x1, y2, x2 = map(math.radians, (lat1, lng1, lat2, lng2))
    dx = x2 - x1
    s, c = math.sin, math.cos
    r = 6378.137
    
    dist = math.acos(
        s(y1)*s(y2) + c(y1)*c(y2)*c(dx)
    )
    return dist*r

# 結果があやしい
def hybeny_distance(lat1, lng1, lat2, lng2):
    y1, x1, y2, x2 = map(math.radians, (lat1, lng1, lat2, lng2))
    dx, dy = x2 - x1, y2 - y1
    my = (y1+y2)/2

    rx = 6378.137            # 赤道半径(m) WGS84
    ry = 6356.752314        # 極半径(m)   WGS84
    e2=(rx*rx-ry*ry)/rx/rx  # 離心率 E^2

    w = math.sqrt(1-e2*math.sin(my)*math.sin(my));
    m = rx*(1-e2)/w**3
    n = rx/w
    return math.sqrt((dy*m)**2 + (dx*n*math.cos(my))**2)
Raw
main.md

緯度軽度でNearestNeighborsを使う

sklearn.neighbors

import pandas as pd
import numpy as np
from sklearn.neighbors import NearestNeighbors

d = {
    '北海道': [43.06417,141.34694],
    '東京': [35.686991,139.539242],
    '愛知': [35.002511,137.208724],
    '大阪': [34.598366,135.545261],
    '沖縄': [26.477084,127.922927],
    'パリ': [49.0083899664, 2.53844117956]
}


df = pd.DataFrame.from_dict(d).T.reset_index()
df.columns = ['CITY', 'lat', 'lng']
CITY lat lng
0 北海道 43.0642 141.347
1 東京 35.687 139.539
2 愛知 35.0025 137.209
3 大阪 34.5984 135.545
4 沖縄 26.4771 127.923
5 パリ 49.0084 2.53844 
R = 6378 # Earth radius
N = 3

neigh = NearestNeighbors(n_neighbors=N+1, metric='haversine')
neigh.fit(np.radians(df[['lat', 'lng']]))

def haversine(x):
    latlng = x[['lat', 'lng']].map(np.radians)
    [dist], [idx] = neigh.kneighbors([latlng])
    return list(df.loc[idx[1:], 'CITY']) + list(dist[1:]*R)

df[['%s%d' % (['neighbor_', 'dist_'][i//N], i%N) for i in range(N*2)]] = df.apply(haversine, axis=1, result_type='expand')
CITY lat lng neighbor_0 neighbor_1 neighbor_2 dist_0 dist_1 dist_2
0 北海道 43.0642 141.347 東京 愛知 大阪 835.747 965.773 1067.6
1 東京 35.687 139.539 愛知 大阪 北海道 224.903 383.186 835.747
2 愛知 35.0025 137.209 大阪 東京 北海道 158.566 224.903 965.773
3 大阪 34.5984 135.545 愛知 東京 北海道 158.566 383.186 1067.6
4 沖縄 26.4771 127.923 大阪 愛知 東京 1161.61 1298.66 1506.94
5 パリ 49.0084 2.53844 北海道 大阪 愛知 9027.44 9633.68 9666.13

On my own

def calc(lat1, lng1, lat2, lng2):
    y1, x1, y2, x2 = map(np.radians, (lat1, lng1, lat2, lng2))
    dx = x2 - x1
    
    dist = np.arccos(
        np.sin(y1)*np.sin(y2) + np.cos(y1)*np.cos(y2)*np.cos(dx)
    )
    theta = np.arctan2(
        np.cos(y1)*np.tan(y2)-np.sin(y1)*np.cos(dx),
        np.sin(dx)
    )
    return dist, np.pi/2-theta
    

a, b = calc(*d['北海道'], *d['東京'])
print('距離 [km]:', a*R, '\n角度:',np.degrees(b))

# 距離 [km]: 835.7472348357185 
# 角度: 191.30841594687126

reference

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment