python矢量化运算_python – 用于计算许多距离的矢量化

如果您可以将坐标投影到局部投影(例如

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.