I have some 2D data-arrays which I need to interpolate. SciPy offers some solutions to this, as nicely explained for example here. My data very roughly looks like a 2D parabola (as a function of x1 and x2), with the opening pointing downwards, as illustrated in Fig. 1 for a horizontal cut along x2=0. As you can see, there are no negative values (the datapoints are all exactly 0 there).
I wanted to perform cubic interpolation as I require smooth data. This gives a problem at the edges, resulting in "wiggling" or "overshooting" of the fit/interpolation, as illustrated in Fig. 2. Negative values are, however, not allowed in the following post-processing of the interpolated data (also the following overshooting to positive values where they should be zero need to be suppressed).
.
I thought that a clever "solution" would be to simply set the values which are 0 (note that those are all exactly the same) to NaN, such that they are ignored by the interpolation. But using SciPy's griddata and the cubic method is not working with NaNs. The linear method can handle this, but I need cubic.
My question is, am I missing something or doing something wrong that results in griddata not working properly with NaNs and the cubic method ?
An example code is as follows:
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()
Any hints are greatly appreciated!
Update: Forgot to include the versions I am using:
numpy: 1.15.1
scipy: 1.1.0
