[SciPy-Dev] bilinear interpolation

[Eric Firing]

On 06/12/2011 04:32 PM, Ole Nielsen wrote:
> Hi again Chuck and thanks for helping me.
>
>  >Can you be a bit more specific about the
>  >grid and how you want to resample it? That is, do you need to sample it a
>  >random points, or just resample it on a different grid, etc.
>
> Our task is not about resampling but (the simpler) task of interpolating
> from a regular grid (e.g. earthquake ground shaking) to a collection of
> arbitrary points (e.g. location of critical infrastructure) - see e.g
> the illustration at http://en.wikipedia.org/wiki/Bilinear_interpolation
>
> so we need
>
>    1. A simple and robust way of interpolating and from a regular grid
>       to points without overshoot and ringing. Smoothing is not important.
>    2. Ability to hand NaN sensibly. E.g. if the is one grid point which
>       is NaN only interpolation points near it should be NaN.
>
> It is interesting how implementations of the hard problems (e.g. high
> order splines on variable resolution) are abundant, whereas something as
> basic as this is harder to find :-)

There is an implementation of bilinear interpolation in the matplotlib 
basemap toolkit, in the module __init__.py.  I don't think it depends on 
anything else in the module.  It optionally uses masked arrays to handle 
missing values; maybe it would work with nans also, but if not, then 
converting to a masked array and back to an ndarray with nans is fast 
and easy.

See
https://github.com/matplotlib/basemap
for the whole package, or
https://github.com/matplotlib/basemap/blob/master/lib/mpl_toolkits/basemap/__init__.py
for interp().

The docstring says the x and y output arrays have to be 2-D, but there 
is actually no such restriction; arbitrary x and y arrays will suffice 
(though they do have to be numpy arrays):


from mpl_toolkits.basemap import interp
xin = np.arange(5)
yin = np.arange(10)
zin = yin[:, np.newaxis] * xin[np.newaxis, :]
xout = [2.2, 3.3]
yout = [4.4, 5.5]
xout = np.array(xout)
yout = np.array(yout)
zout = interp(zin, xin, yin, xout, yout)
print zout

Eric

>
> I am quite keen to have a go at implementing a simple bilinear
> interpolator as a Python C-extension using numpy arrays if there is
> interest in it.
> Meanwhile, if you could tell me where the LinearNDInterpolator package
> can be found?
>
> Cheers and thanks
> Ole