python - 在 numpy 数组 ((x, y, z)...) 中搜索与最近的 x, y 匹配的 z

Question

我有一个非常大的数组，类似于以下格式的高程数据:

triplets = ((x0, y0, z0), 
            (x1, y1, z1), 
            ... ,
            (xn, yn, zn))

其中 x、y、z 都是以米为单位的 float 。您可以创建与此格式匹配的合适测试数据:

x = np.arange(20, 40, dtype=np.float64)
y = np.arange(30, 50, dtype=np.float64)
z = np.random.random(20) * 25.0
triplets = np.hstack((x, y, z)).reshape((len(x),3))

我希望能够高效地找到给定 (x,y) 对的相应 z 值。到目前为止，我的研究引出了更多问题。这是我得到的:

1. 遍历所有三元组:

   query = (a, b) # where a, b are the x and y coordinates we're looking for
   for i in triplets:
     if i[0] == query[0] and i[1] == query[1]:
       result = i[2]
   

   缺点:速度慢； a, b必须存在，这是 float 比较的问题。

2. 使用 scipy.spatial.cKDTree找到最近的:

   points = triplets[:,0:2] # drops the z column
   tree = cKDTree(points)
   idx = tree.query((a, b))[1] # this returns a tuple, we want the index
   query = tree.data[idx]
   result = triplets[idx, 2]
   

   缺点:返回最近点而不是插值。

3. 使用 interp2d根据评论:

   f = interp2d(x, y, z)
   result = f(a, b)
   

   缺点:不适用于大型数据集。在实际数据上运行时，我得到 OverflowError: Too many data points to interpolate。 (我的真实数据是1100万点左右。)

所以问题是:有没有我忽略的直接方法？有没有办法减少上述的弊端？

Answer 1

如果您想对结果进行插值，而不是仅仅找到最近邻居的 z 值，我会考虑执行如下操作:

1. 使用 k-d 树根据它们的 (x, y) 坐标划分数据点
2. 对于给定的(xi, yi)点进行插值，找到它的k个最近的邻居
3. 取它们的 z 值的平均值，根据它们与 (xi, yi) 的距离加权

代码可能是这样的:

import numpy as np
from scipy.spatial import cKDTree

# some fake (x, y, z) data
XY = np.random.rand(10000, 2) - 0.5
Z = np.exp(-((XY ** 2).sum(1) / 0.1) ** 2)

# construct a k-d tree from the (x, y) coordinates
tree = cKDTree(XY)

# a random point to query
xy = np.random.rand(2) - 0.5

# find the k nearest neighbours (say, k=3)
distances, indices = tree.query(xy, k=3)

# the z-values for the k nearest neighbours of xy
z_vals = Z[indices]

# take the average of these z-values, weighted by 1 / distance from xy
dw_avg = np.average(z_vals, weights=(1. / distances))

值得研究一下 k 的值，即要取平均值的最近邻居的数量。这本质上是 kernel density estimation 的粗略形式。 ，其中 k 的值控制您对 z 值的基础分布施加的“平滑度”程度。 k 越大，平滑度越高。

同样，您可能想根据点与 (xi, yi) 的距离对点的贡献进行加权，这取决于您如何看待z 随着 x, y 距离的增加而减小。例如，您可能希望通过 (1 / distances ** 2) 来加权而不是 (1 / distances) .

在性能方面，constructing and searching k-d trees are both very efficient .请记住，您只需为数据集构建一次树，如有必要，您可以通过将 (N, 2) 数组传递给 tree.query() 来一次查询多个点。 .

用于近似最近邻搜索的工具，例如 FLANN ，可能会更快，但这些通常在数据维度非常高的情况下更有帮助。