SciPy интерполяция большой матрицы - PullRequest
Новая
       6

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

2 голосов
/ 16 марта 2011

У меня есть 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?)

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

Ответы [ 2 ]

2 голосов
/ 16 марта 2011

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

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

2 голосов
/ 16 марта 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])
...