2014-07-20 70 views
1

我正在嘗試實現Lat和Lon數據的最近鄰居搜索。這裏是DATA.TXT距離大於零的Python KD樹最近的neigbour

61.3000183105 -21.2500038147 0 
62.299987793 -23.750005722 1 
66.3000488281 -28.7500038147 2 
40.8000183105 -18.250005722 3 
71.8000183105 -35.7500038147 3 
39.3000183105 -19.7500019073 4 
39.8000183105 -20.7500038147 5 
41.3000183105 -20.7500038147 6 

的問題是,當我想要做的近鄰每個緯度和經度上的數據集,它正在尋找它的自我。例如最近鄰(-21.2500038147,61.3000183105)將是(-21.2500038147,61.3000183105),結果距離將爲0.0。我試圖避免這個,但沒有運氣。我試圖做的,如果沒有(array_equal),但仍...

下面是我的Python代碼

import numpy as np 
from numpy import * 
import decimal 
from scipy import spatial 
from scipy.spatial import KDTree 
from math import radians,cos,sin,sqrt,exp 


Lat =[] 
Lon =[] 
Day =[] 

nja = [] 


Data = np.loadtxt('Data.txt',delimiter=" ") 
for i in range(0,len(Data)): 
    Lon.append(Data[i][:][0]) 
    Lat.append(Data[i][:][1]) 
    Day.append(Data[i][:][2]) 

tree =spatial.KDTree(zip(Lon,Lat)) 

print "Lon :",len(Lon) 
print "Tree :",len(tree.data) 

for i in range(0,len(tree.data)): 
    pts = np.array([tree.data[i][0],tree.data[i][1]]) 
    nja.append(pts) 

for i in range(0, len(nja)): 
    if not (np.array_equal(nja,tree.data)): 
    nearest = tree.query(pts,k=1,distance_upper_bound =9) 
    print nearest 

回答

-1

How'bout低科技的解決方案?如果你有大量的點(比如10000個或更多),這是沒有比較合理的,但對於一個小數目這個蠻力解決方案可能是有用的:

import numpy as np 

dist = (Lat[:,None]-Lat[None,:])**2 + (Lon[:,None]-Lon[None,:])**2 

現在你有一個N×N的陣列(N爲點的數量)與所有點對之間的距離(或距離的平方,更精確)。找到每個點的最短距離就是找到每一行的最小值。要排除的點本身,您可以設置對角NaN和使用nanargmax

np.fill_diagonal(dist, np.nan) 
closest = np.nanargmin(dist, axis=1) 

這個辦法很簡單,保證找到最近的點,但有兩個顯著缺點:

  1. 這是爲O(n^2),和在10000點花費約爲一秒
  2. OT消耗大量的存儲器(800 MB用於上述情況中)

後一個問題當然可以通過分段來避免,但第一個問題不包括大點集。


此,可以進行還通過使用scipy.spatial.distance.pdist

dist=scipy.spatial.distance.pdist(np.column_stack((Lon, Lat))) 

這是一個快一點(由半至少),但輸出矩陣是濃縮形式,請參見文檔scipy.spatial.distance.squareform

如果您需要計算真實距離,那麼這是一個很好的選擇,因爲pdist可以處理球體上的距離。


然後,再次,你可以僅通過擴展您的查詢到兩個最接近的點使用KDtree方法:

nearest = tree.query(pts, k=2, distance_upper_bound=9) 

然後nearest[1][0]具有點本身(「我,我自己和我」) ,nearest[1][1]真實最近的鄰居(或inf,如果沒有足夠的東西)。

最好的解決方案取決於你擁有的點數。另外,如果地圖點在地球上彼此不接近,則可能需要使用笛卡爾2D距離以外的其他東西。如果你只是試圖假裝他們是2D笛卡爾點,你弄錯了:


只是一個有關使用經度和緯度尋找距離音符。在北緯60度,一度緯度爲1111公里,而一度經度爲555公里。所以,至少你必須用cos(緯度)來劃分經度。即使有了這個竅門,當經度從東向西變化時,你也會陷入困境。

大概是出於這樣的煩惱的最簡單的方法是計算座標點成直角的三維點:

x = cos(lat) * cos(lon) 
y = cos(lat) * sin(lon) 
z = sin(lat) 

如果再計算這些點之間的最短距離,你會得到正確的結果。 (請注意,距離與地球表面真實的最短距離不一樣。)

1

對於數據集中的每個點P[i],您問的是:「我的數據中哪一點最接近P[i]組?」你會得到答案「這是P[i]」。

如果你問一個不同的問題「這是最接近P[i]?兩點」,即tree.query(pts,k=2) (與您的代碼是s/k=1/k=2/差),你會得到P[i],也是一個P[j],第二最近點,那是你想要的結果。

旁註:

  • 我建議您在構建樹,原因在你的緯度範圍有什麼由經度1度的距離意味着一個大的波動之前,項目數據。