Wrapped (circular) 2D interpolation in python
<p>I have angular data on a domain that is wrapped at pi radians (i.e. 0 = pi). The data are 2D, where one dimension represents the angle. I need to interpolate this data onto another grid in a wrapped way.</p>

<p>In one dimension, the np.interp function takes a period kwarg (for numpy 1.10 and later): 
<a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.interp.html" rel="nofollow">http://docs.scipy.org/doc/numpy/reference/generated/numpy.interp.html</a></p>

<p>This is exactly what I need, but I need it in two dimensions. I'm currently just stepping through columns in my array and using np.interp, but this is of course slow. Anything out there that could achieve this same outcome but faster?</p>

<h2>An explanation of how <code>np.interp</code> works</h2>

<p>Use the <a href="https://github.com/numpy/numpy/blob/v1.11.0/numpy/lib/function_base.py#L1570-L1692" rel="nofollow">source</a>, Luke!</p>

<p>The <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.interp.html" rel="nofollow">numpy doc for <code>np.interp</code></a> makes the source particularly easy to find, since it has the link right there, along with the documentation. Let's go through this, line by line.</p>

<p>First, recall the parameters:</p>

<pre><code>"""
x : array_like
    The x-coordinates of the interpolated values.
xp : 1-D sequence of floats
    The x-coordinates of the data points, must be increasing if argument
    `period` is not specified. Otherwise, `xp` is internally sorted after
    normalizing the periodic boundaries with ``xp = xp % period``.
fp : 1-D sequence of floats
    The y-coordinates of the data points, same length as `xp`.
period : None or float, optional
    A period for the x-coordinates. This parameter allows the proper
    interpolation of angular x-coordinates. Parameters `left` and `right`
    are ignored if `period` is specified.
"""
</code></pre>

<p>Let's take a simple example of a triangular wave while going through this:</p>

<pre><code>xp = np.array([-np.pi/2, -np.pi/4, 0, np.pi/4])
fp = np.array([0, -1, 0, 1])
x = np.array([-np.pi/8, -5*np.pi/8])  # Peskiest points possible }:)
period = np.pi
</code></pre>

<p>Now, I start off with the <code>period != None</code> branch in the source code, after all the type-checking happens:</p>

<pre><code># normalizing periodic boundaries
x = x % period
xp = xp % period
</code></pre>

<p>This just ensures that all values of <code>x</code> and <code>xp</code> supplied are between <code>0</code> and <code>period</code>. So, since the period is <code>pi</code>, but we specified <code>x</code> and <code>xp</code> to be between <code>-pi/2</code> and <code>pi/2</code>, this will adjust for that by adding <code>pi</code> to all values in the range <code>[-pi/2, 0)</code>, so that they effectively appear after <code>pi/2</code>. So our <code>xp</code> now reads <code>[pi/2, 3*pi/4, 0, pi/4]</code>.</p>

<pre><code>asort_xp = np.argsort(xp)
xp = xp[asort_xp]
fp = fp[asort_xp]
</code></pre>

<p>This is just ordering <code>xp</code> in increasing order. This is especially required after performing that modulo operation in the previous step. So, now <code>xp</code> is <code>[0, pi/4, pi/2, 3*pi/4]</code>. <code>fp</code> has also been shuffled accordingly, <code>[0, 1, 0, -1]</code>.</p>

<pre><code>xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period))
fp = np.concatenate((fp[-1:], fp, fp[0:1]))
return compiled_interp(x, xp, fp, left, right)   # Paraphrasing a little
</code></pre>

<p><code>np.interp</code> does linear interpolation. When trying to interpolate between two points <code>a</code> and <code>b</code> present in <code>xp</code>, it only uses the values of <code>f(a)</code> and <code>f(b)</code> (i.e., the values of <code>fp</code> at the corresponding indices). So what <code>np.interp</code> is doing in this last step is to take the point <code>xp[-1]</code> and put it in front of the array, and take the point <code>xp[0]</code> and put it after the array, but after subtracting and adding one period respectively. So you now have a new <code>xp</code> that looks like <code>[-pi/4, 0, pi/4, pi/2, 3*pi/4, pi]</code>. Likewise, <code>fp[0]</code> and <code>fp[-1]</code> have been concatenated around, so <code>fp</code> is now <code>[-1, 0, 1, 0, -1, 0]</code>.</p>

<p>Note that after the modulo operations, <code>x</code> had been brought into the <code>[0, pi]</code> range too, so <code>x</code> is now <code>[7*pi/8, 3*pi/8]</code>. Which lets you easily see that you'll get back <code>[-0.5, 0.5]</code>.</p>

<hr>

<h2>Now, coming to your 2D case:</h2>

<p>Say you have a grid and some values. Let's just take all values to be between <code>[0, pi]</code> off the bat so we don't need to worry about modulos and shufflings.</p>

<pre><code>xp = np.array([0, np.pi/4, np.pi/2, 3*np.pi/4])
yp = np.array([0, 1, 2, 3])
period = np.pi

# Put x on the 1st dim and y on the 2nd dim; f is linear in y
fp = np.array([0, 1, 0, -1])[:, np.newaxis] + yp[np.newaxis, :] 
# &gt;&gt;&gt; fp
# array([[ 0,  1,  2,  3],
#        [ 1,  2,  3,  4],
#        [ 0,  1,  2,  3],
#        [-1,  0,  1,  2]])
</code></pre>

<p>We now know that all you need to do is to add <code>xp[[-1]]</code> in front of the array and <code>xp[[0]]</code> at the end, adjusting for the period. Note how I've indexed using the singleton lists <code>[-1]</code> and <code>[0]</code>. This is a <a href="http://stackoverflow.com/a/18183182/525169">trick</a> to ensure that <a href="http://stackoverflow.com/q/3551242/525169">dimensions are preserved</a>.</p>

<pre><code>xp = np.concatenate((xp[[-1]]-period, xp, xp[[0]]+period))
fp = np.concatenate((fp[[-1], :], fp, fp[[0], :]))
</code></pre>

<p>Finally, you are free to use <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html#scipy.interpolate.interpn" rel="nofollow"><code>scipy.interpolate.interpn</code></a> to achieve your result. Let's get the value at <code>x = pi/8</code> for all <code>y</code>:</p>

<pre><code>from scipy.interpolate import interpn
interp_points = np.hstack(( (np.pi/8 * np.ones(4))[:, np.newaxis], yp[:, np.newaxis] ))
result = interpn((xp, yp), fp, interp_points)
# &gt;&gt;&gt; result
# array([ 0.5,  1.5,  2.5,  3.5])
</code></pre>

<p><code>interp_points</code> has to be specified as an Nx2 matrix of points, where the first dimension is for each point you want interpolation at the second dimension gives the x- and y-coordinate of that point. See <a href="http://stackoverflow.com/a/39357219/525169">this answer</a> for a detailed explanation.</p>

<p>If you want to get the value outside of the range <code>[0, period]</code>, you'll need to modulo it yourself:</p>

<pre><code>x = 21 * np.pi / 8
x_equiv = x % period   # Now within [0, period]
interp_points = np.hstack(( (x_equiv * np.ones(4))[:, np.newaxis], yp[:, np.newaxis] ))
result = interpn((xp, yp), fp, interp_points)
# &gt;&gt;&gt; result
# array([-0.5,  0.5,  1.5,  2.5])
</code></pre>

<p>Again, if you want to generate <code>interp_points</code> for a bunch of x- and y- values, look at <a href="http://stackoverflow.com/a/39357219/525169">this answer</a>.</p>

python numpy interpolation