How to pass arrays into Scipy Interpolate RectBivariateSpline?

Question

I am creating a Scipy Interpolate RectBivariateSpline as follows:

import numpy as np
from scipy.interpolate import RectBivariateSpline

x = np.array([1,2,3,4])
y = np.array([1,2,3,4,5])
vals = np.array([
    [4,1,4,4,2],
    [4,2,3,2,6],
    [3,7,4,3,5],
    [2,4,5,3,4]
])

rect_B_spline = RectBivariateSpline(x, y, vals)

I then try to pass-in an array of x and y points:

a = np.array([3.2, 3.8, 2.2])
b = np.array([2.4, 4.3, 3.3])

print(rect_B_spline(a, b))

To which I get an error as follows:

Traceback (most recent call last):
  File "path/file", line 18, in <module>
    print(rect_B_spline(a, b))
  File "/path/scipy/interpolate/fitpack2.py", line 728, in __call__
    raise ValueError("Error code returned by bispev: %s" % ier)
ValueError: Error code returned by bispev: 10

This error is remedied when I pass in the grid=False parameter to the method.

My impression from the documentation was that if the input grid coordinates form a regular grid, then the grid parameter should be True. Is there something I am missing?

Solution 1

Yes, the documentation is perhaps a bit weak here. The default call that you are using expects that x and y define grid points. Like the original call to create the spline fit, these need to be in strictly ascending order. It will then return the full grid of spline interpolations.

If you instead use RectBivariateSpline.ev(xi,yi), in your case call rect_B_spline.ev(a,b), you get the spline evaluated at each of (xi[0],yi[0]), ..., (xi[j],yi[j]) data pairs.

Not quite sure which you wanted here - if you want a full x by y grid, make the points in each of x,y to be strictly ascending. If you want the results at a series of points, use the .ev(x,y) method.

Solution 2

This works:

RectBivariateSpline(x,y,vals)([2.2,3.2,3.8],[2.4,3.3,4.3],grid=True)

array([[ 3.197056,  2.75356 ,  2.38796 ],
       [ 6.408896,  3.56696 ,  3.46736 ],
       [ 5.768704,  4.33264 ,  3.10824 ]])

RectBivariateSpline(x[:,None],y,vals)(a,b,grid=False)
# array([ 6.408896,  3.10824 ,  2.75356 ])

I sorted the values, as specified in the docs:

If grid is True: evaluate spline at the grid points defined by the coordinate arrays x, y. The arrays must be sorted to increasing order.