当给出1D输入时，scipy interp2d / bisplrep的输出结果出乎意料。

Question

当给出1D输入时，scipy interp2d / bisplrep的输出结果出乎意料。

5

我在使用scipy的interp2d函数时遇到了无效输入错误。原来问题出在bisplrep函数中，如下所示：

import numpy as np
from scipy import interpolate

# Case 1
x = np.linspace(0,1)
y = np.zeros_like(x)
z = np.ones_like(x)

tck = interpolate.bisplrep(x,y,z)  # or interp2d

返回结果：

返回：ValueError: 无效输入

事实证明，我给interp2d的测试数据在第二轴上只包含一个不同的值，就像上面的测试样例一样。 interp2d内部的bisplrep函数将其视为无效输出：这可能被认为是可接受的行为：interp2d和bisplrep需要一个2D网格，而我只提供它们沿着一条线的值。

顺便说一句，我发现错误消息非常不清楚。可以在interp2d中包含一个测试来处理这种情况：类似以下内容的东西

if len(np.unique(x))==1 or len(np.unique(y))==1: 
    ValueError ("Can't build 2D splines if x or y values are all the same")

可能足以检测到这种无效的输入，并引发更明确的错误消息，甚至直接调用更合适的interp1d函数（在这里完美运行）

我以为我已经正确地理解了问题。然而，请考虑以下代码示例：

# Case 2
x = np.linspace(0,1)
y = x
z = np.ones_like(x)

tck = interpolate.bisplrep(x,y,z)

在这种情况下，y与x成比例，我也将数据沿一条线输入bisplrep。但是令人惊讶的是，在这种情况下，bisplrep能够计算出2D样条插值。我绘制了它：

# Plot
def plot_0to1(tck):
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D

    X = np.linspace(0,1,10)
    Y = np.linspace(0,1,10)
    Z = interpolate.bisplev(X,Y,tck)

    X,Y = np.meshgrid(X,Y)

    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot_surface(X, Y, Z,rstride=1, cstride=1, cmap=cm.coolwarm,
                    linewidth=0, antialiased=False)
    plt.show()

plot_0to1(tck)

结果如下所示：

在这里，bisplrep似乎用0填充了缺口，当我扩展下面的图时更好地显示出来：

关于是否期望添加0，我的真正问题是：为什么bisplrep在Case 1中无法工作，但在Case 2中可以工作？ 换句话说：当仅沿一个方向输入2D插值时（Case 1和2失败），我们希望它返回错误吗？还是不希望？（Case 1和2应该返回一些东西，即使是不可预测的）。

- erwanp

哪个函数或子调用返回了 ValueError？我在 bisplrep 中没有看到这样的 raise。这个函数是 FITPACK 的前端，一个 FORTRAN 库。像这样的库并非以用户友好著称。它们是由专业人员为自己和其他专业人员编写的。 - hpaulj

我相信这样做是可以的。考虑双线性插值:您独立地沿x和y进行插值。这意味着输入数据必须对于插值有效。但是，如果在您的输入数据中有一个轴是恒定的，则无法沿该方向插值。在这种情况下，您只是缺少信息。 - Andras Deak -- Слава Україні

1

为了让我的观点更清晰：坐标轴的方向在通常的二维插值中会产生非常强的偏差。考虑双线性情况：插值函数的偏导数不连续性由笛卡尔坐标轴的方向确定。因此，数据是沿着坐标轴还是以某种不同的配置排列都非常重要。这就是为什么使用倾斜的输入行（使用interp2d(x,x,z)）在计算上要好得多的原因：即使几何上仍然使用一条直线，但有许多x和y值。 - Andras Deak -- Слава Україні

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Andras Deak -- Слава Україні · Accepted Answer

我原本想向您展示，如果您的输入数据沿着坐标轴而不是在某个一般方向上定位，2D插值会产生多大的差异，但结果比我预期的要混乱得多。我尝试在插值矩形网格上使用随机数据集进行比较，并将其与相同的 x 和 y 坐标旋转45度进行插值的情况进行比较。结果很糟糕。

然后我尝试使用更平滑的数据集进行比较：结果发现scipy.interpolate.interp2d存在很多问题。因此，我的底线是“使用scipy.interpolate.griddata”。

为了说明目的，这是我的（非常混乱的）代码：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.cm as cm

n = 10                            # rough number of points
dom = np.linspace(-2,2,n+1)       # 1d input grid
x1,y1 = np.meshgrid(dom,dom)      # 2d input grid
z = np.random.rand(*x1.shape)     # ill-conditioned sample
#z = np.cos(x1)*np.sin(y1)        # smooth sample

# first interpolator with interp2d:
fun1 = interp.interp2d(x1,y1,z,kind='linear')

# construct twice finer plotting and interpolating mesh
plotdom = np.linspace(-1,1,2*n+1)             # for interpolation and plotting
plotx1,ploty1 = np.meshgrid(plotdom,plotdom)
plotz1 = fun1(plotdom,plotdom)                # interpolated points


# construct 45-degree rotated input and interpolating meshes
rotmat = np.array([[1,-1],[1,1]])/np.sqrt(2)           # 45-degree rotation
x2,y2 = rotmat.dot(np.vstack([x1.ravel(),y1.ravel()])) # rotate input mesh
plotx2,ploty2 = rotmat.dot(np.vstack([plotx1.ravel(),ploty1.ravel()])) # rotate plotting/interp mesh

# interpolate on rotated mesh with interp2d
# (reverse rotate by using plotx1, ploty1 later!)
fun2 = interp.interp2d(x2,y2,z.ravel(),kind='linear')

# I had to generate the rotated points element-by-element
# since fun2() accepts only rectangular meshes as input
plotz2 = np.array([fun2(xx,yy) for (xx,yy) in zip(plotx2.ravel(),ploty2.ravel())])

# try interpolating with griddata
plotz3 = interp.griddata(np.array([x1.ravel(),y1.ravel()]).T,z.ravel(),np.array([plotx1.ravel(),ploty1.ravel()]).T,method='linear')
plotz4 = interp.griddata(np.array([x2,y2]).T,z.ravel(),np.array([plotx2,ploty2]).T,method='linear')


# function to plot a surface
def myplot(X,Y,Z):
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot_surface(X, Y, Z,rstride=1, cstride=1,
                    linewidth=0, antialiased=False,cmap=cm.coolwarm)
    plt.show()


# plot interp2d versions
myplot(plotx1,ploty1,plotz1)                    # Cartesian meshes
myplot(plotx1,ploty1,plotz2.reshape(2*n+1,-1))  # rotated meshes

# plot griddata versions
myplot(plotx1,ploty1,plotz3.reshape(2*n+1,-1))  # Cartesian meshes
myplot(plotx1,ploty1,plotz4.reshape(2*n+1,-1))  # rotated meshes

这里是结果展示。使用随机输入 z 数据和 interp2d，笛卡尔（左）与旋转插值（右）：

请注意右侧的可怕比例尺，输入点位于 0 到 1 之间。即使母亲也不会认出数据集。请注意，在评估旋转数据集时会有运行时警告，因此我们被警告一切都是垃圾。
现在让我们使用 griddata 做同样的事情：

我们应该注意到这些图像更接近，并且它们似乎比 interp2d 的输出更有意义。例如，请注意第一个图像中的超调现象。
这些伪影总是在输入数据点之间出现。由于它仍然是插值，因此插值函数必须复制输入点，但线性插值函数在数据点之间超调是相当奇怪的。很明显，griddata 不会遭受这个问题。
考虑一个更清晰的示例：其他一组 z 值，它们是平滑且确定的。与 interp2d 相关的表面：

救命啊！请召唤插值警察！笛卡尔输入案例中已经出现了不可解释的（至少对我来说是这样）虚假特征，而旋转输入案例则可能会出现s͔̖̰͕̞͖͇ͣ́̈̒ͦ̀̀ü͇̹̞̳ͭ̊̓̎̈m̥̠͈̣̆̐ͦ̚m̻͑͒̔̓ͦ̇oͣ̐ͣṉ̟͖͙̆͋i͉̓̓ͭ̒͛n̹̙̥̩̥̯̭ͤͤͤ̄g͈͇̼͖͖̭̙ ̐z̻̉ͬͪ̑ͭͨ͊ä̼̣̬̗̖́̄ͥl̫̣͔͓̟͛͊̏ͨ͗̎g̻͇͈͚̟̻͛ͫ͛̅͋͒o͈͓̱̥̙̫͚̾͂的威胁。

因此，让我们使用griddata进行相同的操作：

多亏了~~超能力女孩~~scipy.interpolate.griddata，我们挽救了这一天。作业：使用cubic插值检查是否相同。

顺便说一下，在help(interp.interp2d)中可以得到你最初问题的非常简短的答案：

 |  Notes
 |  -----
 |  The minimum number of data points required along the interpolation
 |  axis is ``(k+1)**2``, with k=1 for linear, k=3 for cubic and k=5 for
 |  quintic interpolation.

线性插值需要沿插值轴至少有4个点，也就是至少有4个独特的x和y值才能得到有意义的结果。请检查以下内容：

nvals = 3  # -> RuntimeWarning
x = np.linspace(0,1,10)
y = np.random.randint(low=0,high=nvals,size=x.shape)
z = x
interp.interp2d(x,y,z)

nvals = 4  # -> no problem here
x = np.linspace(0,1,10)
y = np.random.randint(low=0,high=nvals,size=x.shape)
z = x
interp.interp2d(x,y,z)

当然，这一切都与您的问题有关：如果您的几何1D数据集沿着笛卡尔坐标轴之一，则会产生巨大差异，或者如果它以一般方式存在，使得坐标值假定各种不同的值。尝试从几何1D数据集进行2D插值可能是没有意义的（或者至少非常模糊），但是如果您的数据沿着平面的一般方向，则算法至少不应该中断。