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

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

On Wed, Apr 24, 2013 at 2:22 AM, denis <denis-bz-py@t-online.de> wrote:

> Antonio <tritemio <at> gmail.com> writes:
> > I'm trying (with no success) to compute the interpolated values of a 2D
> > function that is known on a regular grid X,Y.  Although the function is
> > defined on a regular grid I want to compute the interpolation on
> arbitrary
> > positions.
> See the nice example of ndimage.map_coordinates under
> http://stackoverflow.com/questions/6238250/multivariate-spline-interpolation-in-python-scipy

Thanks, now I see how to use it. However, in order to use
map_coordinates(), we need to translate the coordinates at which data is
known to a grid starting at (0.0) and having step 1 in both dimensions.

Here I attach a reworked example using map_coordinates. The translation of
coordinates is a bit convoluted, I don't know if there is a cleaner way to
do it.

# 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
                                    #  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()

## Interpolation based on map_coordinates()
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])

    # If one coordinate size is 1 assume it to be constant
    if size(xq) == 1 and size(yq) > 1:
        xq = xq*ones(yq.size)
    elif size(yq) == 1 and size(xq) > 1:
        yq = yq*ones(xq.size)

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

For 1e6 points the last lines runs in 212ms on my laptop.

> (Some knowledgeable people like Fitpack, some map_coordinates;
> both need better doc.)

Totally agree on that! It's frustrating and time consuming to search a
solutions for such a simple problem.

In matlab this interpolation can be performed in a single line by interp2.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20130424/b5e8da78/attachment.html 

More information about the SciPy-User mailing list