Interpolación ciencia de matriz grande.

Tengo un ndarray (Z) con unos 500000 elementos en una cuadrícula rectangular (X, Y).

Ahora quiero interpolar valores en unas 100 ubicaciones en x, y que no están necesariamente en la cuadrícula.

Tengo un código trabajando en Matlab:

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

Sin embargo, cuando trato de usar el mismo enfoque con scipy.interpolate, obtengo varios errores según el método. Por ejemplo, interp2d falla con MemoryError si especifico kind = 'linear' y “OverflowError: hay demasiados puntos de datos para interpolar” si especifico kind='cubic' . También probé Rbf y bisplev pero también fallan.

Entonces, la pregunta es: ¿existe una función de interpolación que permita interpolaciones de matrices grandes? ¿Hay otra solución al problema? (¿O tengo que codificar una función que elija el área adecuada alrededor de los puntos para interpolar y luego las llamadas interp2d?)

Además: ¿Cómo hacer esto con números complejos?

Como sus datos están en una cuadrícula, puede usar RectBivariateSpline .

Para trabajar con números complejos, puede interpolar data.real y data.imag separado (las rutinas FITPACK IIRC no manejan datos complejos).

editar: Whoops. ¡Apenas me di cuenta de que el OP sugirió esta solución en la pregunta!

No sé por qué las rutinas de interpolación tardan tanto tiempo y memoria en encontrar los nudos de datos estructurados, pero como solo está utilizando una pequeña parte de la cuadrícula completa, podría dividir su interpolación en parches para hacer las cosas más eficientes. .

 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])