python - Python 中的环绕(圆形)2D 插值

Question

我有一个域的角度数据，该数据以 pi 弧度(即 0 = pi)包裹。数据是 2D，其中一维表示角度。我需要以包装方式将这些数据插入到另一个网格中。

在一个维度上，np.interp 函数采用句点 kwarg(对于 NumPy 1.10 及更高版本): http://docs.scipy.org/doc/numpy/reference/generated/numpy.interp.html

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

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

Answer 1

如何解释 np.interp作品

使用 source ，卢克!

numpy doc for 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/2 之后.所以我们的xp现在显示 [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) (即 fp 在相应索引处的值)。那又怎样np.interp这最后一步做的就是取点xp[-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] .这是 trick确保dimensions are preserved .

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 坐标。参见 this answer以获得详细的解释。

如果你想得到[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 值，请查看 this answer .