我有一个包含一些给定数据的网格。该数据由其角度给出(从
0
到π
)。
在这个网格中,我有另一个较小的网格。
这可能看起来像这样:
现在我想在该网格上插入角度。
我通过使用
scipy.interpolate.griddata
来尝试这个,结果很好。但是当角度从几乎0
变成几乎π
时就有问题了(因为中间是π/2
……)
这是结果,很容易看出哪里出了问题。
我该如何处理这个问题?谢谢你! :)
这里是重现的代码:
import numpy as np
from matplotlib import pyplot as plt
from scipy.interpolate import griddata
ax = plt.subplot()
ax.set_aspect(1)
# Simulate some given data.
x, y = np.meshgrid(np.linspace(-10, 10, 20), np.linspace(-10, 10, 20))
data = np.arctan(y / 10) % np.pi
u = np.cos(data)
v = np.sin(data)
ax.quiver(x, y, u, v, headlength=0.01, headaxislength=0, pivot='middle', units='xy')
# Create a smaller grid within.
x1, y1 = np.meshgrid(np.linspace(-1, 5, 15), np.linspace(-6, 2, 20))
# ax.plot(x1, y1, '.', color='red', markersize=2)
# Interpolate data on grid.
interpolation = griddata((x.flatten(), y.flatten()), data.flatten(), (x1.flatten(), y1.flatten()))
u1 = np.cos(interpolation)
v1 = np.sin(interpolation)
ax.quiver(x1, y1, u1, v1, headlength=0.01, headaxislength=0, pivot='middle', units='xy',
color='red', scale=3, width=0.03)
plt.show()
编辑:
感谢@bubble,有一种方法可以在插值之前调整给定的角度,从而使结果符合要求。 因此:
定义一个整流函数:
def RectifyData(data):
for j in range(len(data)):
step = data[j] - data[j - 1]
if abs(step) > np.pi / 2:
data[j] += np.pi * (2 * (step < 0) - 1)
return data
插值如下:
interpolation = griddata((x.flatten(), y.flatten()),
RectifyData(data.flatten()),
(x1.flatten(), y1.flatten()))
u1 = np.cos(interpolation)
v1 = np.sin(interpolation)
我尝试直接插值
cos(angle)
和 sin(angle)
值,但这仍然会中断,导致错误的线方向。主要思想在于减少中断,例如[2.99,3.01, 0.05,0.06]
应该转换为这样的东西:[2.99, 3.01, pi+0.05, pi+0.06]
。这是正确应用 2D 插值算法所必需的。几乎相同的问题出现在以下post.
def get_rectified_angles(u, v):
angles = np.arcsin(v)
inds = u < 0
angles[inds] *= -1
# Direct approach of removing discontinues
# for j in range(len(angles[1:])):
# if abs(angles[j] - angles[j - 1]) > np.pi / 2:
# sel = [abs(angles[j] + np.pi - angles[j - 1]), abs(angles[j] - np.pi - angles[j-1])]
# if np.argmin(sel) == 0:
# angles[j] += np.pi
# else:
# angles[j] -= np.pi
return angles
ax.quiver(x, y, u, v, headlength=0.01, headaxislength=0, pivot='middle', units='xy')
# # Create a smaller grid within.
x1, y1 = np.meshgrid(np.linspace(-1, 5, 15), np.linspace(-6, 2, 20))
angles = get_rectified_angles(u.flatten(), v.flatten())
interpolation = griddata((x.flatten(), y.flatten()), angles, (x1.flatten(), y1.flatten()))
u1 = np.cos(interpolation)
v1 = np.sin(interpolation)
ax.quiver(x1, y1, u1, v1, headlength=0.01, headaxislength=0, pivot='middle', units='xy',
color='red', scale=3, width=0.03)
可能,
numpy.unwrap
功能可用于修复中断。对于一维数据,numpy.interp
有关键字 period
来处理周期性数据。
还可以使用
cmath
模块将角度转换为复数,进行插值,然后在插值后转换回角度。
import cmath
# represent angles with unit vectors using complex numbers (radius=1, angle=angles).
cmplx = cmath.rect( 1+ 0*angles, angles)
interpolation = griddata((x.flatten(), y.flatten(), cmplx, (x1.flatten(), y1.flatten()))
# imaginary part of the complex number is the angle
angles_interpolated = interpolation.imag