How to improve performance when interpolating on 3d data with SciPy

Question

Hello there,

I have 3d-data representing the atmosphere. Now I want to interpolate this data to a common Z coordinate (what I mean by that should be clear from the function's doctring). The following code works fine, but I was wondering if there were a way to improve the performance ...

def interpLevel(grid,value,data,interp='linear'):
    """
    Interpolate 3d data to a common z coordinate.

    Can be used to calculate the wind/pv/whatsoever values for a common
    potential temperature / pressure level.

    grid : numpy.ndarray
       The grid. For example the potential temperature values for the whole 3d
       grid.

    value : float
       The common value in the grid, to which the data shall be interpolated.
       For example, 350.0

    data : numpy.ndarray
       The data which shall be interpolated. For example, the PV values for
       the whole 3d grid.

    kind : str
       This indicates which kind of interpolation will be done. It is directly
       passed on to scipy.interpolate.interp1d().

    returs : numpy.ndarray
       A 2d array containing the *data* values at *value*.

    """
    ret = np.zeros_like(data[0,:,:])
    # we need to copy the grid to a new one, because otherwise the flipping
    # done below will be messed up
    gr = np.zeros_like(grid)
    da = np.zeros_like(data)
    for latIdx in xrange(grid.shape[1]):
        for lonIdx in xrange(grid.shape[2]):
            # check if we need to flip the column
            if grid[0,latIdx,lonIdx] > grid[-1,latIdx,lonIdx]:
                gr[:,latIdx,lonIdx] = grid[::-1,latIdx,lonIdx]
                da[:,latIdx,lonIdx] = data[::-1,latIdx,lonIdx]
            else:
                gr[:,latIdx,lonIdx] = grid[:,latIdx,lonIdx]
                da[:,latIdx,lonIdx] = data[:,latIdx,lonIdx]
            f = interpolate.interp1d(gr[:,latIdx,lonIdx], \
                    da[:,latIdx,lonIdx], \
                    kind=interp)
            ret[latIdx,lonIdx] = f(value)
    return ret

Answer 1

Well, this might give a small speed-up just because it uses less memory.

ret = np.zeros_like(data[0,:,:])
for latIdx in xrange(grid.shape[1]):
    for lonIdx in xrange(grid.shape[2]):
        # check if we need to flip the column
        if grid[0,latIdx,lonIdx] > grid[-1,latIdx,lonIdx]:
            ind = -1
        else:
            ind = 1
        f = interpolate.interp1d(grid[::ind,latIdx,lonIdx], \
                data[::ind,latIdx,lonIdx], \
                kind=interp)
        ret[latIdx,lonIdx] = f(value)
return ret

All I've done is get rid of gr and da really.

Other than that, are you calling this function with a whole lot of different values(i.e. value being different but other parameters the same)? If so, you might want to make the function be able to handle multiple values (add another dimension to ret in other words that is as long as the length of values). Then you are making better use of the interpolation function that you've created.

The last suggestion is to try a profiler. It will allow you to see what is taking the most time.