Python/Scipy 2D插值(非均匀数据)

18
这是对我之前发布的帖子Python/Scipy Interpolation (map_coordinates)的跟进问题。
假设我想在一个二维矩形区域内进行插值。我的变量“z”包含如下所示的数据。每一列都处于一个常数值,然而,数组的每一行可能处于不同的值,如下面的注释所示。
from scipy import interpolate
from numpy import array
import numpy as np
#                                               # 0.0000, 0.1750, 0.8170, 1.0000
z = array([[-2.2818,-2.2818,-0.9309,-0.9309],   # 0.0000, 0.0000, 0.0000, 0.0000
           [-2.2818,-2.2818,-0.9309,-0.9309],   # 0.2620, 0.2784, 0.3379, 0.3526
           [-1.4891,-1.4891,-0.5531,-0.5531],   # 0.6121, 0.6351, 0.7118, 0.7309
           [-1.4891,-1.4891,-0.5531,-0.5531]])  # 1.0000, 1.0000, 1.0000, 1.0000
# Rows, Columns = z.shape

cols = array([0.0000, 0.1750, 0.8170, 1.0000])
rows = array([0.0000, 0.2620, 0.6121, 1.0000])

sp = interpolate.RectBivariateSpline(rows, cols, z, kx=1, ky=1, s=0)

xi = np.array([0.00000, 0.26200, 0.27840, 0.33790, 0.35260, 0.61210, 0.63510,
               0.71180, 0.73090, 1.00000], dtype=np.float)
yi = np.array([0.000, 0.167, 0.815, 1.000], dtype=np.float)
print sp(xi, yi)

作为另一种可视化方式,我知道的值的数组将是:
rows = array([0.0000, 0.2620, 0.2784, 0.3379, 0.3526,
                      0.6121, 0.6351, 0.7118, 0.7309, 1.0000])
#          # 0.0000, 0.1750, 0.8170, 1.0000
z = array([[-2.2818,-2.2818,-0.9309,-0.9309],   # 0.0000
           [-2.2818,      ?,      ?,      ?],   # 0.2620,
           [      ?,-2.2818,      ?,      ?],   # 0.2784
           [      ?,      ?,-0.9309,      ?],   # 0.3379
           [      ?      ,?,      ?,-0.9309],   # 0.3526
           [-1.4891,      ?,      ?,      ?],   # 0.6121
           [      ?,-1.4891,      ?,      ?],   # 0.6351
           [      ?,      ?,-0.5531,      ?],   # 0.7118
           [      ?,      ?,      ?,-0.5531],   # 0.7309
           [-1.4891,-1.4891,-0.5531,-0.5531]])  # 1.0000

我不知道'?'的值,它们需要被插值。我试图用None替换它们,但结果全部变成了'nan'。

编辑:

我认为我需要使用'griddata'或者'interp2'。 'griddata'似乎可以产生我期望的结果,但'interp2'则不行。

from scipy import interpolate
from numpy import array
import numpy as np

z = array([[-2.2818,-2.2818,-0.9309,-0.9309],
           [-2.2818,-2.2818,-0.9309,-0.9309],
           [-1.4891,-1.4891,-0.5531,-0.5531],
           [-1.4891,-1.4891,-0.5531,-0.5531]])

rows = array([0.0000, 0.0000, 0.0000, 0.0000,
              0.2620, 0.2784, 0.3379, 0.3526,
              0.6121, 0.6351, 0.7118, 0.7309,
              1.0000, 1.0000, 1.0000, 1.0000])

cols = array([0.0000, 0.1750, 0.8180, 1.0000,
              0.0000, 0.1750, 0.8180, 1.0000,
              0.0000, 0.1750, 0.8180, 1.0000,
              0.0000, 0.1750, 0.8180, 1.0000])

xi = array([0.0000, 0.2620, 0.2784, 0.3379, 0.3526, 0.6121, 0.6351, 0.7118,
               0.7309, 1.0000], dtype=np.float)
yi = array([0.000, 0.175, 0.818, 1.000], dtype=np.float)

GD = interpolate.griddata((rows, cols), z.ravel(),
                          (xi[None,:], yi[:,None]), method='linear')
I2 = interpolate.interp2d(rows, cols, z, kind='linear')

print GD.reshape(4, 10).T
print '\n'
print I2(xi, yi).reshape(4, 10).T

import matplotlib.pyplot as plt
import numpy.ma as ma

plt.figure()
GD = interpolate.griddata((rows.ravel(), cols.ravel()), z.ravel(),
                          (xi[None,:], yi[:,None]), method='linear')
CS = plt.contour(xi,yi,GD,15,linewidths=0.5,colors='k')
CS = plt.contourf(xi,yi,GD,15,cmap=plt.cm.jet)
plt.colorbar()
plt.scatter(rows,cols,marker='o',c='b',s=5)

plt.figure()
I2 = I2(xi, yi)
CS = plt.contour(xi,yi,I2,15,linewidths=0.5,colors='k')
CS = plt.contourf(xi,yi,I2,15,cmap=plt.cm.jet)
plt.colorbar()
plt.scatter(rows,cols,marker='o',c='b',s=5)
plt.show()

1
interp2在非结构化或非矩形数据上将无法给出预期结果。正如文档http://docs.scipy.org/doc/scipy/reference/interpolate.html中所述,它适用于结构化数据。它甚至进一步说明它是更为严格的规则(结构化)网格。我并不信服,因为`interp2d`的`__doc__`有一个使用矩形(结构化,而不是规则)网格的示例。我详细阐述了我的答案,以尽可能清晰地定义结构化、非结构化和规则网格。 - Paul
1个回答

19

看起来你已经掌握了。

在你的上面的代码示例和之前的问题(链接)中,你拥有结构化数据。这些数据可以使用 RectBivariateSplineinterp2d 进行插值。这意味着你拥有可以在网格上描述的数据(网格上的所有点都有已知值)。网格的 dx 和 dy 不一定相同。(如果所有的 dx 和 dy 相等,你会得到一个规则网格。)

现在,你当前的问题是如果不是所有的点都已知怎么办。这被称为非结构化数据。你只有一个字段中的一些点的选择,不能 necessarily 构造所有顶点都具有已知值的矩形。对于这种类型的数据,你可以使用(像你已经使用的)griddataBivariateSpline 的一种口味。

那么应该选择哪一个呢?

与结构化的 RectBivariateSpline 最接近的类之一是 unstructuredBivariateSpline 中的一个:SmoothBivariateSplineLSQBivariateSpline。如果你想使用样条插值数据,那么就使用这些。这会使你的函数平滑且可微分,但你可能会得到一个超出 Z.max() 或 Z.min() 范围的表面。

由于您设置了ky=1kx=1,并且在structured数据上得到的结果很可能只是线性插值,我个人建议您从RectBivariateSpline样条方案切换到interp2d structured网格插值方案。我知道文档说它适用于regular grids,但是__doc__中的示例仅限于structured而不是规则的。

如果您最终决定切换,我很想知道这些方法之间是否存在显着差异。欢迎使用SciPy。


使用interp2d而不是griddata有什么优势吗?它们似乎基本上做相同的事情。我认为主要区别实际上只是对于griddata,您提供要获取信息的点并获得结果数组。对于interp2d,我认为您会得到一个interp类,然后可以在任何点进行插值。这基本上是唯一的区别吗? - Scott B
我发布了另一个答案,比较了这两种方法。 - Scott B
interp2d 的文档实际上是错误的——它只是在底层使用 BivariateSpline,并且对于 FITPACK 能处理的离散数据也有效。然而,Fitpack 中的二维节点拟合在实践中有点棘手,经常失败…… - pv.
@Paul:来自原始帖子的例子对我有效。插值器(interpolants)有一些微小差异,但这只是因为interp2d是一个张量样条(Tensor Spline),而griddata是基于三角剖分的方法。对于一些数据集,FITPACK 2D例程的效果较差(会产生警告并产生不良的插值器(interpolants)),但这似乎不是这种情况。 - pv.
3
如何使用结构化方法处理NaN值,特别是对于大型数组? - regeirk
显示剩余3条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接