assert y.shape[:2] == new_x.shape and x.shape == y.shape[2:]
nx, ny = y.shape[:2]
new_x = new_x.ravel()
j = np.arange(len(new_x))
k = np.searchsorted(x, new_x).clip(1, len(x) - 1)
y = y.reshape(-1, x.shape[0])
p = (new_x - x[k - 1]) / (x[k] - x[k - 1])
result = (1 - p) * y[j, k - 1] + p * y[j, k]
return result.reshape(nx, ny)
Fast interpolation over 3D array