稀疏(规则)数据的 matplotlib 轮廓显示人工制品

matplotlib contour of sparse (regular) data shows artefacts

我想绘制非常稀疏且最大值沿对角线穿过图片的数据轮廓; matplotlib 轮廓函数在采样的最大值之间发明了最小值。

从密集采样的情况开始,一切看起来都符合预期:

import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np

x_1D = np.linspace(0., 10., 100)
y_1D = np.linspace(0., 10., 100)
x, y = np.meshgrid(x_1D, y_1D)
z = np.empty_like(x)
def peak(y, y0):
    return np.exp(-(y-y0)**2)
for i in range(x_1D.size):
    z[:,i] = peak(y_1D, i/x_1D.size*y_1D.max())

fig, ax = plt.subplots(ncols=3)

ax[0].set_title('measured data')
ax[0].scatter(x, y, marker='s', c=z, cmap=plt.cm.jet, s=25)

ax[1].set_title('contour')
ax[1].contourf(x, y, z, levels=14, cmap=plt.cm.jet)

# define grid
xi = np.linspace(x_1D.min()-0.1, x_1D.max()+0.1, 1000)
yi = np.linspace(y_1D.min()-0.1, y_1D.max()+0.1, 1000)

# grid the data
triang = tri.Triangulation(x.flatten(), y.flatten())
interpolator = tri.LinearTriInterpolator(triang, z.flatten())
Xi, Yi = np.meshgrid(xi, yi)
zi = interpolator(Xi, Yi)

ax[2].set_title('interpolated')
ax[2].contourf(xi, yi, zi, levels=14, cmap=plt.cm.jet)

plt.show()

产量

当 x 的样本减少 10 倍时,即 x_1D = np.linspace(0., 10., 10),最小值出现在等高线图中的样本最大值之间。

有没有办法避免这种伪像并使稀疏采样数据的轮廓看起来像密集采样数据之一?

编辑:感谢您提供的答案在我提供的示例中效果很好。不幸的是,我把问题简化得太过了。我不应该谈论一条对角线,我应该询问在图片中向任意方向移动的任意数量的峰;例如用

替换峰值生成
z = np.zeros_like(x)
def peak(y, y0):
    return np.exp(-(y-y0)**2)
for i in range(x_1D.size):
    z[:,i] += peak(y_1D, np.cos(i/x_1D.size*np.pi)*y_1D.max()*0.05+y_1D.max()*0.8)
for i in range(x_1D.size):
    z[:,i] += peak(y_1D, np.sin(i/x_1D.size*np.pi/2.)*y_1D.max()*0.5)

导致

您的方法的主要问题是三角测量算法不知道峰值应该在 "x-slices"(您的常数 x 的密集数据点线)之间相互连接。

稍微简化一下,三角剖分算法将查看 x 和 y 方向上的邻居并连接到那些邻居。然后,当尝试使用此三角测量进行插值时,峰之间的点将大致是 x 方向上最近点的平均值,因此会出现最小值。最好的解决方案是进行自己的三角测量,将峰直接连接起来。

幸运的是,我们实际上可以破解三角测量,通过在 y 方向上移动坐标使其连接到峰,从而使峰全部水平对齐。这是可行的,因为三角测量算法使用您传递给它的坐标。在您的示例中,这很容易实现,因为我们只需应用简单的移位 y_s = y - x。通常,您必须获得峰值的方程式(称之为 y_p(x)),然后从 y 中减去它以获得 y_s

现在您已经有了移动的三角测量,您可以制作一个新的更密集的网格(就像您所做的那样)并应用相同的移动。然后,您使用移位的密集网格在移位的网格中进行插值,以获得正确插值的 z 值。最后,您 un-shift 密集网格以获得正确的 y 值并绘制它。

以下是将此概念应用到您的代码和最终结果的代码。如你看到的。它在这种情况下效果很好。

import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np

def peak(y, y0):
    return np.exp(-(y-y0)**2)

x_1D = np.linspace(0., 10., 10)
y_1D = np.linspace(0., 10., 100)
x, y = np.meshgrid(x_1D, y_1D)
z = np.empty_like(x)

for i in range(x_1D.size):
    z[:,i] = peak(y_1D, i/x_1D.size*y_1D.max())

fig, ax = plt.subplots(ncols=4)

ax[0].set_title('measured data')
ax[0].scatter(x, y, marker='s', c=z, cmap=plt.cm.jet, s=25)

ax[1].set_title('contour')
ax[1].contourf(x, y, z, levels=14, cmap=plt.cm.jet)

# define output grid
xi_1D = np.linspace(x_1D.min()-0.1, x_1D.max()+0.1, 1000)
yi_1D = np.linspace(y_1D.min()-0.1, y_1D.max()+0.1, 1000)
xi, yi = np.meshgrid(xi_1D, yi_1D)

# Old Linear Interpolation
triang = tri.Triangulation(x.flatten(), y.flatten())
interpolator = tri.LinearTriInterpolator(triang, z.flatten())
zi = interpolator(xi, yi)

ax[2].set_title('interpolated')
ax[2].contourf(xi, yi, zi, levels=14, cmap=plt.cm.jet)

# === SHIFTED LINEAR INTERPOLATION ===

# make shifted interpolating mesh for the data
y_s=y-x
triang_s = tri.Triangulation(x.flatten(), y_s.flatten())
interpolator_s = tri.LinearTriInterpolator(triang_s, z.flatten())

# interpolate in the shifted state
yi_s = yi-xi
zi_s = interpolator_s(xi, yi_s)

# unshift the fine mesh
yi_us = yi_s+xi

ax[3].set_title('interpolated (shifted)')
ax[3].contourf(xi, yi_us, zi_s, levels=14, cmap=plt.cm.jet)


plt.show()