作者piacere (Beol)
看板Python
標題[問題] 如何使用np.array的經緯度計算兩兩的距離
時間Tue Jan 10 11:33:34 2023
各位先進們好
我是一名基礎不太好的python新手QQ
現有一批大量的經緯度座標要做KNN
但我想添加一個距離限制(例如小於1000m才會被cluster)
參考網路上的做法，我目前是這樣寫的:
假設我有十組經緯度資料，陣列為:
data_point=[[120.228986 22.92753 ]
[120.222007 22.9854525]
[120.21645 22.99625 ]
[120.221625 22.99833 ]
[120.1566975 22.987169 ]
[120.281875 23.106358 ]
[120.314682 23.319719 ]
[120.219985 22.998485 ]
[120.215055 22.99942 ]
[120.20783 22.999415 ]]
def np_getDistance(data, data_point, i, j):
ra = 6378140
rb = 6356755
flatten = 0.003353
midLatA=data_point[i,1]
midLonA=data_point[i,0]
midLatB=data_point[j,1]
midLonB=data_point[j,0]
radLatA = np.radians(midLatA)
radLonA = np.radians(midLonA)
radLatB = np.radians(midLatB)
radLonB = np.radians(midLonB)
pA = np.arctan(rb / ra * np.tan(radLatA))
pB = np.arctan(rb / ra * np.tan(radLatB))
x = np.arccos( np.multiply(np.sin(pA),np.sin(pB)) +
np.multiply(np.multiply(np.cos(pA),np.cos(pB)),np.cos(radLonA - radLonB)))
c1 = np.multiply((np.sin(x) - x) , np.power((np.sin(pA) + np.sin(pB)),2))
/ np.power(np.cos(x / 2),2)
c2 = np.multiply((np.sin(x) + x) , np.power((np.sin(pA) - np.sin(pB)),2))
/ np.power(np.sin(x / 2),2)
dr = flatten / 8 * (c1 - c2)
distance = 0.001 * ra * (x + dr)
return distance
for i in range(dataLen)):
knn = KNN(i, K, data_point, tree) #計算KNN
for j in knn:
if l[i] != l[j]:
if clusterSim(c[l[i]], c[l[j]], data, alpha) <= 1:
if np_getDistance(data, data_point,l[i], l[j]) <= 1000:
merge(c, l[i], l[j], l)
但是輸出的數值會直接不見=口=
請問我是計算上哪裡有問題，還是code寫錯了，還是array取錯了呢?
若不加np_getDistance的話code跑起來是沒問題
或者還有什麼有效率的方法可以從np.array兩兩計算經緯度的距離呢?
謝謝大家的指教QQ
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 1.174.118.160 (臺灣)
※ 文章網址: https://webptt.com/m.aspx?n=bbs/Python/M.1673321616.A.AB6.html
1F：→ DavisX: 你func的 def 'data'和'data_mid_point' 和裡面的名稱不同01/10 14:07
啊對的，感謝我改正了，但還是跑不出來QQ
※ 編輯: piacere (223.138.15.215 臺灣), 01/10/2023 14:26:04
2F：→ DavisX: 你code不完整 建議你data大概留3點 每個有新算東西就print01/10 14:50
3F：→ DavisX: 你這算球面多點距離 太遠二點誤差太大又要踢掉01/10 14:53
大大我改成功了，目前看起來是正常，我的CODE改成：
def rad(d):
return d * math.pi / 180.0
def getDistance(data_point, ci, cj):
EARTH_REDIUS = 6378.137
radLat1 = rad(data_point[ci,1])
radLat2 = rad(data_point[cj,1])
a = radLat1 - radLat2
b = rad(data_point[ci,0]) - rad(data_point[cj,0])
s = 2 * math.asin(math.sqrt(math.pow(sin(a/2), 2) + cos(radLat1) *
cos(radLat2) * math.pow(sin(b/2), 2)))
s = s * EARTH_REDIUS
print("distance=",s*1000)
return s*1000 #meters
目前又發現一個問題...
我在縮小距離範圍後，卻產生錯誤訊息QQ
https://i.imgur.com/L3DKgGA.jpg
這個是什麼意思呢....資料量變多他也處理不了了QQ
※ 編輯: piacere (1.174.118.160 臺灣), 01/10/2023 15:25:45
※ 編輯: piacere (223.138.15.215 臺灣), 01/10/2023 16:20:50
4F：推 wuyiulin: 你距離是怎麼抓的？話說你要的那個功能通常叫做 Mask01/10 16:36
5F：→ DavisX: 你pop()裡怎還有引數啊, 那pop不就沒效?01/10 17:03
6F：→ lycantrope: 沒有完整code只能算命，dict用pop但沒有key的error01/10 17:07
因為原始code太長不知道該怎麼縮減放上來....orz
※ 編輯: piacere (223.138.15.215 臺灣), 01/10/2023 17:16:07
8F：→ chang1248w: 我記得sklearn裡面有現成的可以用 01/10 21:12
9F：→ lycantrope: 通常是sklearn 把metrics換成經緯度距離 01/10 21:31