Using scipy.interpolate.interpn to interpolate a N-Dimensional array
<p>Suppose I have data that depends on 4 variables: a, b, c and d. I want interpolate to return a 2D array which corresponds to a single value of a and b, and an array of values for c and d. However, The array size need not be the same. To be specific, my data comes from a transistor simulation. The current depends on 4 variables here. I want to plot a parametric variation. The number of points on the parameter is a lot less than the number of points for the horizontal axis. </p>

<pre><code>import numpy as np
from scipy.interpolate import interpn
arr = np.random.random((4,4,4,4))
x1 = np.array([0, 1, 2, 3])
x2 = np.array([0, 10, 20, 30])
x3 = np.array([0, 10, 20, 30])
x4 = np.array([0, .1, .2, .30])
points = (x1, x2, x3, x4)
</code></pre>

<p>The following works:</p>

<pre><code>xi = (0.1, 9, np.transpose(np.linspace(0, 30, 4)), np.linspace(0, 0.3, 4))
result = interpn(points, arr, xi)
</code></pre>

<p>and so does this:</p>

<pre><code>xi = (0.1, 9, 24, np.linspace(0, 0.3, 4))
result = interpn(points, arr, xi)
</code></pre>

<p>but not this:</p>

<pre><code>xi = (0.1, 9, np.transpose(np.linspace(0, 30, 3)), np.linspace(0, 0.3, 4))
result = interpn(points, arr, xi)
</code></pre>

<p>As you can see, in the last case, the size of the last two arrays in <code>xi</code> is different. Is this kind of functionality not supported by scipy or am I using <code>interpn</code> incorrectly? I need this create a plot where one of the <code>xi</code>'s is a parameter while the other is the horizontal axis.</p>

<p>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.</p>

<pre><code>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)
</code></pre>

<p>This gives us:
<a href="http://i.stack.imgur.com/kOAqE.png" rel="nofollow"><img src="http://i.stack.imgur.com/kOAqE.png" alt="surface plot of values"></a></p>

<p>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 <code>(x1, x2)</code> obviously. Suppose we want to interpolate to a point <code>x1 = 3.5</code>, which is in between two points on the x1-axis.</p>

<pre><code>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()
</code></pre>

<p>This gives you the result as desired: note that the black points correctly lie on the plane, at an x1 value of <code>3.5</code>.
<a href="http://i.stack.imgur.com/AKDzK.png" rel="nofollow"><img src="http://i.stack.imgur.com/AKDzK.png" alt="surface plot of interpolated points"></a></p>

<p>Note that most of the "magic", and the answer to your question, lies in these two lines:</p>

<pre><code>interp_mesh = np.array(np.meshgrid(interp_x, interp_y))
interp_points = np.rollaxis(interp_mesh, 0, 3).reshape((10, 2))
</code></pre>

<p>I have explained the working of this <a href="http://stackoverflow.com/a/27286794/525169">elsewhere</a>. In short, what it does is to create an array of size 10x2, containing the coordinates of the 10 points you want to interpolate <code>arr</code> at. (The only difference between that post and this one is that I've written that explanation for <code>np.mgrid</code>, which is a shortcut to writing <code>np.meshgrid</code> for a bunch of <code>arange</code>s.)</p>

<p>For your 4x4x4x4 case, you will probably need something like this:</p>

<pre><code>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)
</code></pre>

<p>Hope that helps!</p>

python python-2.7 numpy scipy