在Python中包装（圆形）二维插值

Question

我在以pi弧度包裹的域（即0 = pi）上有角度数据。数据为2D，其中一维代表角度。我需要以包装方式将此数据插值到另一个网格上。

在一个维度上，该np.interp函数需要一个时间段kwarg（适用于NumPy1.10及更高版本）：http : //docs.scipy.org/doc/numpy/reference/generated/numpy.interp.html

这正是我需要的，但是我需要两个维度。我目前只是逐步浏览数组中的列并使用np.interp，但这当然很慢。

有什么可以达到相同结果但更快的结果吗？

Answer 1

np.interp工作原理的解释

使用来源，卢克！

numpy 文档np.interp使得源代码特别容易找到，因为它有链接以及文档。让我们逐行讨论一下。

首先回忆一下参数：

"""
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.
"""

让我们以一个简单的三角波为例来进行说明：

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

period != None现在，在所有类型检查发生之后，我从源代码中的分支开始：

# normalizing periodic boundaries
x = x % period
xp = xp % period

这只是确保提供的所有值x和xp都在0和之间period。因此，由于周期为pi，但我们指定x和xp位于-pi/2和 之间，因此将通过添加范围 中的所有值pi/2来进行调整，以便它们有效地出现在 后。所以我们现在读的是.pi[-pi/2, 0)pi/2xp[pi/2, 3*pi/4, 0, pi/4]

asort_xp = np.argsort(xp)
xp = xp[asort_xp]
fp = fp[asort_xp]

这只是按xp递增顺序排序。在执行上一步中的模运算后尤其需要这样做。那么，现在xp是[0, pi/4, pi/2, 3*pi/4]。fp也相应地进行了洗牌[0, 1, 0, -1]。

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

np.interp进行线性插值。当尝试在两点之间插值a并b出现在 中时，它仅使用和xp的值（即相应索引处的值）。所以最后一步所做的就是取出点并将其放在数组前面，并将点放在数组后面，但分别减去和添加一个句点。所以你现在有了一个看起来像的新的。同样，和已串联起来，现在 也是如此。f(a)f(b)fpnp.interpxp[-1]xp[0]xp[-pi/4, 0, pi/4, pi/2, 3*pi/4, pi]fp[0]fp[-1]fp[-1, 0, 1, 0, -1, 0]

请注意，在模运算之后，也x已进入[0, pi]范围，x现在 也是如此[7*pi/8, 3*pi/8]。这让您很容易看到您会回来[-0.5, 0.5]。

现在，来看您的 2D 案例：

假设您有一个网格和一些值。让我们立即将所有值取为介于[0, pi]两者之间，这样我们就不需要担心模数和洗牌。

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, :] 
# >>> fp
# array([[ 0,  1,  2,  3],
#        [ 1,  2,  3,  4],
#        [ 0,  1,  2,  3],
#        [-1,  0,  1,  2]])

xp[[-1]]我们现在知道您需要做的就是在数组前面和xp[[0]]末尾添加，并根据句点进行调整。请注意我如何使用单例列表[-1]和[0]. 这是确保保留尺寸的技巧。

xp = np.concatenate((xp[[-1]]-period, xp, xp[[0]]+period))
fp = np.concatenate((fp[[-1], :], fp, fp[[0], :]))

最后，您可以自由使用scipy.interpolate.interpn来实现您的结果。x = pi/8让我们得到所有的值y：

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)
# >>> result
# array([ 0.5,  1.5,  2.5,  3.5])

interp_points必须指定为 Nx2 点矩阵，其中第一个维度是您想要在第二个维度进行插值的每个点，给出该点的 x 和 y 坐标。请参阅此答案以获取详细说明。

如果你想获得超出 range 的值[0, period]，你需要自己对其取模：

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)
# >>> result
# array([-0.5,  0.5,  1.5,  2.5])

同样，如果您想生成interp_points一堆 x 和 y 值，请查看此答案。