添加具有不同轴数组的 2D numpy 数组：如何用 RectBivariateSpline 正确替换已弃用的 interp2d？

Question

我需要添加两个可能不同形状的 2D numpy 数组和不同的对应轴数组。

我的意思是这样的：让我们定义两组不同的 x 轴和 y 轴，并根据某个函数计算 z 值（此处：二维高斯分布）：

import numpy as np
import matplotlib.pyplot as plt


def gaussian_2D(X, Y, amplitude, mu_x, mu_y, sigma_x, sigma_y, theta, offset=0):
    """
    X,Y: are expected to be numpy meshgrids
    """
    a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)
    b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)
    c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)
    return offset + amplitude*np.exp( - (a*((X-mu_x)**2) + 2*b*(X-mu_x)*(Y-mu_y) + c*((Y-mu_y)**2)))


x1 = np.linspace(10, 100, num=100)
y1 = np.linspace(0, 200, num=100)
X1,Y1 = np.meshgrid(x1,y1)
Z1 = gaussian_2D(X1,Y1, 10, 35, 100, 10, 20, 0)

x2 = np.linspace(0, 150, num=120)
y2 = np.linspace(30, 220, num=120)
X2,Y2 = np.meshgrid(x2,y2)
Z2 = gaussian_2D(X2,Y2, 10, 75, 150, 5, 4, 12)

上面的代码生成 2 个长度为 (100,) 和 (120,) 的不同 x 数组，2 个长度为 (100,) 和 (120,) 的不同 y 数组以及 2 个形状为 (100,100) 和 ( 120,120).

x 和 y 数组定义了一些物理尺寸，这里的空间以毫米为单位，因此将两者相加是有意义的，即使底层数组不同。我选择这些轴重叠，但实际上没有必要，数组可以完全分开。绘制它们可以显示正在发生的事情：

fig1, ax1 = plt.subplots()
ax1.pcolormesh(X1,Y1,Z1)
ax1.set_xlim([0,150])
ax1.set_ylim([0,220])
ax1.set_xlabel("x [mm]")
ax1.set_ylabel("y [mm]")

fig2, ax2 = plt.subplots()
ax2.pcolormesh(X2,Y2,Z2)
ax2.set_xlim([0,150])
ax2.set_ylim([0,220])
ax2.set_xlabel("x [mm]")
ax2.set_ylabel("y [mm]")

结果：

我们可以通过简单地将两个数组用一些

alpha

值绘制在彼此之上来了解添加两个数组时结果图应该是什么样子：

fig3, ax3 = plt.subplots()
ax3.pcolormesh(X1, Y1, Z1, alpha=0.5)
ax3.pcolormesh(X2, Y2, Z2, alpha=0.3)
ax3.set_xlabel("x [mm]")
ax3.set_ylabel("y [mm]")

---

现在，显然因为数组的底层结构不同，我们必须使用一些插值将两个数组“标准化”到某个公共网格上。为了方便起见，让结果数组始终为 (100, 100) （即与两个输入数组的形状大致相同。这是一个通常对我来说成立的假设......让我们忽略边缘情况...... .)

def add_arrays_with_axes(array_1, array_1_x, array_1_y, array_2, array_2_x, array_2_y, method="interp2d"):
    import scipy.interpolate as interp

    # Define the new x and y axis ranges that cover both array_1 and array_2 axes
    new_x = np.linspace(min(array_1_x[0], array_2_x[0]), max(array_1_x[-1], array_2_x[-1]), num=100)
    new_y = np.linspace(min(array_1_y[0], array_2_y[0]), max(array_1_y[-1], array_2_y[-1]), num=100)
    
    # Interpolate array_1 and array_2 onto the new grid
    match method:
        case "interp2d":
            interp_array_1 = interp.interp2d(array_1_x, array_1_y, array_1, kind='cubic', bounds_error=False, fill_value=np.nan)
            interp_array_2 = interp.interp2d(array_2_x, array_2_y, array_2, kind='cubic', bounds_error=False, fill_value=np.nan)

        case "RectBivariateSpline":
            interp_array_1 = interp.RectBivariateSpline(array_1_x, array_1_y, array_1)
            interp_array_2 = interp.RectBivariateSpline(array_2_x, array_2_y, array_2)

        case _:
            raise Exception("Unsupported interp method.")
    
    # Evaluate the interpolations on the new x and y grid
    array_1_resampled = interp_array_1(new_x, new_y)
    array_2_resampled = interp_array_2(new_x, new_y)
    
    # Replace NaNs with 0 in each array where data was not originally defined
    array_1_resampled = np.nan_to_num(array_1_resampled, nan=0.0)
    array_2_resampled = np.nan_to_num(array_2_resampled, nan=0.0)
    
    # Sum the resampled arrays
    summed_array = array_1_resampled + array_2_resampled
    
    return summed_array, new_x, new_y

在上面我允许使用两种不同的插值方法：

interp2d

和

RectBivariateSpline

。

让我们看看结果是什么：

interp2d

new_array, new_x, new_y = add_arrays_with_axes(Z1, x1, y1, Z2, x2, y2, method="interp2d")
new_mesh = np.meshgrid(new_x, new_y)

fig4, ax4 = plt.subplots()
ax4.pcolormesh(X1, Y1, Z1, alpha=0.5)
ax4.pcolormesh(X2, Y2, Z2, alpha=0.3)
ax4.set_xlabel("x [mm]")
ax4.set_ylabel("y [mm]")

ax4.pcolormesh(*new_mesh, new_array)

---

矩形双变量样条线

new_array, new_x, new_y = add_arrays_with_axes(Z1, x1, y1, Z2, x2, y2, method="RectBivariateSpline")
new_mesh = np.meshgrid(new_x, new_y)

fig5, ax5 = plt.subplots()
ax5.pcolormesh(X1, Y1, Z1, alpha=0.5)
ax5.pcolormesh(X2, Y2, Z2, alpha=0.3)
ax5.set_xlabel("x [mm]")
ax5.set_ylabel("y [mm]")

ax5.pcolormesh(*new_mesh, new_array)

显然这两种方法显示出截然不同的结果

interp2d

展示了我期望发生的情况：您只能在定义原始数组的地方获得值，并在未定义的地方获得零，从而导致

x=10

处的可见边缘。 另一方面，

RectBivariateSpline

似乎不仅拉伸了数据，而且还移动了数据......例如比较更广泛的高斯的中心点。它被定义为 (x=35, y=100)，并且

interp2d

正确地再现了这一点，但在

RectBivariateSpline

情况下，中心移动到更像 (x=35, y=80) 的位置。 为什么要这样做？ 看起来

interp2d

是可行的方法，但该功能实际上已已停止。

RectBivariateSpline

确实还有一个可能相关的参数：

bbox

，但我不确定它的作用，也不确定我会将其设置为除原始边界框之外的原始数组，无论如何它都默认为..另外，这个论点似乎存在“一些问题”。 知道

RectBivariateSpline

发生了什么吗？

Answer 1

__call__

 的 

RectBivariateSpline

 方法，我们发现代码需要 

indexing="ij" 格式的点，而 np.meshgrid

 默认为 

indexing="xy"。要解决此问题，您需要进行两项更改：

创建

RectBivariateSpline
1. 对象时转置数据数组。

case "RectBivariateSpline":
    interp_array_1 = interp.RectBivariateSpline(array_1_x, array_1_y, array_1.T)
    interp_array_2 = interp.RectBivariateSpline(array_2_x, array_2_y, array_2.T)

使用

new_mesh
2. 订购创建您的

   "ij"

   。
new_mesh = np.meshgrid(new_x, new_y, indexing="ij")

这将产生与 

interp2d

的匹配结果。

您还应该记住，

RectBivariateSpline

会进行推断，而

interp2d

会填充未知值（在您的情况下，您将使用

np.nan

填充）。仅当您在一个数据集的边界附近/处有某种结构并且该边界在两个数据集之间共享时，这才会成为问题。