Scipy interp2d / bisplrep unexpected exit given 1D input

Question

I had the wrong input errors while working with the scipy interp2d function. It turns out the problem comes from the bisplrep function, as shown 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 

Returns: ValueError: Invalid inputs

It turned out that the test data that I gave interp2d contained only one excellent value for the 2nd axis, as in the above test sample. The bisplrep function inside interp2d considers it unacceptable: This can be considered as acceptable behavior: interp2d and bisplrep expect a 2D mesh, and I give them values along the same line .

The error message, on the other hand, is rather obscure. It would be possible to conduct a test in interp2d to examine such cases: something like strings

 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") 

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

---

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

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

In this case, y proportional to x , I also span bisplrep data along one line. But, surprisingly, bisplrep can calculate 2D spline interpolation in this case. I built 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) 

The result is the following:

where bisplrep seems to fill in the blanks 0, as best shown when I expand the chart below:

Regarding the expected addition of 0, my real question is: why does bisplrep work in case 2, but not in case 1?

Or, in other words: do we want it to return an error when two-dimensional interpolation is applied only with input in one direction (case 1 and 2 fail) or not? (Case 1 and 2 should return something, even if it is unpredictable).

Answer 1

I was originally going to show you what the difference is for 2d interpolation if your input data is oriented along the coordinate axes, and not in some general direction, but it turns out that the result will be even more messy than I expected. I tried using a random dataset over an interpolated rectangular grid and compared this to the case when the same x and y coordinates were rotated 45 degrees for interpolation. The result was terrible.

Then I tried to make a comparison with a smoother dataset: it turns out scipy.interpolate.interp2d has a lot of problems. Therefore, my bottom line will be "use scipy.interpolate.griddata ".

For instructive purposes, here is my (rather dirty) 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 

So here is the results gallery. Using data from random input z and interp2d , Cartesian (left) and rotational interpolation (right):

Pay attention to the terrible scale on the right side, noting that the input points are between 0 and 1 . Even his mother does not recognize the data set. Note that while evaluating a rotating dataset there are warnings at runtime, so we warn that all this shit.

Now do the same with griddata :

It should be noted that these numbers are much closer to each other, and they seem to make more sense than the output of interp2d . For example, pay attention to overshoot on the scale of the very first digit.

These artifacts always occur between input data points. Since it is still interpolation, the input points must be reproduced using the interpolation function, but it is rather strange that the linear interpolation function is overshoot between the data points. It is clear that griddata does not suffer from this problem.

Consider an even clearer case: another set of z values ​​that are smooth and deterministic. Surfaces with interp2d :

HELP! Call the interpolation police! Already in the Cartesian input case, inexplicable (well, at least by me) false signs in it, and a rotating input case poses a threat s g ̻͇͈͚̟̻͛ͫ͛̅͋͒ o ͈͓̱̥̙̫͚̾͂.

So do the same with griddata :

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

---

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. 

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

 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) 

And, of course, all this connects you with the following question: it is of great importance if your geometrically 1d-data set is located along one of the Cartesian axes or if it is generally such that the coordinate values ​​take different values. It is probably pointless (or at least very poorly defined) to try 2d interpolation from a geometrically 1d data set, but at least the algorithm should not be interrupted if your data is along the general direction of the x,y .