二维插值问题

Question

我有一些数据，x轴和y轴上都有数值，输出结果在z轴上。

比如说：

y = 10
x = [1,2,3,4,5,6]
z = [2.3,3.4,5.6,7.8,9.6,11.2]

y = 20 
x = [1,2,3,4,5,6]
z = [4.3,5.4,7.6,9.8,11.6,13.2]

y = 30 
x = [1,2,3,4,5,6]
z = [6.3,7.4,8.6,10.8,13.6,15.2]

我想知道当y等于15，x等于3.5时，z的值是多少。

我试着用scipy这个工具，但我对它还很陌生。

非常感谢你的帮助。

vibhor

Answer 1

我觉得可以简单地取周围数值的平均值。比如说，如果你需要的坐标是X=3.5和Y=15（也就是(3.5,15)），你可以计算(3,10)、(3,20)、(4,10)和(4,20)这几个点的平均值。因为我不知道你正在处理的数据是什么，所以不太确定距离的准确性是否重要。如果不重要的话，你就可以直接用平均值；如果重要的话，你可能需要考虑一些反距离加权的方法。

Answer 2

假设有这样的情况（不是Python代码，因为第二次赋值会覆盖第一次的结果，当然这是显而易见的；-）：

y = 10
x = [1,2,3,4,5,6]
z = [2.3,3.4,5.6,7.8,9.6,11.2]

y = 20 
x = [1,2,3,4,5,6]
z = [4.3,5.4,7.6,9.8,11.6,13.2]

你问：“当y=15，x=3.5时，如何找到z的值？”

因为你在看一个点，它在给定的“网格”中，x和y的距离是一样的，所以你只需要在网格值之间找一个中间值（如果值的距离不一样，你就要找一个按比例的中间值，稍后会讲到）。所以对于y=10，x为3和4时，z的值分别是5.6和7.8，那么对于x=3.5，你可以估算它们的中间值，结果是6.7；同样，对于y=20，你可以估算7.6和9.8之间的中间值，也就是8.7。最后，因为你有y=15，所以6.7和8.7之间的中间值就是你最终估算的z值：7.7。

假设你有y=13和x=3.8。那对于x，你需要取80%的值，也就是：

• 对于y=10，0.2*5.6+0.8*7.8 -> 7.36

• 对于y=20，0.2*7.6+0.8*9.8 -> 9.46

现在你想要z值在这两个值之间的30%的位置，计算方式是0.3*7.36 + 0.7*9.46 -> 8.83，这就是z。

这就是线性插值，其实非常简单。你想手动计算，还是找一些可以帮你计算的工具（比如用numpy数组作为“网格”）？即使是后者，我希望这个“手动”的解释（用最基础的算术术语展示你在做什么）能帮助你理解你在做的事情...;-）。

当然，还有更高级的插值方法——你需要那些吗，还是线性插值就足够满足你的需求？

Answer 3

scipy.interpolate.bisplrep

参考资料： http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate/bisplrep.html

import scipy
import math
import numpy
from scipy import interpolate


x= [1,2,3,4,5,6]
y= [10,20,30]

Y = numpy.array([[i]*len(x) for i in y])
X = numpy.array([x for i in y])
Z = numpy.array([[2.3,3.4,5.6,7.8,9.6,11.2],
                 [4.3,5.4,7.6,9.8,11.6,13.2],
                 [6.3,7.4,8.6,10.8,13.6,15.2]]) 

tck = interpolate.bisplrep(X,Y,Z)
print interpolate.bisplev(3.5,15,tck) 


7.84921875

编辑：

上面的解决方案并不能完美拟合数据。请检查

print interpolate.bisplev(x,y,tck)

[[  2.2531746    4.2531746    6.39603175]
 [  3.54126984   5.54126984   7.11269841]
 [  5.5031746    7.5031746    8.78888889]
 [  7.71111111   9.71111111  10.9968254 ]
 [  9.73730159  11.73730159  13.30873016]
 [ 11.15396825  13.15396825  15.2968254 ]]

为了克服这个问题，可以在x方向使用5次多项式，在y方向使用2次多项式进行插值。

tck = interpolate.bisplrep(X,Y,Z,kx=5,ky=2)
print interpolate.bisplev(x,y,tck) 

[[  2.3   4.3   6.3]
 [  3.4   5.4   7.4]
 [  5.6   7.6   8.6]
 [  7.8   9.8  10.8]
 [  9.6  11.6  13.6]
 [ 11.2  13.2  15.2]]

这样可以得到

print interpolate.bisplev(3.5,15,tck)

7.88671875

绘图：
参考资料 http://matplotlib.sourceforge.net/examples/mplot3d/surface3d_demo.html

fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z,rstride=1, cstride=1, cmap=cm.jet)
plt.show()