[問題] 如何使用np.array的經緯度計算兩兩的距離

看板Python作者piacere (Beol)時間2年前 (2023/01/10 11:33)推噓1(1推 0噓 8→)

留言9則, 4人參與討論串1/1

各位先進們好我是一名基礎不太好的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://www.ptt.cc/bbs/Python/M.1673321616.A.AB6.html

→

DavisX

01/10 14:07, 2年前 , 1^F

01/10 14:07, 1^F

啊對的，感謝我改正了，但還是跑不出來QQ ※ 編輯: piacere (223.138.15.215 臺灣), 01/10/2023 14:26:04

→

DavisX

01/10 14:50, 2年前 , 2^F

01/10 14:50, 2^F

→

DavisX

01/10 14:53, 2年前 , 3^F

01/10 14:53, 3^F

大大我改成功了，目前看起來是正常，我的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