Python中的包裹(圆形)2D插值

2023年3月31日 226次阅读

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

在一个维度中,np.interp函数采用周期kwarg(对于NumPy 1.10及更高版本)：
http://docs.scipy.org/doc/numpy/reference/generated/numpy.interp.html

这正是我需要的,但我需要它在两个方面.我目前只是单步执行数组中的列并使用np.interp,但这当然很慢.

那里有什么可以实现同样的结果但更快？

最佳答案 解释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

现在,在所有类型检查发生之后,我从源代码中的句点！=无分支开始：

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

这只是确保提供的x和xp的所有值都在0和周期之间.因此,由于周期是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执行线性插值.当试图在xp中存在的两个点a和b之间进行插值时,它仅使用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来获得结果.让我们得到所有y的x = pi / 8的值：

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,句号]的值,你需要自己模数：

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])

同样,如果要为一堆x值和y值生成interp_points,请查看this answer.

点赞