Evaluate the output from scipy 2D interpolation along a curve

user_3261485 0
Viewed 2447

I have data z sampled from a 2D function f at grid points x, y, as in z = f(x, y).

It is easy to interpolate f with scipy.interp2d via f = interp2d(x, y, z).

However, evaluating f(x, y) returns an entire 2D grid as if I had done

xx, yy = np.meshgrid(x, y)
f(xx, yy)

The behaviour I want is to simply return the values [f(x[i], y[i]) for i in range(len(x))], which I believe is the behaviour for pretty much any other method in numpy.

The reason I want this behaviour is that I'm looking for the path traced out along the surface of f over "time" by the pair (t, u(t)).

It is also surprising that np.diag(f(t, u(t))) differs from np.array([f(ti, u(ti)) for ti in t]), so It's not clear to me how to get at the path f(t, u(t)) from what is returned via interp2d.

EDIT: About diag, I just thought it seemed we should have np.diag(f(t, u(t))) == np.array([f(ti, u(ti)) for ti in t]), but that isn't the case.

Complete example:

def f(t, u):
    return (t**2) * np.exp((u**2) / (1 + u**2))

x = np.linspace(0, 1, 250)
xx, yy = np.meshgrid(x, x)

z = f(xx, yy)
f = scipy.interpolate.interp2d(x, y, z)

print(f(x, y))
print(np.array([f(xi, yi)[0] for xi, yi in zip(x, y)]))

I would like the output of both print statements to be the same.

scipy numpy python
3 Answers

The method interp2d returns an object whose call method expects the x, y vectors to be coordinates of a rectangular grid. And the reason you don't get the desired values from the diagonal of the returned array is that it sorts x, y first.

But there is a workaround, which I also used in Querying multiple points on a bivariate spline in the B-spline basis. After executing 

import scipy.interpolate as si
f = si.interp2d(x, y, z)

evaluate f not by calling it, but by passing its tck properties, followed by your x, y coordinates, to internal bispeu method. Like this:

print(si.dfitpack.bispeu(f.tck[0], f.tck[1], f.tck[2], f.tck[3], f.tck[4], x, y)[0])

The above returns the same as the slow loop

print(np.array([f(xi, yi)[0] for xi, yi in zip(x, y)]))

Explanation

The object f is secretly a B-spline of order 1. The spline parameters (knots, coefficients, order) are contained in its tck property and can be used directly by lower-order routines to the desired effect.

(Ideally, the call method of f would have a Boolean parameter grid which we'd set to False to let it know we don't want grid evaluation. Alas, it's not implemented.)

Admin1

It seems that interp2d() behaves they way it does because that is how the corresponding Fortran function was conceived. The only workaround to this (that I can think of) is to call f on pairs of coordinates:

[f(*p)[0] for p in zip(x, y)]
user_80335850

Based on user6655984's suggestion, I've posted the following wrapper function in another thread:

import scipy.interpolate as si
def interp2d_pairs(*args,**kwargs):
    """ Same interface as interp2d but the returned interpolant will evaluate its inputs as pairs of values.
    """
    # Internal function, that evaluates pairs of values, output has the same shape as input
    def interpolant(x,y,f):
        x,y = np.asarray(x), np.asarray(y)
        return (si.dfitpack.bispeu(f.tck[0], f.tck[1], f.tck[2], f.tck[3], f.tck[4], x.ravel(), y.ravel())[0]).reshape(x.shape)
    # Wrapping the scipy interp2 function to call out interpolant instead
    return lambda x,y: interpolant(x,y,si.interp2d(*args,**kwargs))

# Create the interpolant (same interface as interp2d)
f = interp2d_pairs(X,Y,Z,kind='cubic')
# Evaluate the interpolant on each pairs of x and y values
z=f(x,y)
user_6279110
Related