[SciPy-User] 2D interpolation on random points of a function defined on a regular grid

Antonio tritemio@gmail....
Wed Apr 24 14:49:49 CDT 2013


On Wed, Apr 24, 2013 at 12:20 AM, Pauli Virtanen <pav@iki.fi> wrote:
>
> Antonio <tritemio <at> gmail.com> writes:
> [clip]
> > import scipy.interpolate as SI
> > interp_fun = SI.interp2d(xd,yd,data)
> > z = interp_fun(x,y)
>
> Use
>
> interp_fun = SI.RectBivariateSpline(xd, yd, data.T)
> z = interp_fun.ev(x, y)
>
> Note that this expects the data be in (nx, ny) shape,
> hence the transpose.


Here we are! That's much better than map_coordinates. It's also
faster: a linear interpolation on 1e6 points run on 132ms in my laptop
(map_coordinates runs the same interpolation in 212ms). Thank you
Pauli!

In the meanwhile I also tried to use
scipy.interpolate.LinearNDInterpolator. The documentation is scarce
but I managed to make it work. The only caveat is that it requires a
meshgrid of the coordinates in which data is computed (I don't know if
is possible to avoid this). However the interpolation on 1e6 points
runs in 414ms in this case.

Here I attach the complete example that perform the same interpolation
with RectBivariateSpline, map_coordinates and LinearNDInterpolator.

Cheers,
Antonio

------
# Generate a simple 2D function
nx, ny = 5,4 # size of x and y axis
xd = arange(nx) + 3. # x axis
yd = arange(ny) - 4. # y axis
data = arange(nx*ny).reshape(ny,nx) # assume this is the result of
                                    #  a function evaluation on the Cartesian
                                    #  grid [xd , yd]
data = data.astype(float64)

# Generate some points where to compute the interpolation
N = 1e6
x = rand(N)*(xd.max()-xd.min()) + xd.min()
y = rand(N)*(yd.max()-yd.min()) + yd.min()

## Let try different interpolation methods

# interp2d doesn't work because (x,y) are unstructured
#import scipy.interpolate as SI
#interp_fun = SI.interp2d(xd,yd,data)
#z = interp_fun(x,ty) # Raises ValueError

## Method 1: map_coordinates (212ms on 1e6 points)
from scipy.ndimage import map_coordinates

def interp2d_map_c(xd,yd,data,xq,yq, **kwargs):
    nx, ny = xd.size, yd.size
    x_step, y_step = (xd[1]-xd[0]), (yd[1]-yd[0])
    assert (ny, nx) == data.shape
    assert (xd[-1] > xd[0]) and (yd[-1] > yd[0])

    # Translate to a unit-cell coordinate system (so that indexes are
    # coordinates)
    # This mapping requires a regular (uniform) grid for (xd,yd)
    # In principle could be extended to a non-uniform rectilinear grid
    xp = (xq-xd[0])*(nx-1)/(xd[-1]-xd[0])
    yp = (yq-yd[0])*(ny-1)/(yd[-1]-yd[0])
    coord = vstack([yp,xp])
    zq = map_coordinates(data, coord, **kwargs)
    return zq

z = interp2d_map_c(xd,yd,data, x,y, order=1)

## Method 2: LinearNDInterpolator (414ms on 1e6 points)
Xd,Yd = meshgrid(xd,yd)
coord = hstack([Xd.ravel().reshape(Xd.size,1), Yd.ravel().reshape(Yd.size,1)])
interp_fun = SI.LinearNDInterpolator(coord,data.ravel())
z = interp_fun(x,y).ravel()

## Method 3: RectBivariateSpline (131ms on 1e6 points) BEST
interp_fun = SI.RectBivariateSpline(xd, yd, data.T, kx=1, ky=1)
z = interp_fun.ev(x,y)


More information about the SciPy-User mailing list