我是numpy / pandas和矢量化计算的新手.我正在做一个数据任务,我有两个数据集.数据集1包含具有经度和纬度的位置列表以及变量A.数据集2还包含具有其经度和纬度的位置列表.对于数据集1中的每个位置,我想计算它到数据集2中所有位置的距离,但我只想得到数据集2中小于变量A值的位数.另请注意数据集非常大,因此我需要使用矢量化操作来加速计算.
例如,我的dataset1可能如下所示：
id lon lat vara
1 20.11 19.88 100
2 20.87 18.65 90
3 18.99 20.75 120
我的数据集2可能如下所示：
placeid lon lat
a 18.75 20.77
b 19.77 22.56
c 20.86 23.76
d 17.55 20.74
然后对于dataset1中的id == 1,我想计算它到数据集2中所有四个点(a,c,c,d)的距离,我希望计算出有多少距离小于相应的距离vara的值.例如,计算的四个距离是90,70,120,110,vara是100.那么该值应该是2.
我已经有了一个矢量化函数来计算两对坐标之间的距离.假设函数(hasrsine(x,y))已正确实现,我有以下代码.
dataset2['count'] = dataset1.apply(lambda x:
haversine(x['lon'],x['lat'],dataset2['lon'], dataset2['lat']).shape[0], axis
= 1)
但是,这会给出总行数,但不会满足我的要求.
有人能够指出我如何使代码工作？
解决方法:
如果您可以将坐标投影到局部投影(例如UTM),这对于pyproj非常直接并且通常比lon / lat更有利于测量,那么使用scipy.spatial有更多更快的方法. df [‘something’] = df.apply(…)和np.vectorize()都没有真正矢量化,在引擎盖下,它们使用循环.
ds1
id lon lat vara
0 1 20.11 19.88 100
1 2 20.87 18.65 90
2 3 18.99 20.75 120
ds2
placeid lon lat
0 a 18.75 20.77
1 b 19.77 22.56
2 c 20.86 23.76
3 d 17.55 20.74
from scipy.spatial import distance
# gey coordinates of each set of points as numpy array
coords_a = ds1.values[:,(1,2)]
coords_b = ds2.values[:, (1,2)]
coords_a
#out: array([[ 20.11, 19.88],
# [ 20.87, 18.65],
# [ 18.99, 20.75]])
distances = distance.cdist(coords_a, coords_b)
#out: array([[ 1.62533074, 2.70148108, 3.95182236, 2.70059253],
# [ 2.99813275, 4.06178532, 5.11000978, 3.92307278],
# [ 0.24083189, 1.97091349, 3.54358575, 1.44003472]])
距离实际上是每对点之间的距离. coords_a.shape是(3,2),coords_b.shape是(4,2),所以结果是(3,4). np.distance的默认度量标准是eculidean,但也有其他指标.
为了这个例子,让我们假设vara是：
vara = np.array([2,4.5,2])
(而不是100 90 120).我们需要确定第一行中距离中哪个值小于2,第二行中哪个值小于4.5,…,解决此问题的一种方法是从相应行中减去vara中的每个值(注意我们必须调整vara的大小) ：
vara.resize(3,1)
res = res - vara
#out: array([[-0.37466926, 0.70148108, 1.95182236, 0.70059253],
# [-1.50186725, -0.43821468, 0.61000978, -0.57692722],
# [-1.75916811, -0.02908651, 1.54358575, -0.55996528]])
然后将正值设置为零并将负值设为正值将为我们提供最终数组：
res[res>0] = 0
res = np.absolute(res)
#out: array([[ 0.37466926, 0. , 0. , 0. ],
# [ 1.50186725, 0.43821468, 0. , 0.57692722],
# [ 1.75916811, 0.02908651, 0. , 0.55996528]])
现在,总结每一行：
sum_ = res.sum(axis=1)
#out: array([ 0.37466926, 2.51700915, 2.34821989])
并计算每行中的项目：
count = np.count_nonzero(res, axis=1)
#out: array([1, 3, 3])
这是一个完全矢量化(自定义)的解决方案,您可以根据自己的喜好进行调整,并且应该适应任何级别的复杂性.另一个解决方案是cKDTree.代码来自文档.将它用于你的问题应该相当容易,但如果你需要帮助,请不要犹豫.
x, y = np.mgrid[0:4, 0:4]
points = zip(x.ravel(), y.ravel())
tree = spatial.cKDTree(points)
tree.query_ball_point([2, 0], 1)
[4, 8, 9, 12]
query_ball_point()查找点(x)的距离r内的所有点,并且速度惊人.
最后一点注意：不要将这些算法与lon / lat输入一起使用,特别是如果您感兴趣的区域远离赤道,因为错误会变得很大.
更新：
要投影坐标,您需要将wgs84(lon / lat)转换为适当的UTM.要找出你应该投射哪个区域使用epsg.io.
lon = -122.67598
lat = 45.52168
wgs84 = "+init=epsg:4326"
epsg3740 = "+init=epsg:3740"
Proj_to_epsg3740 = pyproj.Proj(epsg3740)
Proj_to_epsg3740(lon,lat)
# out: (525304.9265963673, 5040956.147893889)
你可以做df.apply()并使用Proj_to _…来投射df.
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 [email protected] 举报，一经查实，本站将立刻删除。
