非長方形の2Dグリッド上の最も近い点のインデックスを効率的に見つける

Question

不規則な（長方形ではない）lon/latグリッドと、lon/lat座標の点群があります。数値的な理由から少しずれている可能性があります）。今、対応するlon/lat点のインデックスが必要です。

私はこれを行う関数を書いていますが、それは本当に遅いです。

def find_indices(lon,lat,x,y): 
    lonlat = np.dstack([lon,lat]) 
    delta = np.abs(lonlat-[x,y]) 
    ij_1d = np.linalg.norm(delta,axis=2).argmin() 
    i,j = np.unravel_index(ij_1d,lon.shape) 
    return i,j 

ind = [find_indices(lon,lat,p*) for p in points] 

numpy/scipyには、優れた（さらに速い）ソリューションがあると確信しています。私はすでにかなりたくさんのグーグルで探偵してきましたが、答えは今まで私を逃しています。

対応する（最も近い）ポイントのインデックスを効率的に見つける方法を教えてください。

PS：この質問は別のもの（click）から現れました。

編集：ここでは、大局的にこのソリューションともDivakarの答えから1を配置するには

def find_indices(points,lon,lat,tree=None): 
    if tree is None: 
     lon,lat = lon.T,lat.T 
     lonlat = np.column_stack((lon.ravel(),lat.ravel())) 
     tree = sp.spatial.cKDTree(lonlat) 
    dist,idx = tree.query(points,k=1) 
    ind = np.column_stack(np.unravel_index(idx,lon.shape)) 
    return [(i,j) for i,j in ind] 

です：@Cong馬の答えに基づくソリューション

、私は次の解決策を見つけました私は（上記のリンクを参照）find_indicesを使用して（と、それはスピードの面でのボトルネックだ）いている機能のいくつかのタイミング：

spatial_contour_frequency/pil0    : 331.9553 
spatial_contour_frequency/pil1    : 104.5771 
spatial_contour_frequency/pil2    :  2.3629 
spatial_contour_frequency/pil3    :  0.3287 

私の最初のアプローチは、pil1Divakar'sであり、pil2/pil3である。上記の最終的な解決策では、ツリーは、オンザフライでpil2で作成される。ループの反復ごとにfind_indicesが呼び出されます）、一度だけpil3（詳細は他のスレッドを参照してください）。私の最初のアプローチのDivakarの洗練は私に3倍のスピードアップをもたらしますが、cKDTreeはこれを全く新しいレベルまで50倍のスピードアップで実現します！そして、関数から木の作成を動かすことは事をより速くします。

Answer 1

ポイントは十分に局在化している場合は、自分in another postで説明したように、あなたは、のscipy.spatial直接cKDTree実施をしようとします。そのポストは補間に関するものでしたが、あなたはそれを無視してクエリ部分を使うことができます。

TL; DRのバージョン：

がscipy.sptial.cKDTreeのドキュメントを読んでください。初期化子に(n, m)字numpy ndarrayオブジェクトを渡すことによって、ツリーを作成し、ツリーがnm次元座標から作成されます。

tree = scipy.spatial.cKDTree(array_of_coordinates) 

その後、（おそらく近似と並列して、ドキュメントを参照してください）k番目の最近傍を検索、または指定された距離の許容範囲内のすべてのネイバーを見つけるためにtree.query_ball_point()を使用するtree.query()を使用しています。

ポイントがうまくローカライズされていない場合、および球面曲率/非密着空間キックでは、ローカルとみなされるのに十分なそれぞれが小さい、複数の部分にマニホールドを破壊してみてください。

Answer 2

ここscipy.spatial.distance.cdist使用して、一般的なベクトル化されたアプローチだ - サンプル実行

import scipy 

# Stack lon and lat arrays as columns to form a Nx2 array, where is N is grid**2 
lonlat = np.column_stack((lon.ravel(),lat.ravel())) 

# Get the distances and get the argmin across the entire N length 
idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0) 

# Get the indices corresponding to grid's shape as the final output 
ind = np.column_stack((np.unravel_index(idx,lon.shape))).tolist() 

を -

In [161]: lon 
Out[161]: 
array([[-11. , -7.82 , -4.52 , -1.18 , 2.19 ], 
     [-12. , -8.65 , -5.21 , -1.71 , 1.81 ], 
     [-13. , -9.53 , -5.94 , -2.29 , 1.41 ], 
     [-14.1 , -0.04 , -6.74 , -2.91 , 0.976]]) 

In [162]: lat 
Out[162]: 
array([[-11.2 , -7.82 , -4.51 , -1.18 , 2.19 ], 
     [-12. , -8.63 , -5.27 , -1.71 , 1.81 ], 
     [-13.2 , -9.52 , -5.96 , -2.29 , 1.41 ], 
     [-14.3 , -0.06 , -6.75 , -2.91 , 0.973]]) 

In [163]: lonlat = np.column_stack((lon.ravel(),lat.ravel())) 

In [164]: idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0) 

In [165]: np.column_stack((np.unravel_index(idx,lon.shape))).tolist() 
Out[165]: [[0, 4], [0, 4], [0, 4], [0, 4], [0, 4], [0, 4], [3, 3]] 

ランタイムテスト -

関数を定義します。

def find_indices(lon,lat,x,y): 
    lonlat = np.dstack([lon,lat]) 
    delta = np.abs(lonlat-[x,y]) 
    ij_1d = np.linalg.norm(delta,axis=2).argmin() 
    i,j = np.unravel_index(ij_1d,lon.shape) 
    return i,j 

def loopy_app(lon,lat,pts): 
    return [find_indices(lon,lat,pts[i,0],pts[i,1]) for i in range(pts.shape[0])] 

def vectorized_app(lon,lat,points): 
    lonlat = np.column_stack((lon.ravel(),lat.ravel())) 
    idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0) 
    return np.column_stack((np.unravel_index(idx,lon.shape))).tolist() 

のタイミング：もっとパフォーマンスを搾り出すために

In [179]: lon = np.random.rand(100,100) 

In [180]: lat = np.random.rand(100,100) 

In [181]: points = np.random.rand(50,2) 

In [182]: %timeit loopy_app(lon,lat,points) 
10 loops, best of 3: 47 ms per loop 

In [183]: %timeit vectorized_app(lon,lat,points) 
10 loops, best of 3: 16.6 ms per loop 

、np.concatenateはnp.column_stackの代わりに使用することができます。