给定1D输入时scipy interp2d/bisplrep意外输出 [英] scipy interp2d/bisplrep unexpected output when given 1D input

问题描述

使用scipy interp2d函数时，我一直遇到无效的输入错误.事实证明，问题出在bisplrep函数，如下所示:

I've been having invalid input errors when working with scipy interp2d function. It turns out the problem comes from the bisplrep function, as showed here:

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网格，而我只给它们赋值沿一行.

It turned out the test data I was giving interp2d contained only one distinct value for the 2nd axis, as in the test sample above. The bisplrep function inside interp2d considers it as an invalid output: This may be considered as an acceptable behaviour: interp2d & bisplrep expect a 2D grid, and I'm only giving them values along one line.

另一方面，我发现错误消息还不清楚.可以在interp2d中包含一个测试来处理此类情况:类似

On a side note, I find the error message quite unclear. One could include a test in interp2d to deal with such cases: something along the lines of

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函数(在这里非常有效)

may be enough to detect this kind of invalid input, and raise a more explicit error message, or even directly call the more appropriate interp1d function (which works perfectly here)

我以为我已经正确理解了这个问题.但是，请考虑以下代码示例:

I thought I had correctly understood the problem. However, consider the following code sample:

# 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样条插值.我绘制了它:

In that case, y being proportional to x, I'm also feeding bisplrep with data along one line. But, surprisingly, bisplrep is able to compute a 2D spline interpolation in that case. I plotted it:

# 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填补了空白，当我扩展下面的图时更好地显示了这一点:

where bisplrep seems to fill the gaps with 0's, as better showed when I extend the plot below:

关于是否期望加0，我真正的问题是:为什么bisplrep在案例2中起作用而在案例1中不起作用?

Regarding of whether adding 0 is expected, my real question is: why does bisplrep work in Case 2 but not in Case 1?

或者，换句话说:当仅沿一个方向输入2D插补时(是否出现情况1和2)，我们是否希望它返回错误? (案例1和案例2应该返回某些内容，即使是不可预测的).

Or, in other words: do we want it to return an error when 2D interpolation is fed with input along one direction only (Case 1 & 2 fail), or not? (Case 1 & 2 should return something, even if unpredicted).

推荐答案

我本来将向您展示如果输入数据沿坐标轴而不是沿某个通用方向定向，则它对2d插值有多大影响，但结果却比我预期的还要混乱.我尝试在插值的矩形网格上使用随机数据集，并将其与相同的x和y坐标旋转45度进行插值的情况进行比较.结果很糟糕.

I was originally going to show you how much of a difference it makes for 2d interpolation if your input data are oriented along the coordinate axes rather than in some general direction, but it turns out that the result would be even messier than I had anticipated. I tried using a random dataset over an interpolated rectangular mesh, and comparing that to a case where the same x and y coordinates were rotated by 45 degrees for interpolation. The result was abysmal.

然后我尝试与更平滑的数据集进行比较:事实证明scipy.interpolate.interp2d存在很多问题.所以我的底线将是"use scipy.interpolate.griddata".

I then tried doing a comparison with a smoother dataset: turns out scipy.interpolate.interp2d has quite a few issues. So my bottom line will be "use scipy.interpolate.griddata".

出于指导目的，这是我的代码(很乱):

For instructive purposes, here's my (quite messy) code:

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，笛卡尔(左)与旋转插值(右):

So here's a gallery of the results. Using random input z data, and interp2d, Cartesian (left) vs rotated interpolation (right):

请注意右侧的可怕范围，请注意输入点在0和1之间.即使是它的母亲也无法识别数据集.请注意，在评估旋转数据集的过程中会出现运行时警告，因此我们被警告说这都是废话.

Note the horrible scale on the right side, noting that the input points are between 0 and 1. Even its mother wouldn't recognize the data set. Note that there are runtime warnings during the evaluation of the rotated data set, so we're being warned that it's all crap.

现在让我们对griddata做同样的事情:

Now let's do the same with griddata:

我们应该注意，这些数字彼此之间非常接近，并且它们似乎比interp2d的输出更有意义.例如，请注意第一个数字的比例过冲.

We should note that these figures are much closer to each other, and they seem to make way more sense than the output of interp2d. For instance, note the overshoot in the scale of the very first figure.

这些伪像总是出现在输入数据点之间.由于仍然是插值，因此必须通过插值函数来再现输入点，但是线性插值函数在数据点之间过冲是很奇怪的.显然，griddata不受此问题的困扰.

These artifacts always arise between input data points. Since it's still interpolation, the input points have to be reproduced by the interpolating function, but it's pretty weird that a linear interpolating function overshoots between data points. It's clear that griddata doesn't suffer from this issue.

考虑一个更明确的情况:另一组z值是平滑且确定的. interp2d:

Consider an even more clear case: the other set of z values, which are smooth and deterministic. The surfaces with interp2d:

帮助！致电插值警察！笛卡尔输入箱中已经存在莫名其妙的伪造特征(至少对我而言)，旋转的输入箱构成了s͔̖̰͕̞͖͇ͣ́̈̒ͦ̀̀ü̈̒ͦ̀̀س͇̹̞̳ͭm̥̠͈̣̆̐ͦ̚m̻͑͒̔̓ͦ̇mͣ̐ͣṉ̟͖͙̆͋oͣ̐ͣṉ̟͖͙̆͋i͉̓̓ͭ̒͛n̹̙̥̩̥̯̭ͤͤͤ̄g͈͇̼͖͖̭̙̐z̻̉ͬͪ̑ͭͨ͊ä̼̣̬̗̖́̄ͥl̫̣͔͓̟͛͊̏ͨ͗̎g̻͇͈͚̟̻͛ͫ͛̅͋͒o͈͓̱̥̙̫͚̾͂的威胁.

HELP! Call the interpolation police! Already the Cartesian input case has inexplicable (well, at least by me) spurious features in it, and the rotated input case poses the threat of s͔̖̰͕̞͖͇ͣ́̈̒ͦ̀̀ü͇̹̞̳ͭ̊̓̎̈m̥̠͈̣̆̐ͦ̚m̻͑͒̔̓ͦ̇oͣ̐ͣṉ̟͖͙̆͋i͉̓̓ͭ̒͛n̹̙̥̩̥̯̭ͤͤͤ̄g͈͇̼͖͖̭̙ ̐z̻̉ͬͪ̑ͭͨ͊ä̼̣̬̗̖́̄ͥl̫̣͔͓̟͛͊̏ͨ͗̎g̻͇͈͚̟̻͛ͫ͛̅͋͒o͈͓̱̥̙̫͚̾͂.

所以我们对griddata做同样的事情:

So let's do the same with griddata:

感谢飞天小女警 scipy.interpolate.griddata，这一天得以保存.作业:使用cubic插值进行检查.

The day is saved, thanks to The Powerpuff Girls scipy.interpolate.griddata. Homework: check the same with cubic interpolation.

顺便说一句，您对原始问题的简短回答是在help(interp.interp2d)中:

By the way, a very short answer to your original question is in 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值才能获得有意义的结果.检查这些:

For linear interpolation you need at least 4 points along the interpolation axis, i.e. at least 4 unique x and y values have to be present to get a meaningful result. Check these:

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平面的大致方向，则至少算法不会中断.

And of course this all ties in to you question like this: it makes a huge difference if your geometrically 1d data set is along one of the Cartesian axes, or if it's in a general way such that the coordinate values assume various different values. It's probably meaningless (or at least very ill-defined) to try 2d interpolation from a geometrically 1d data set, but at least the algorithm shouldn't break if your data are along a general direction of the x,y plane.

这篇关于给定1D输入时scipy interp2d/bisplrep意外输出的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！