### Question :

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
```

Returns: `ValueError: Invalid inputs`

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**.

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

*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)
```

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

The result is the following:

where `bisplrep`

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

Regarding of whether adding 0 is expected, 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 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).

##
Answer #1:

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.

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
```

So here’s a gallery of the results. Using random input `z`

data, and `interp2d`

, Cartesian (left) vs rotated interpolation (right):

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.

Now let’s do the same with `griddata`

:

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.

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.

Consider an even more clear case: the other set of `z`

values, which are smooth and deterministic. The surfaces with `interp2d`

:

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??????????????u??????????m????????m???????o???n??????i??????n???????????g??????? ?z????????a????????l???????????g????????????o?????????.

So let’s 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 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)
```

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.