[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
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()
## 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.
Cheers,
Antonio
-------------- 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