I'll try to explain this to you in 2D so that you get a better idea of what's happening. First, let's create a linear array to test with.
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
# Set up grid and array of values
x1 = np.arange(10)
x2 = np.arange(10)
arr = x1 + x2[:, np.newaxis]
# Set up grid for plotting
X, Y = np.meshgrid(x1, x2)
# Plot the values as a surface plot to depict
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, arr, rstride=1, cstride=1, cmap=cm.jet,
linewidth=0, alpha=0.8)
fig.colorbar(surf, shrink=0.5, aspect=5)
This gives us:
Then, let's say you want to interpolate along a line, i.e., one point along the first dimension, but all points along the second dimension. These points are not in the original arrays (x1, x2)
obviously. Suppose we want to interpolate to a point x1 = 3.5
, which is in between two points on the x1-axis.
from scipy.interpolate import interpn
interp_x = 3.5 # Only one value on the x1-axis
interp_y = np.arange(10) # A range of values on the x2-axis
# Note the following two lines that are used to set up the
# interpolation points as a 10x2 array!
interp_mesh = np.array(np.meshgrid(interp_x, interp_y))
interp_points = np.rollaxis(interp_mesh, 0, 3).reshape((10, 2))
# Perform the interpolation
interp_arr = interpn((x1, x2), arr, interp_points)
# Plot the result
ax.scatter(interp_x * np.ones(interp_y.shape), interp_y, interp_arr, s=20,
c='k', depthshade=False)
plt.xlabel('x1')
plt.ylabel('x2')
plt.show()
This gives you the result as desired: note that the black points correctly lie on the plane, at an x1 value of 3.5
.
Note that most of the "magic", and the answer to your question, lies in these two lines:
interp_mesh = np.array(np.meshgrid(interp_x, interp_y))
interp_points = np.rollaxis(interp_mesh, 0, 3).reshape((10, 2))
I have explained the working of this elsewhere. In short, what it does is to create an array of size 10x2, containing the coordinates of the 10 points you want to interpolate arr
at. (The only difference between that post and this one is that I've written that explanation for np.mgrid
, which is a shortcut to writing np.meshgrid
for a bunch of arange
s.)
For your 4x4x4x4 case, you will probably need something like this:
interp_mesh = np.meshgrid([0.1], [9], np.linspace(0, 30, 3),
np.linspace(0, 0.3, 4))
interp_points = np.rollaxis(interp_mesh, 0, 5)
interp_points = interp_points.reshape((interp_mesh.size // 4, 4))
result = interpn(points, arr, interp_points)
Hope that helps!