Добавление двумерных массивов numpy с разными массивами осей: как правильно заменить устаревший interp2d на RectBivariat ⇐ Python

[Anonymous]

Мне нужно добавить два двумерных массива, возможно, разных форм и разных соответствующих массивов осей.
Я имею в виду следующее: давайте определим два разных набора осей 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 разных 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 [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]")

результат:




И мы можем получить представление о том, как должен выглядеть итоговый график при добавлении двух массивов, просто поместив их друг на друга с некоторым значением альфа:

Код: Выделить всё

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)

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 [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?

Подробнее здесь: https://stackoverflow.com/questions/791 ... lace-the-d