Return 2-D array from scipy Regular Grid Interpolator

Question

I'm using Scipy to interpolate data in 6 dimensions, and would like a way to return a two dimensional array from my interpolator object instead of a one-D array. At the moment, I can only accomplish this by using a for loop and calling the interpolator object many times- I'm hoping there's a better method out there.

For example, in 3D:

#Random data
points=(np.arange(10), np.arange(5), np.arange(5))
values=np.random.rand(10, 5, 5)

interp_object=scipy.interpolate.RegularGridInterpolator(points, values)

I'd like to be able to do:

#Pick some points along each axis
x_points=np.arange(10)
y_points=[1.2, 3.1, 1.4, 4.8, 0.1]
z_points=3.5


xi=(x_points, y_points, z_points)

out_array=interp_object(xi)

But this leads to

ValueError: shape mismatch: objects cannot be broadcast to a single shape

and I have to end up doing:

out_array=np.empty((10, 5))
for i, y in enumerate(y_points):
    out_array[:, i]=interp_object((x_points, y, z_points))

This loop is the main bottleneck in my code and I'd like to avoid it if possible! Is what I'd like to do possible, using RegularGridInterpolator or another interpolation method?

Answer 1

This part of the code:

#Pick some points along each axis
x_points=np.arange(10)
y_points=[1.2, 3.1, 1.4, 4.8, 0.1]
z_points=3.5

xi=(x_points, y_points, z_points)

out_array=interp_object(xi)

will not work, as you have to give the inputs as an array of all the points you want. So we need to generate a matrix with all the combinations. For this we can use meshgrid. We also need to include some manipulation with the array dimensions to get everything working, so it might look a bit messy. Below is an example:

#Pick some points along each axis
x_points = np.arange(9) + 0.5  # Points must be inside
y_points = np.array([1.2, 3.1, 1.4, 3.9, 0.1])
z_points = np.array([3.49])

points = np.meshgrid(x_points, y_points, z_points)
flat = np.array([m.flatten() for m in points])
out_array = interp_object(flat.T)
result = out_array.reshape(*points[0].shape)

Also notice I've changed the points you gave a bit. The RegularGridInterpolator can only be used for interpolation, and will not work if used outside the ranges given when creating the interpolation function, which is called interp_object in this case.