给定一维输入时，scipy interp2d/bisplrep 意外输出

Question

在使用 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: Invalid inputs

结果我给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 在案例 2 中有效，但在案例 1 中无效？

或者，换句话说：当二维插值仅沿一个方向输入时（案例 1 和 2 失败），我们是否希望它返回错误？ （情况 1 和 2 应该返回一些东西，即使是不可预测的）。

Answer 1

我最初打算向您展示如果您的输入数据沿坐标轴而不是某个一般方向定向，它对 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͔̖̰͕̞͖͇ͣ̈̒ͦü͇̹̞̳ͭ̊̓̈m̥̠͈̣̆̐ͦ̚m̻͑͒̔̓ͦ̇oği͉̓̓ͭ̒͛n̹̙̥̩̥̯̭ͤͤͤ̄g͈͇̼͖͖̭̙͈͇̼͖͖̭̙̻̐̐̉ͬͪ̑ͭͨ͊ͣ̐ͣṉ̟͖͙̆͋ͣ̐ͣṉ̟͖͙̆͋p͈͓̥̙̫͚̾的威胁。

让我们对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 插值可能毫无意义（或至少定义非常不明确），但如果您的数据沿着x,y 平面的大致方向，至少该算法不应该中断。