作者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/cn.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