griddata scipy interpolation not working (giving nan)

Question

I was trying out the 2d example given in the scipy.interpolation.griddata help file. It works for interpolation with 'nearest'. But it gives a matrix filled with nan while using any other interpolation like 'linear' or 'cubic'. If I give the argument fill_value=5 , it gives the matrix filled with 5.

Is this due to some installation problem?

I was trying the exact same thing they have given in the help document. But somehow it is giving the result as if the points i asked to interpolating is lying outside the input points. (which is not!! I followed the example)

I shall post the example to reproduce the error (taken form doc)

def func(x, y):
    return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2

grid_x, grid_y = np.mgrid[0:1:10j, 0:1:10j]
points = np.random.rand(100, 2)
values = func(points[:,0], points[:,1])

from scipy.interpolate import griddata

grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest')
grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear')
grid_z2 = griddata(points, values, (grid_x, grid_y), method='cubic')

I am getting grid_z1 and grid_z2 to be a matrix filled with nan.

UPDATE : I installed all the packages in another Ubuntu 11.10 machine. And the same script gave perfectly correct answer. Previously I was trying on Porteus distro (live slackware family). Hence I think i can safely conclude that this was some problem in my installation. Anybody have any idea what could have gone wrong? Does any library conflict result in this kind of behavior? Since my main machine is Portues, i have no other option than to repair the scipy in it.

Answer 1

You say "filled with nan", but it's not really filled. Using your code but adding

np.random.seed(7)

at the start so that we're working with the same dataset, I find

>>> np.isnan(grid_z1).sum()
744
>>> np.isnan(grid_z2).sum()
744

And these NaNs occur all on a band on the outside:

>>> np.isnan(grid_z1[5:-5,5:-5]).sum()
0

which makes it likely what the problem is. The points which are giving NaN are outside the specified points, so it doesn't know what to do with them. For the special case of "nearest" interpolation, you can still find something that's near, so you don't get any NaNs out.

So when you say the points to be interpolated at aren't lying outside the input points, I beg to differ:

# brute force, because I'm too lazy
from collections import Counter
d = Counter()
for x, y, val in zip(grid_x.flat, grid_y.flat, grid_z1.flat):
    pg = (points >= [x, y])
    boxed = len(set(tuple(p) for p in pg)) == 4
    d[np.isnan(val), boxed] += 1

produces

>>> d
Counter({(False, True): 19189, (True, False): 744, (False, False): 67})

And there are no (True, True) cases. IOW, every NaN lacks a bounding box in the points. There are some (False, False) cases, where the value doesn't have a bounding box but doesn't wind up a NaN, which is mildly surprising, but if they've assumed that everything is contained it would probably depend upon boring implementation details what happens if they're not. Short version: I think everything here is probably working correctly, in the sense of as expected.