SciPy интерполяция большой матрицы

У меня есть ndarray (Z) с примерно 500000 элементов на прямоугольной сетке (X, Y).

Теперь я хочу интерполировать значения примерно в 100 местах по x, y, которые не обязательно находятся в сетке.

У меня есть код, работающий в Matlab:

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

Однако, когда я пытаюсь использовать тот же подход с scipy.interpolate, я получаю различные ошибки в зависимости от метода. Например, interp2d завершается ошибкой MemoryError, если я указываю kind = 'linear', и "OverflowError: слишком много точек данных для интерполяции", если я указываю kind='cubic'. Я также пробовал Rbf и bisplev, но они тоже не работают.

Итак, вопрос: существует ли функция интерполяции, которая позволяет интерполировать большие матрицы? Есть ли другое решение проблемы? (Или мне нужно написать функцию, которая выбирает подходящую область вокруг точек для интерполяции и затем вызывает interp2d?)

Кроме того: Как это сделать с комплексными числами?


person Joma    schedule 16.03.2011    source источник
comment
Хотите показать свой код? 500000 не так уж и много. Спасибо   -  person eat    schedule 16.03.2011


Ответы (2)


Поскольку ваши данные находятся в сетке, вы можете использовать RectBivariateSpline.

Чтобы работать с комплексными числами, вы можете интерполировать data.real и data.imag отдельно (процедуры FITPACK IIRC не обрабатывают сложные данные).

person pv.    schedule 16.03.2011

изменить: Упс. Только что понял, что ОП предложил это решение в вопросе!

Я не знаю, почему процедурам интерполяции требуется так много времени и памяти, чтобы найти узлы структурированных данных, но, поскольку вы используете только небольшие части полной сетки, вы можете разбить интерполяцию на патчи, чтобы сделать работу более эффективной. .

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])
person Paul    schedule 16.03.2011