Numpy griddata插值到一定半径

Question

我正在使用griddata() 插入我的（不规则的）二维深度测量； x,y,depth。该方法做得很好 - 但它在整个网格上进行插值，它可以找到相反的点。我不想要那种行为。我想在现有测量值周围进行插值，比如在一定半径范围内。

是否可以告诉 numpy/scipy：如果您离现有测量太远，请不要插值？导致 NODATA 值？理想=griddata(.., .., .., radius=5.0)

编辑示例： 在下图中；黑点是测量值。蓝色阴影是 numpy 插值的单元格。以绿色标记的区域实际上是图片的一部分，但被 numpy 视为 NODATA（因为中间没有点）。现在，红色区域已被插值，但我想摆脱它们。有任何想法吗？

Answer 1

好的，酷。我认为griddata() 没有内置选项可以满足您的需求，因此您需要自己编写。

这归结为计算N 输入数据点和M 插值点之间的距离。这很简单，但如果你有很多点，在 ~O(M*N) 时可能会很慢。但这里有一个例子，计算每个插值点到所有N 数据点的距离。如果半径内的数据点数至少为neighbors，则保留该值。否则就是写入NODATA的值。

neighbors 为 4，因为 griddata() 将使用双线性插值，这需要在每个维度 (2*2 = 4) 中限定插值的点。

#invec - input points Nx2 numpy array
#mvec - interpolation points Mx2 numpy array

#just some random points for example
N=100
invec = 10*np.random.random([N,2])

M=50
mvec = 10*np.random.random([M,2])

# --- here you would put your griddata() call, returning interpolated_values
interpolated_values = np.zeros(M)
NODATA=np.nan

radius = 5.0
neighbors = 4

for m in range(M):
    data_in_radius = np.sqrt(np.sum( (invec - mvec[m])**2, axis=1)) <= radius

    if np.sum(data_in_radius) < neighbors :
        interpolated_values[m] = NODATA

编辑： 好的，重新阅读并注意到输入确实是 2D 的。示例已修改。

作为补充说明，如果您首先构建从每个点 mvec[m] 到相关数据点的子集的粗略映射，这可以大大加快速度。 循环中成本最高的步骤将从

np.sqrt(np.sum( (invec - mvec[m])**2, axis=1))

类似

np.sqrt(np.sum( (invec[subset[m]] - mvec[m])**2, axis=1))

有很多方法可以做到这一点，例如使用四叉树、散列函数或二维索引。但这是否会带来性能优势取决于应用程序、数据的结构等。