[Python] Is there a clean pythonic way of querying multiple points on a bivariate spline in...

I need to efficiently evaluate a bivariate spline on a B-spline basis. I have already calculated the knot positions and spline coefficients (independently of scipy classes/methods such as RectBivariateSpline which calculate these for you). Hence I only want to evaluate the spline function itself using the knots and coefficients I provide. This spline function needs to be evaluated on a cloud of points rather than a grid formed of a Cartesian product. I came across the following stackoverflow query which explains how to do it: Querying multiple points on a bivariate spline in the B-spline basis . For convenience, I reproduce the code which does what I want:

from scipy.interpolate import dfitpack

x = x.ravel(order="C")
y = y.ravel(order="C")
# NOTE: specify the order parameter to avoid leaving this to chance

# NOTE: access low-level fortran routine to carry out calculations (this
# runs much faster than any python equivalent, but the documentation
# is only available in the fortran routine or online)
psf_grid, ier = dfitpack.bispeu(
yknots, xknots, psf_spline.ravel(order="C"), degree, degree, y, x
)
if not ier == 0:
raise ValueError("Error in spline2grid. Error code from bispeu: %s" % ier)


However, accessing the low level Fortran routine bispeu directly no longer works in later versions of scipy. I get DeprecationWarning in scipy-1.14.1, and the program fails in later versions. Therefore I have resorted to doing this instead:

from scipy.interpolate import BivariateSpline

# set up BivariateSpline object (to access bispeu indirectly)
spline_object = BivariateSpline() # is this legal ?
spline_object.degrees = (degree, degree)
spline_object.tck = (yknots, xknots, psf_spline.ravel(order="C"))
x = x.ravel(order="C")
y = y.ravel(order="C")
psf_grid = spline_object(y, x, grid=False)


The above solution works, but I'm wondering if it's legal or if there is a better way of doing this. For instance, the documentation to BivariateSpline state "This class is meant to be subclassed, not instantiated directly." However, if I try using a subclass such as RectBivariateSpline, it typically wants to recalculate the knots and coefficients (which is a waste of time and doesn't correspond to what I want), and may also try to force me to calculate the spline on a Cartesian grid rather than a cloud of points.

Continue reading...