给定1D输入时,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仅包含第二轴的一个不同值,如上面的测试样本.该bisplrep函数内部interp2d认为它作为一个无效的输出:这可以被认为是可接受的行为: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中没有？

或者,换句话说:当2D插值仅沿一个方向输入(情况1和2失败)时,我们是否希望它返回错误？(案例1和2应该返回一些东西,即使是不可预测的).

Answer 1

如果您的输入数据沿着坐标轴而不是一般方向定向,我原本会向您展示它对二维插值有多大差异,但事实证明结果会比我预期的更加混乱.我尝试在插值矩形网格上使用随机数据集,并将其与相同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͈͇̼͖͖̭̙͈͇̼͖͖̭̙z̻̉ͬͪ̑ͭͨ͊ǟ̼̣̬̗̖ͥl̫̣͔͓̟͛͊̏ͨ͗g̻͇͈͚̟̻͛ͫ͛̅͋͒o͈͓̥̙̫͚̾的威胁.

所以我们也这样做griddata:

感谢The Powerpuff Girls ,这一天得救了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数据集中尝试二维插值可能毫无意义(或至少非常不明确),但如果您的数据沿着x,y平面的大致方向,则至少算法不应该破坏.