interpolate time series of 2d gridded data to a new 2d grid - efficiently

Viewed 39

I have an numpy array of 2D gridded data with a temporal axis, so my array has a shape of (nsteps, ny, nx)

I'm trying to interpolate data from this grid to a very slightly different grid (different resolution and thus node points).

I was able to do this fine via:

import numpy as np
from scipy.interpolate import RectBivariateSpline

#some example arrays
p_dat = np.random.random((10, 182, 361)) #old grid, 182rows, 361cols
w_dat = np.random.random((10, 200, 400)) #new grid, 200rows, 400cols

#the grids
x0  = np.linspace(0, 360, 361) #old
y0  = np.linsapce(-90, 90, 182) #old
x   = np.linspace(0, 360, 400) #new
y   = np.linspace(-90, 90 , 200) #new

#new array with 2d shape of w_dat
out = np.full((10, 200, 400), np.nan)

#interpolate one timestep at a time
for i in range(out.shape[0]):
    interp = RectBivariateSpline(y0, x0, p_dat[i])
    dat = interp(y, x)
    out[i,:,:] = dat

Is there a way I can avoid this loop and vectorize this interpolation over the 0th axis?

1 Answers

Do you have any data you could add to your question? That would enable us to have a closer look.

I think np.apply_along_axis might work here.

Right now i can only provide this untested approach:

def helper(idx: int):
    interp = RectBivariateSpline(pres.geometry.y, pres.geometry.x, p_dat[idx])
    dat = interp(wind.geometry.y, wind.geometry.x)
    return dat

out = np.apply_along_axis(helper, axis=0, arr=np.arange(out.shape[0]))

Edit: I can't comment on speedups due to the lack of test data, but this is very similar to map(helper, np.arange(out.shape[0])). In essence, it's just hiding the loop.

Related