Мне нужно сложить два двумерных массива 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: должны быть сетками numpy
    """
    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 различных массива x длиной (100,) и (120,), 2 различных массива y длиной (100,) и (120,) и 2 массива z формы (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 [мм]")
ax1.set_ylabel("y [мм]")

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

результатом будет:

И мы можем получить представление о том, как должен выглядеть результирующий график при сложении двух массивов, просто наложив их друг на друга с некоторым значением 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 [мм]")
ax3.set_ylabel("y [мм]")

Теперь, очевидно, поскольку структура исходных массивов различна, мы должны использовать некоторую интерполяцию, чтобы “стандартизировать” два массива на общей сетке. Для удобства сделаем так, чтобы результирующий массив всегда имел форму (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

    # Определим новые диапазоны осей x и y, которые охватывают обе оси array_1 и array_2
    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)

    # Интерполируем array_1 и array_2 на новой сетке
    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("Неподдерживаемый метод интерполяции.")

    # Оцениваем интерполяции на новой сетке x и y
    array_1_resampled = interp_array_1(new_x, new_y)
    array_2_resampled = interp_array_2(new_x, new_y)

    # Заменяем NaN на 0 в каждом массиве, где данные изначально не определены
    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)

    # Складываем переинтерполированные массивы
    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 [мм]")
ax4.set_ylabel("y [мм]")

ax4.pcolormesh(*new_mesh, new_array)

RectBivariateSpline

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 [мм]")
ax5.set_ylabel("y [мм]")

ax5.pcolormesh(*new_mesh, new_array)

Очевидно, что два метода показывают очень разные результаты

interp2d показывает именно то, что я ожидал: значения получаются только там, где исходные массивы были определены, и нули там, где они не были определены, что приводит к видимому краю при x=10.
RectBivariateSpline, с другой стороны, похоже, не только растягивает данные, но и изменяет их… Сравните, например, центральную точку более широкой гауссовой функции. Она была определена как (x=35, y=100), и interp2d воспроизводит это правильно, но в случае RectBivariateSpline центр сдвинулся куда-то ближе к (x=35, y=80). Почему это происходит? Кажется, что interp2d — это правильный путь, но эта функция на самом деле устарела. RectBivariateSpline имеет еще один, возможно, релевантный аргумент: bbox, но я не уверен, что он делает, и что мне следует установить, кроме оригинального ограничивающего прямоугольника исходного массива, который у него по умолчанию… Плюс к тому же, кажется, есть некоторая проблема с этим аргументом.

Есть ли у вас идеи относительно того, что происходит с RectBivariateSpline?

Посмотрев на __call__ метод RectBivariateSpline, мы видим, что код ожидает точки в формате indexing="ij", в то время как np.meshgrid по умолчанию использует indexing="xy" (не спрашивайте меня, почему). Чтобы исправить это, вам нужно будет внести два изменения:

1. Транспонировать массивы данных при создании объектов RectBivariateSpline.

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)

2. Создать new_mesh с использованием порядка "ij".

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

Это даст результаты, сопоставимые с interp2d.

Также следует учитывать, что RectBivariateSpline экстраполирует, в то время как interp2d заполняет неизвестные значения (в вашем случае вы заполняете np.nan). Это станет проблемой только в том случае, если у вас есть какая-либо структура около/на границе одного набора данных, и эта граница является общей для двух наборов данных.