[SciPy-user] Bilinear Interpolation

[Ryan May]

Stefan van der Walt wrote:
> Hi Ryan
> 
> On Wed, Sep 27, 2006 at 08:09:18AM -0500, Ryan May wrote:
>>> That's what I kinda thought, but using 'linear' for interp2d didn't give 
>>> me the same answer as when I performed the calculation manually.  I'll 
>>> have to see what I can find on this...
>> Ok, if I select the 4 points surrounding the location of interest, 
>> interp2d gives me the value I get with manually calculating a bilinear 
>> interpolation.  However if I use the whole field (or even 9 points 
>> instead of 4), I get a different answer.  I _know_ I wouldn't expect 
>> this for bilinear interpolation, and it would seem to imply that even 
>> _linear_ B-splines uses the information from additional points.  Am I 
>> missing something here, or are the methods just not truly
>> equivalent?
> 
> Please send some code so we can see what you are doing.  Without that,
> we are just discussing things in the air.
> 

Here's some example code that demonstrates what I mean, but I think 
Scott Sinclair hit the nail on the head in another reply.  The methods 
are just different except in the case of 4 points, since the splines are 
performing a global fit and the bilinear interpolation is only local.

from scipy.interpolate import interpolate
from numpy.random import randn
from numpy import *

data = randn(16).reshape(4,4)
x = y = arange(4)
X,Y = meshgrid(x,y)

x_loc = 1.4
y_loc = 2.7

#Do the interpolation with all points, a partial set, and only 4 points,
#making sure that the deisred location is within the domain
full_bilin = interpolate.interp2d(X,Y,data,kind='linear')
partial_bilin = 
interpolate.interp2d(X[1:,1:],Y[1:,1:],data[1:,1:],kind='linear')
single_bilin = 
interpolate.interp2d(X[2:4,1:3],Y[2:4,1:3],data[2:4,1:3],kind='linear')

#If interp2d was truly doing bilinear interpolation, these should be equal
print full_bilin(x_loc,y_loc)
print partial_bilin(x_loc,y_loc)
print single_bilin(x_loc,y_loc)

#This manual calculation of the interpolation is only equal to the last one
xw = 0.4
yw = 0.7
xws = array([1-xw, xw])
yws = array([1-yw, yw])
print dot(yws,dot(data[2:4,1:3],xws))


Ryan May
Graduate Research Assistant
School of Meteorology
University of Oklahoma
rmay at rossby.ou.edu