2011-03-16 83 views
2

我有一個在矩形網格(X,Y)上有大約500000個元素的ndarray(Z)。大矩陣的SciPy插值

現在我想在x,y中的大約100個位置插值,這些位置不一定在網格上。

我有一些代碼在Matlab工作:

data = interp2(X,Y,Z, x,y); 

然而,當我嘗試使用與scipy.interpolate同樣的方法,我得到取決於方法的各種錯誤。例如,如果我指定kind = 'linear'和「溢出錯誤:太多的數據點進行插值」,如果我指定了kind='cubic',則interp2d將失敗並顯示MemoryError。我也試過Rbfbisplev,但他們也失敗了。

所以問題是:是否有插值函數允許插值大矩陣?有沒有解決問題的另一種方法? (或者我必須編碼一個函數,選擇點周圍的適當區域來插入並調用interp2d?)

另外:如何用複數來做到這一點?

+0

小心顯示您的代碼? 500000不是那麼大。謝謝 – eat 2011-03-16 16:18:32

回答

2

由於您的數據位於網格上,因此您可以使用RectBivariateSpline

要處理複雜數字,可以分別插入data.realdata.imag(FITPACK例程IIRC不處理複雜數據)。

1

編輯:哎呦。剛剛意識到OP在問題中提出了這個解決方案!

我不知道爲什麼插值例程花費很多時間和內存來查找結構化數據的節點,但由於您只使用整個網格的一小部分,因此可以將插值分解爲多個補丁讓事情更有效率。

from scipy import interpolate 
import numpy as np 

def my_interp(X, Y, Z, x, y, spn=3): 
    xs,ys = map(np.array,(x,y)) 
    z = np.zeros(xs.shape) 
    for i,(x,y) in enumerate(zip(xs,ys)): 
     # get the indices of the nearest x,y 
     xi = np.argmin(np.abs(X[0,:]-x)) 
     yi = np.argmin(np.abs(Y[:,0]-y)) 
     xlo = max(xi-spn, 0) 
     ylo = max(yi-spn, 0) 
     xhi = min(xi+spn, X[0,:].size) 
     yhi = min(yi+spn, Y[:,0].size) 
     # make slices of X,Y,Z that are only a few items wide 
     nX = X[xlo:xhi, ylo:yhi] 
     nY = Y[xlo:xhi, ylo:yhi] 
     nZ = Z[xlo:xhi, ylo:yhi] 
     intp = interpolate.interp2d(nX, nY, nZ) 
     z[i] = intp(x,y)[0] 
    return z 

N = 1000 
X,Y = np.meshgrid(np.arange(N), np.arange(N)) 
Z = np.random.random((N, N)) 

print my_interp(X, Y, Z, [13.2, 999.9], [0.01, 45.3]) 
+0

感謝您的示例代碼! – Joma 2011-04-05 15:18:01