Python/Scipy Interpolation (map_coordinates)
I'm trying to do some interpolation with scipy. I've gone through many examples, but I'm not finding exactly what I want.
Let's say I have some data where the row and column variable can vary from 0 to 1. The delta changes between each row and column is not always the same (see below).
| 0.00 0.25 0.80 1.00
------|----------------------------
0.00 | 1.40 6.50 1.50 1.80
0.60 | 8.90 7.30 1.10 1.09
1.00 | 4.50 9.20 1.80 1.20
Now I want to be able to take a set of x,y points and determine the interpolated values. I know I can do this with map_coordinates. I'm wondering if there is any easy/clever way to make an x,y value to the appropriate index in the data array.
For example, if I input x,y = 0.60, 0.25, then I should get back the correct index to be interpolated. In this case, that would be 1.0, 1.0 since 0.60, 0.25 would map exactly to the second row and second column. x=0.3 would map to 0.5 since it is halfway between 0.00 and 0.60.
I know how to get the result I want, but I feel certain that there is a very quick/clear one or two-liner (or function that already exists) that can do this to make my code more clear. Basically it needs to piecewise interpolate between some array.
Here is an example (based heavily on the code from Scipy interpolation on a numpy array) - I put TODO where this new function would go:
from scipy.ndimage import map_coordinates
from numpy import arange
import numpy as np
# 0.000, 0.175, 0.817, 1.000
z = array([ [ 3.6, 6.5, 9.1, 11.5], # 0.0000
[ 3.9, -7.3, 10.0, 13.1], # 0.2620
[ 1.9, 8.3, -15.0, -12.1], # 0.6121
[-4.5, 9.2, 12.2, 14.8] ]) # 1.0000
ny, nx = z.shape
xmin, xmax = 0., 1.
ymin, ymax = 0., 1.
xrange = array([0.000, 0.175, 0.817, 1.000 ])
yrange = array([0.0000, 0.2620, 0.6121, 1.0000])
# Points we want to interpolate at
x1, y1 = 0.20, 0.45
x2, y2 = 0.30, 0.85
x3, y3 = 0.95, 1.00
# To make our lives easier down the road, let's
# turn these into arrays of x & y coords
xi = np.array([x1, x2, x3], dtype=np.float)
yi = np.array([y1, y2, y3], dtype=np.float)
# Now, we'll set points outside the boundaries to lie along an edge
xi[xi > xmax] = xmax
xi[xi < xmin] = xmin
yi[yi > ymax] = ymax
yi[yi < ymin] = ymin
# We need to convert these to (float) indicies
# (xi should range from 0 to (nx - 1), etc)
xi = (nx - 1) * (xi - xmin) / (xmax - xmin)
yi = (ny - 1) * (yi - ymin) / (ymax - ymin)
# TODO: Instead, xi and yi need to be mapped as described. This can only work with
# even spacing...something like:
#xi = SomeInterpFunction(xi, xrange)
#yi = SomeInterpFunction(yi, yrange)
# Now we actually interpolate
# map_coordinates does cubic interpolation by default,
# use "order=1" to preform bilinear interpolation instead...
print xi
print yi
z1, z2, z3 = map_coordinates(z, [yi, xi], order=1)
# Display the results
for X, Y, Z in zip((x1, x2, x3), (y1, y2, y3), (z1, z2, z3)):
print X, ',', Y, '-->', Z
Answer
I think you want a bivariate spline on a rectangular structured mesh:
import numpy
from scipy import interpolate
x = numpy.array([0.0, 0.60, 1.0])
y = numpy.array([0.0, 0.25, 0.80, 1.0])
z = numpy.array([
[ 1.4 , 6.5 , 1.5 , 1.8 ],
[ 8.9 , 7.3 , 1.1 , 1.09],
[ 4.5 , 9.2 , 1.8 , 1.2 ]])
# you have to set kx and ky small for this small example dataset
# 3 is more usual and is the default
# s=0 will ensure this interpolates. s>0 will smooth the data
# you can also specify a bounding box outside the data limits
# if you want to extrapolate
sp = interpolate.RectBivariateSpline(x, y, z, kx=2, ky=2, s=0)
sp([0.60], [0.25]) # array([[ 7.3]])
sp([0.25], [0.60]) # array([[ 2.66427408]])