I need a faster solution to interpolate a 3D array in one dimension to a 2D surface. My data is gridded (k, j, i): the horizontal grid is regular (j, i), but the vertical grid is not (i.e depth points are different at each i, j point). Additionally, my array has a 2D horizontal mask. My current solution is:
from scipy.interpolate import interp1d
def interpz(depths, var, z_levs):
var_z = np.zeros_like(var[0,:,:])
for j in range(var_z.shape[0]):
for i in range(var_z.shape[1]):
if var_z.mask[j,i]==False:
f = interp1d(depths[:,j,i], var[:,j,i])
var_z[j,i] = f(z_levs[j,i])
return var_z
It works, but it's too slow (my horizontal grid is about 200x600, and the number of interpolations reduces quite a lot with the mask but it still takes about one hour).
I've looked into griddata, but there is one issue, and one problem: the issue, it uses all the values to interpolate the point, when I actually only require 1D interpolation in the z direction (which, whatever), BUT the problem is that apparently, it doesn't like nans (in fact. not actual NaNs, but masked values):
z, y, x = depths.flatten(), y.flatten(), x.flatten()
data = var.flatten()
var_z = griddata((z, y, x), data, (z_levs.flatten(), y[0,:,:].flatten(),
...: x[0].flatten()), method='linear')
I get:
in griddata(points, values, xi, method, fill_value, rescale)
215 elif method == 'linear':
216 ip = LinearNDInterpolator(points, values, fill_value=fill_value,
--> 217 rescale=rescale)
218 return ip(xi)
219 elif method == 'cubic' and ndim == 2:
scipy/interpolate/interpnd.pyx in scipy.interpolate.interpnd.LinearNDInterpolator.__init__ (scipy/interpolate/interpnd.c:5530)()
scipy/spatial/qhull.pyx in scipy.spatial.qhull.Delaunay.__init__ (scipy/spatial/qhull.c:18174)()
scipy/spatial/qhull.pyx in scipy.spatial.qhull._Qhull.__init__ (scipy/spatial/qhull.c:4788)()
ValueError: Points cannot contain NaN
Suggestions are deeply appreciated