У меня есть несколько 2D-массивов данных, которые мне нужно интерполировать. SciPy предлагает некоторые решения для этого, как хорошо объяснено, например, здесь . Мои данные очень грубо выглядят как двумерная парабола (в зависимости от x1 и x2) с отверстием, направленным вниз, как показано на рис. 1 для горизонтального разреза вдоль x2 = 0. Как видите, отрицательных значений нет (все точки данных равны нулю).
Я хотел выполнить кубическую c интерполяцию как мне требуется гладкие данные. Это создает проблему на краях, что приводит к «покачиванию» или «перерегулированию» подгонки / интерполяции, как показано на рисунке 2. Однако отрицательные значения не допускаются в последующей последующей обработке интерполированных данных (также следующее подавление положительных значений, где они должны быть равны нулю, должно быть подавлено).
.
Я думал, что разумным «решением» было бы просто установить значения, которые равны 0 (обратите внимание, что все они в точности совпадают), равными NaN, так что они игнорируются интерполяцией. Но использование SciPy griddata и метод cubic не работают с NaN с. Метод linear может справиться с этим, но мне нужно cubic.
Мой вопрос: я что-то упускаю или делаю что-то неправильно, что приводит к griddata неправильной работе с NaN s и методу cubic?
Пример кода такой следует:
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate as interp
def some_funct( x1, x2 ):
result = -x1**2 - 2.*x2**2 + 60.
result[result < .0] = .0
return result
# define original (sparse) grid
N_x1, N_x2 = 20, 20
x1_old = np.linspace( -9, 10, N_x1 )
x2_old = np.linspace( -9, 10, N_x2 )
X1_old, X2_old = np.meshgrid( x1_old, x2_old )
# datapoints to be interpolated
z_old = some_funct( X1_old, X2_old )
# set 0 datapoints to nan
z_old[z_old==0] = np.nan
# grid for interpolation
x1_new = np.linspace( -9, 10, 10*N_x1 )
x2_new = np.linspace( -9, 10, 10*N_x2 )
X1_new, X2_new = np.meshgrid( x1_new, x2_new )
# perform interpolation
z_new = interp.griddata( np.array([X1_old.ravel(),X2_old.ravel()]).T, z_old.ravel(),
(X1_new, X2_new),
method='cubic', fill_value=.0 # only works for 'linear'
)
# plot horizonal cut along x2=0
fig1 = plt.figure( figsize=(8,6) )
ax1 = fig1.add_subplot( 1,1,1 )
x2_old_0_id = (np.abs(x2_old - .0)).argmin()
x2_new_0_id = (np.abs(x2_new - .0)).argmin()
ax1.plot( x1_old, z_old[ x2_old_0_id , : ], marker='x', linestyle='None', label='data' )
ax1.plot( x1_new, z_new[ x2_new_0_id , : ], label='interpolation' )
ax1.legend()
ax1.set_xlabel( 'x1' )
ax1.set_ylabel( 'z' )
plt.show()
Любые советы приветствуются!
Обновление: Забыл указать версии, которые я использую:
numpy: 1.15.1
scipy: 1.1.0