如何解决2D网格的1D插值...
我知道这可能会造成混乱，所以请告诉我该解释是否需要修改。
假设我输入的数据格式如下：
对于给定压力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]])
很容易猜到，这对于大型数组而言相当慢，而且似乎与麻木的处理方式背道而驰。
编辑：我使用interp1d也是因为它也允许外推，这是我想保留的一个选项。
解决方法
您可以自己计算插值。在这里，我假设您有两个以上的T值，并且p不一定是均匀间隔的。另外，该代码假定您具有多个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))
我为尺寸添加了一些参数；根据您选择的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,y):
"""https://stackoverflow.com/a/47233198/7570817"""
return dfitpack.bispeu(f.tck[0],f.tck[1],f.tck[2],f.tck[3],f.tck[4],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，则此自定义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(),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)
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。