2010-02-22 63 views
3

我擁有代表大氣的3d數據。現在我想將這些數據插入到一個普通的Z座標中(我的意思是應該從函數的doctring中清楚)。下面的代碼工作正常,但我不知道是否有是提高性能的一種方式......如何在用SciPy插值3d數據時提高性能

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 

回答

2

好,只是因爲它使用較少的內存,這可能給小加速。

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 

我所做的就是擺脫gr和da真的。

除此之外,你是否用很多不同的值調用這個函數(即,值是不同的,但其他參數相同)?如果是這樣,你可能想讓這個函數能夠處理多個值(添加另一個維度來進行ret,換句話說就是值的長度)。然後,您正在更好地使用您創建的插值函數。

最後一個建議是嘗試a profiler。它可以讓你看到最需要花費的時間。

+0

您還可以使用np.ndenumerate減少兩個外部循環中的一個。這也應該稍微加快一點。 – Jose 2011-04-23 17:55:54