12

緯度と経度のデータがあり、位置を含む 2 つの配列間の距離行列を計算する必要があります。これを使用して、緯度と経度が指定された 2 つの場所の間の距離を取得しました。

これが私のコードの例です：

import numpy as np
import math

def get_distances(locs_1, locs_2):
    n_rows_1 = locs_1.shape[0]
    n_rows_2 = locs_2.shape[0]
    dists = np.empty((n_rows_1, n_rows_2))
    # The loops here are inefficient
    for i in xrange(n_rows_1):
        for j in xrange(n_rows_2):
            dists[i, j] = get_distance_from_lat_long(locs_1[i], locs_2[j])
    return dists


def get_distance_from_lat_long(loc_1, loc_2):

    earth_radius = 3958.75

    lat_dif = math.radians(loc_1[0] - loc_2[0])
    long_dif = math.radians(loc_1[1] - loc_2[1])
    sin_d_lat = math.sin(lat_dif / 2)
    sin_d_long = math.sin(long_dif / 2)
    step_1 = (sin_d_lat ** 2) + (sin_d_long ** 2) * math.cos(math.radians(loc_1[0])) * math.cos(math.radians(loc_2[0])) 
    step_2 = 2 * math.atan2(math.sqrt(step_1), math.sqrt(1-step_1))
    dist = step_2 * earth_radius

    return dist

私の期待される出力はこれです：

>>> locations_1 = np.array([[34, -81], [32, -87], [35, -83]])
>>> locations_2 = np.array([[33, -84], [39, -81], [40, -88], [30, -80]])
>>> get_distances(locations_1, locations_2)
array([[ 186.13522573,  345.46610882,  566.23466349,  282.51056676],
       [ 187.96657622,  589.43369894,  555.55312473,  436.88855214],
       [ 149.5853537 ,  297.56950329,  440.81203371,  387.12153747]])

パフォーマンスは私にとって重要でありCython、ループを高速化するために使用できることの 1 つですが、そこに行く必要がなければいいと思います。

このようなことができるモジュールはありますか？または他の解決策はありますか？

4

4 に答える 4

12

あなたが使用しているHaversine方程式には、最適ではないものがたくさんあります。その一部をトリムして、計算する必要があるサイン、コサイン、および平方根の数を最小限に抑えることができます。以下は私が思いついた中で最高のもので、私のシステムでは、1000 要素と 2000 要素の 2 つのランダムな配列で、Ophion のコード (ベクトル化に関してはほとんど同じことを行います) よりも約 5 倍速く実行されます。

def spherical_dist(pos1, pos2, r=3958.75):
    pos1 = pos1 * np.pi / 180
    pos2 = pos2 * np.pi / 180
    cos_lat1 = np.cos(pos1[..., 0])
    cos_lat2 = np.cos(pos2[..., 0])
    cos_lat_d = np.cos(pos1[..., 0] - pos2[..., 0])
    cos_lon_d = np.cos(pos1[..., 1] - pos2[..., 1])
    return r * np.arccos(cos_lat_d - cos_lat1 * cos_lat2 * (1 - cos_lon_d))

2 つの配列を「そのまま」フィードすると、問題が発生しますが、これはバグではなく、機能です。基本的に、この関数は最後の次元で球の距離を計算し、残りの次元でブロードキャストします。したがって、次のように目的のものを取得できます。

>>> spherical_dist(locations_1[:, None], locations_2)
array([[ 186.13522573,  345.46610882,  566.23466349,  282.51056676],
       [ 187.96657622,  589.43369894,  555.55312473,  436.88855214],
       [ 149.5853537 ,  297.56950329,  440.81203371,  387.12153747]])

しかし、2 つの点リスト間の距離を計算するためにも使用できます。

>>> spherical_dist(locations_1, locations_2[:-1])
array([ 186.13522573,  589.43369894,  440.81203371])

または 2 つの単一点の間:

>>> spherical_dist(locations_1[0], locations_2[0])
186.1352257300577

これは gufuncs の仕組みに着想を得たもので、一度慣れると、さまざまな設定で単一の機能を再利用できる素晴らしい「スイス アーミー ナイフ」コーディング スタイルであることがわかりました。

于 2013-10-16T21:32:20.290 に答える
5

meshgrid を使用して double for ループを置き換えると、より効率的です。

import numpy as np

earth_radius = 3958.75

def get_distances(locs_1, locs_2):
   lats1, lats2 = np.meshgrid(locs_1[:,0], locs_2[:,0])
   lons1, lons2 = np.meshgrid(locs_1[:,1], locs_2[:,1])

   lat_dif = np.radians(lats1 - lats2)
   long_dif = np.radians(lons1 - lons2)

   sin_d_lat = np.sin(lat_dif / 2.)
   sin_d_long = np.sin(long_dif / 2.)

   step_1 = (sin_d_lat ** 2) + (sin_d_long ** 2) * np.cos(np.radians(lats1[0])) * np.cos(np.radians(lats2[0])) 
   step_2 = 2 * np.arctan2(np.sqrt(step_1), np.sqrt(1-step_1))

   dist = step_2 * earth_radius

   return dist
于 2013-10-16T20:56:01.173 に答える
5

これは単にコードをベクトル化するだけです:

def new_get_distances(loc1, loc2):
    earth_radius = 3958.75

    locs_1 = np.deg2rad(loc1)
    locs_2 = np.deg2rad(loc2)

    lat_dif = (locs_1[:,0][:,None]/2 - locs_2[:,0]/2)
    lon_dif = (locs_1[:,1][:,None]/2 - locs_2[:,1]/2)

    np.sin(lat_dif, out=lat_dif)
    np.sin(lon_dif, out=lon_dif)

    np.power(lat_dif, 2, out=lat_dif)
    np.power(lon_dif, 2, out=lon_dif)

    lon_dif *= ( np.cos(locs_1[:,0])[:,None] * np.cos(locs_2[:,0]) )
    lon_dif += lat_dif

    np.arctan2(np.power(lon_dif,.5), np.power(1-lon_dif,.5), out = lon_dif)
    lon_dif *= ( 2 * earth_radius )

    return lon_dif

locations_1 = np.array([[34, -81], [32, -87], [35, -83]])
locations_2 = np.array([[33, -84], [39, -81], [40, -88], [30, -80]])
old = get_distances(locations_1, locations_2)

new = new_get_distances(locations_1,locations_2)

np.allclose(old,new)
True

タイミングを見ると：

%timeit new_get_distances(locations_1,locations_2)
10000 loops, best of 3: 80.6 µs per loop

%timeit get_distances(locations_1,locations_2)
10000 loops, best of 3: 74.9 µs per loop

小さな例では実際には遅くなります。ただし、より大きな例を見てみましょう。

locations_1 = np.random.rand(1000,2)

locations_2 = np.random.rand(1000,2)

%timeit get_distances(locations_1,locations_2)
1 loops, best of 3: 5.84 s per loop

%timeit new_get_distances(locations_1,locations_2)
10 loops, best of 3: 149 ms per loop

これで 40 倍のスピードアップが実現しました。おそらく、いくつかの場所でもう少し速度を上げることができます。

編集: 冗長な場所を削除し、元の場所の配列を変更していないことを明確にするために、いくつかの更新を行いました。

于 2013-10-16T20:56:26.903 に答える
4

Haversine 式は、使用するのに十分な精度を提供しますか? かなりずれることがあります。proj.4、特に python バインディングpyprojを使用すると、精度と速度の両方を得ることができると思います。pyproj は、座標の numpy 配列に対して直接動作できることに注意してください。

于 2013-10-16T22:07:10.593 に答える