python - Numpy 网格数据插值到一定半径

Question

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

是否可以告诉 numpy/scipy:如果您距离现有测量值太远，则不要进行插值？产生 NODATA 值？ ideal = 网格数据(.., .., .., 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

编辑:好的，重新阅读并注意到输入实际上是二维的。示例已修改。

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

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

类似于

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

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