二维网格的一维插值 ...

Question

我知道这可能会让人很困惑，所以如果这个解释需要一些编辑请告诉我。

假设我有这种格式的输入数据：

对于给定压力 p_0 --> 参考此压力值的 2x2 温度网格 (T_0)

对于给定压力 p_1 --> 参考此压力值的 2x2 温度网格 (T_1)

p_0 = 0
T_0 = np.array([[1, 4], [3, 2]])

p_1 = 1
T_1 = np.array([[1, 6], [4, 4]])

p = np.array([p_0, p_1])
T = np.array([T_0, T_1])

现在，我得到了一个 2x2 的新压力值网格

p_target = np.array([[0.1, 0.4], [0.3, 0.2]])

我想使用输入数据得到一个 2x2 的插值温度值网格。

我这样做的方法是针对网格的每个点，构建一个插值函数，然后使用它为该网格点获取新的插值温度值：

from scipy.interpolate import interp1d

T_new = np.empty(p_target.shape)

for ix,iy in np.ndindex(p_target.shape):
    f = interp1d(p, T[:,ix,iy])
    T_new[ix,iy] = f(p_target[ix,iy])

T_new

array([[1. , 4.8],
       [3.3, 2.4]])

很容易猜到，这对于大型数组来说相当慢，而且这似乎与 numpy 的处理方式背道而驰。

编辑：我使用 interp1d 也是因为它也允许外推，这是我想保留的一个选项。

Answer 1

我为维度添加了一些参数；从您选择的 n_x = n_y = n_p = 2 来看，依赖关系不是很清楚。

from scipy.interpolate import interp1d, interp2d, dfitpack

n_x = 30
n_y = 40
n_p = 50
T = np.random.random((n_p, n_x, n_y)) * 100
p = np.random.random(n_p)
p[np.argmin(p)] = 0
p[np.argmax(p)] = 1
p_target = np.random.random((n_x, n_y))

T_new = np.empty(p_target.shape)

for ix, iy in np.ndindex(p_target.shape):
    f = interp1d(p, T[:, ix, iy])
    T_new[ix, iy] = f(p_target[ix, iy])

比一句话给你的造型。如果我理解正确的话，你想要 temperature_xy = fun_xy(pressure)，一个单独的函数用于空间网格上的每个坐标。另一种选择可能是将空间分量包含在组合函数中 temperature_xy = fun(pressure, x, y)。对于第二种方法，请查看 scipy.interpolate.griddata.

您可以重新安排第一种方法，使其适用于 interp2d()。为此，第一个维度是压力 x=pressure，第二个维度表示组合空间维度 y=product(x, y)。为了使其表现为 n_x * n_y 压力值的独立插值，我只在创建插值和评估插值时对空间分量使用相同的虚拟值 0、1、2...。因为 interp2d() 的评估通常只适用于网格坐标，我使用 user6655984 提供的方法仅在一组特定的点上评估函数。

def evaluate_interp2d(f, x, y):
    """
    return dfitpack.bispeu(f.tck[0], f.tck[1], f.tck[2], f.tck[3], f.tck[4], x, y)[0]

f2 = interp2d(x=p, y=np.arange(n_x*n_y), z=T.reshape(n_p, n_x*n_y).T)

T_new2 = evaluate_interp2d(f=f2, x=p_target.ravel(), y=np.arange(n_x*n_y))
T_new2 = T_new2.reshape(n_x, n_y)

print(np.allclose(T_new, T_new2))
# True

通过这些设置，我的时间缩短了将近 10x。但是，如果您使用更大的值，例如 n_x=n_y=1000，此 custom interp2d 方法的内存使用量会变得太大，您将使用迭代方法。

# np=50
#    nx*ny      1e2      1e4      1e5      1e6
# interp1d  0.0056s  0.3420s  3.4133s  33.390s
# interp2d  0.0004s  0.0388s  2.0954s  191.66s

有了这些知识，您可以遍历一个大的 1000x1000 网格并按顺序处理 100x100 个片段，然后您将在 3 秒左右结束，而不是 30 秒。

def interpolate2d_flat(p, p_target_flat, T_flat):
    n_p, n_xy = T_flat.shape
    f2 = interp2d(x=p, y=np.arange(n_xy), z=T_flat.T)
    return evaluate_interp2d(f=f2, x=p_target_flat, y=np.arange(n_xy))


n_splits = n_x * n_y // 1000  # So each patch has size n_p*1000, can be changed 

# Flatten and split the spatial dimensions
T_flat_s = np.array_split(T.reshape(n_p, n_x*n_y), n_splits, axis=1)
p_target_flat_s = np.array_split(p_target.ravel(), n_splits, axis=0)

# Loop over the patches
T_new_flat = np.concatenate([interpolate2d_flat(p=p, p_target_flat=ptf, T_flat=Tf)
                             for (ptf, Tf) in zip(p_target_flat_s, T_flat_s)])
T_new2 = T_new_flat.reshape(n_x, n_y)

Answer 2

您可以自己计算插值。这里我假设你有两个以上的 T 值并且 p 不一定是 evenly-spaced。此外，代码假定您有多个 p_target 值，但显然只适用于一个值。

import numpy as np

p_0 = 0
T_0 = np.array([[1., 4.], [3., 2.]])
p_1 = 1
T_1 = np.array([[1., 6.], [4., 4.]])
p = np.array([p_0, p_1])
T = np.array([T_0, T_1])
p_target = np.array([[0.1, 0.4], [0.3, 0.2]])
# Assume you may have several of p_target values
p_target = np.expand_dims(p_target, 0)

# Find the base index for each interpolated value (assume p is sorted)
idx_0 = (np.searchsorted(p, p_target) - 1).clip(0, len(p) - 2)
# And the next index
idx_1 = idx_0 + 1
# Get p values for each interpolated value
a = p[idx_0]
b = p[idx_1]
# Compute interpolation factor
alpha = ((p_target - a) / (b - a)).clip(0, 1)
# Get interpolation values
v_0 = np.take_along_axis(T, idx_0, axis=0)
v_1 = np.take_along_axis(T, idx_1, axis=0)
# Compute interpolation
out = (1 - alpha) * v_0 + alpha * v_1
print(out)
# [[[1.  4.8]
#   [3.3 2.4]]]

编辑：如果您想要线性外推，只需不要剪裁 alpha 值：

alpha = ((p_target - a) / (b - a))

二维网格的一维插值 ...

1D interpolation ... of 2D grids

python

arrays

interpolation

numpy