【问题标题】:2D linear interpolation: data and interpolated points2D 线性插值:数据和插值点
【发布时间】:2017-07-05 09:51:39
【问题描述】:

考虑这个 y(x) 函数:

我们可以在文件中生成这些散点的位置:dataset_1D.dat:

# x   y   
0   0
1   1
2   0
3   -9
4   -32

以下是这些点的一维插值代码:

  1. 加载这个散点

  2. 创建x_mesh

  3. 执行一维插值

代码:

import numpy as np
from scipy.interpolate import interp2d, interp1d, interpnd
import matplotlib.pyplot as plt


# Load the data:    
x, y  = np.loadtxt('./dataset_1D.dat', skiprows = 1).T

# Create the function Y_inter for interpolation:
Y_inter = interp1d(x,y)

# Create the x_mesh:    
x_mesh = np.linspace(0, 4, num=10)
print x_mesh

# We calculate the y-interpolated of this x_mesh :   
Y_interpolated = Y_inter(x_mesh)
print Y_interpolated

# plot:

plt.plot(x_mesh, Y_interpolated, "k+")
plt.plot(x, y, 'ro')
plt.legend(['Linear 1D interpolation', 'data'], loc='lower left',  prop={'size':12})
plt.xlim(-0.1, 4.2)
plt.grid()
plt.ylabel('y')
plt.xlabel('x')
plt.show()

这绘制了以下内容:

现在,考虑这个 z(x,y) 函数:

我们可以在文件中生成这些散点:dataset_2D.dat

# x    y    z
0   0   0
1   1   0
2   2   -4
3   3   -18
4   4   -48

在这种情况下,我们必须执行 2D 插值:

import numpy as np
from scipy.interpolate import interp1d, interp2d, interpnd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Load the data:
x, y, z  = np.loadtxt('./dataset_2D.dat', skiprows = 1).T

# Create the function Z_inter for interpolation:
Z_inter = interp2d(x, y, z)

# Create the x_mesh and y_mesh :
x_mesh = np.linspace(1.0, 4, num=10)
y_mesh = np.linspace(1.0, 4, num=10)
print x_mesh
print y_mesh

# We calculate the z-interpolated of this x_mesh and y_mesh :
Z_interpolated = Z_inter(x_mesh, y_mesh)
print Z_interpolated
print type(Z_interpolated)
print Z_interpolated.shape

# plot: 
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(x, y, z, c='r', marker='o')
plt.legend(['data'], loc='lower left',  prop={'size':12})
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

plt.show()

这绘制了以下内容:

分散的数据再次以红点显示,与二维图保持一致。

  1. 我不知道如何解释Z_interpolated结果:

    根据上述代码的打印行, Z_interpolated 是一个 n 维 numpy 数组,形状为 (10,10)。换句话说,一个 10 行 10 列的二维矩阵。

我希望 x_mesh[i]y_mesh[i] 的每个值都有一个插入的 z[i] 值,为什么我没有收到这个?

  1. 如何在 3D 图中也绘制插值数据(就像 2D 图中的黑色十字)?

【问题讨论】:

    标签: python numpy matplotlib linear-interpolation mplot3d


    【解决方案1】:

    Z_interpolated的解释:你的一维x_meshy_mesh定义了一个mesh on which to interpolate。因此,您的二维插值返回 z 是一个二维数组,其形状 (len(y), len(x)) 与 np.meshgrid(x_mesh, y_mesh) 匹配。如您所见,您的 z[i, i] 是 x_mesh[i]y_mesh[i] 的预期值,而不是 z[i]。它还有更多,网格上的所有值。

    显示所有插值数据的潜在图:

    from mpl_toolkits.mplot3d import Axes3D
    import matplotlib.pyplot as plt
    import numpy as np
    from scipy.interpolate import interp2d
    
    # Your original function
    x = y = np.arange(0, 5, 0.1)
    xx, yy = np.meshgrid(x, y)
    zz = 2 * (xx ** 2) - (xx ** 3) - (yy ** 2)
    
    # Your scattered points
    x = y = np.arange(0, 5)
    z = [0, 0, -4, -18, -48]
    
    # Your interpolation
    Z_inter = interp2d(x, y, z)
    x_mesh = y_mesh = np.linspace(1.0, 4, num=10)
    Z_interpolated = Z_inter(x_mesh, y_mesh)
    
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    # Plot your original function
    ax.plot_surface(xx, yy, zz, color='b', alpha=0.5)
    # Plot your initial scattered points
    ax.scatter(x, y, z, color='r', marker='o')
    # Plot your interpolation data
    X_real_mesh, Y_real_mesh = np.meshgrid(x_mesh, y_mesh)
    ax.scatter(X_real_mesh, Y_real_mesh, Z_interpolated, color='g', marker='^')
    plt.show()
    

    【讨论】:

    • 非常感谢您的回答。在我的帖子中,我使用interp2d 来使用更好的三次插值。通过在您的代码中实现这一点,Z_inter = interp2d(x, y, z, kind='cubic') 它不起作用,您知道为什么吗?谢谢
    • 在我的示例代码中仅更改行 Z_inter = interp2d(x, y, z) 是行不通的。根据scipy.interpolate.interp2d()的“注释”,“沿插值轴所需的最小数据点数为(k+1)**2,其中k=1表示线性,k=3表示三次,k=5表示五次插值”。为了使它工作,将xy 更改为x = y = np.arange(0, 20) 并将z 更改为z = 2 * (x ** 2) - (x ** 3) - (y ** 2)。如果这不是原因,请告诉我您遇到了什么错误。
    【解决方案2】:

    您需要两步插值。第一个在 y 数据之间进行插值。第二个在 z 数据之间进行插值。然后用两个插值数组绘制x_mesh

    x_mesh = np.linspace(0, 4, num=16)
    
    yinterp = np.interp(x_mesh, x, y)
    zinterp = np.interp(x_mesh, x, z)
    
    ax.scatter(x_mesh, yinterp, zinterp, c='k', marker='s')
    

    在下面的完整示例中,我还在 y 方向添加了一些变化,以使解决方案更通用。

    u = u"""# x    y    z
    0   0   0
    1   3   0
    2   9   -4
    3   16   -18
    4   32   -48"""
    
    import io
    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    
    # Load the data:
    x, y, z  = np.loadtxt(io.StringIO(u), skiprows = 1, unpack=True)
    
    x_mesh = np.linspace(0, 4, num=16)
    
    yinterp = np.interp(x_mesh, x, y)
    zinterp = np.interp(x_mesh, x, z)
    
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.scatter(x_mesh, yinterp, zinterp, c='k', marker='s')
    ax.scatter(x, y, z, c='r', marker='o')
    plt.legend(['data'], loc='lower left',  prop={'size':12})
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    
    plt.show()
    

    对于使用scipy.interpolate.interp1d,解决方案基本相同:

    u = u"""# x    y    z
    0   0   0
    1   3   0
    2   9   -4
    3   16   -18
    4   32   -48"""
    
    import io
    import numpy as np
    from scipy.interpolate import interp1d
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    
    # Load the data:
    x, y, z  = np.loadtxt(io.StringIO(u), skiprows = 1, unpack=True)
    
    x_mesh = np.linspace(0, 4, num=16)
    
    fy = interp1d(x, y, kind='cubic')
    fz = interp1d(x, z, kind='cubic')
    
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.scatter(x_mesh, fy(x_mesh), fz(x_mesh), c='k', marker='s')
    ax.scatter(x, y, z, c='r', marker='o')
    plt.legend(['data'], loc='lower left',  prop={'size':12})
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    
    plt.show()
    

    【讨论】:

    • 非常感谢您的回答。不幸的是,np.interp 只允许线性插值。这就是我使用scipyinterp1d 的原因,以便更好地适应z(x,y),它在x 中是立方的。按照您的指示,并将其翻译为interp1d,我们将:1)在y 数据之间进行插值:Y_inter = interp1d(x,y, kind='cubic') 和2)在z 数据之间进行插值:Z_inter = interp1d(x,z, kind='cubic')。然后,我们可以做Z_interpolated = Z_inter(x_mesh, y_mesh)(我的帖子),但这并没有给出正确的答案......我们怎么能做到这一点?谢谢
    • 谢谢,我正在研究这个
    猜你喜欢
    • 1970-01-01
    • 2011-05-20
    • 2016-11-11
    • 1970-01-01
    • 2013-07-27
    • 1970-01-01
    • 2017-02-22
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多