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

я собираюсь сравнить три вида многомерных методов интерполяции (interp2d/сплайны, griddata и Rbf). Я подвергну их двум видам интерполяционных задач и двум видам базовых функций (точки, из которых должны быть интерполированы). Конкретные примеры продемонстрируют двумерную интерполяцию, но жизнеспособные методы применимы в произвольных измерениях. Каждый метод предоставляет различные виды интерполяции; во всех случаях я буду использовать кубическую интерполяцию (или что-то близкое¹). Важно отметить, что всякий раз, когда вы используете интерполяцию, вы вводите смещение по сравнению с вашими необработанными данными, и конкретные используемые методы влияют на артефакты, которые вы в конечном итоге получите. Всегда помните об этом, и интерполируйте ответственно.

две задачи интерполяции будет

1. upsampling (входные данные на прямоугольнике сетка, выходные данные находятся на более плотной сетке)
2. интерполяция рассеянных данных на регулярную сетку

две функции (над доменом [x,y] in [-1,1]x[-1,1]) будет

1. гладкая и дружественная функция:cos(pi*x)*sin(pi*y); в диапазоне [-1, 1]
2. злая (и в частности, не непрерывная) функция:x*y/(x^2+y^2) со значением 0,5 вблизи начала координат; диапазон в [-0.5, 0.5]

вот как они смотри:

сначала я продемонстрирую, как три метода ведут себя в этих четырех тестах, а затем я подробно расскажу о синтаксисе всех трех. Если вы знаете, что вы должны ожидать от метода, вы не можете тратить свое время на изучение его синтаксиса (глядя на вас, interp2d).

тестовых данных

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

import numpy as np
import scipy.interpolate as interp

# auxiliary function for mesh generation
def gimme_mesh(n):
    minval = -1
    maxval =  1
    # produce an asymmetric shape in order to catch issues with transpositions
    return np.meshgrid(np.linspace(minval,maxval,n), np.linspace(minval,maxval,n+1))

# set up underlying test functions, vectorized
def fun_smooth(x, y):
    return np.cos(np.pi*x)*np.sin(np.pi*y)

def fun_evil(x, y):
    # watch out for singular origin; function has no unique limit there
    return np.where(x**2+y**2>1e-10, x*y/(x**2+y**2), 0.5)

# sparse input mesh, 6x7 in shape
N_sparse = 6
x_sparse,y_sparse = gimme_mesh(N_sparse)
z_sparse_smooth = fun_smooth(x_sparse, y_sparse)
z_sparse_evil = fun_evil(x_sparse, y_sparse)

# scattered input points, 10^2 altogether (shape (100,))
N_scattered = 10
x_scattered,y_scattered = np.random.rand(2,N_scattered**2)*2 - 1
z_scattered_smooth = fun_smooth(x_scattered, y_scattered)
z_scattered_evil = fun_evil(x_scattered, y_scattered)

# dense output mesh, 20x21 in shape
N_dense = 20
x_dense,y_dense = gimme_mesh(N_dense)

плавная функция и апсэмплинг

давайте начнем с самой простой задачи. Вот как апсэмплинг из сетки формы [6,7] в одном из [20,21] работает вне для ровного теста функция:

несмотря на то, что это простая задача, уже есть тонкие различия между выходами. На первый взгляд все три выхода разумны. Есть две особенности, которые следует отметить, основываясь на наших предварительных знаниях о базовой функции: средний случай griddata искажает данные больше всего. Обратите внимание на y==-1 граница участка (ближайшая к x label): функция должна быть строго нулевой (так как y==-1 - это Узловая линия для гладкой функции), но это не так для griddata. Также обратите внимание на x==-1 граница участков (сзади, слева): базовая функция имеет локальный максимум (подразумевающий нулевой градиент вблизи границы) при [-1, -0.5], но griddata выход показывает явно ненулевой градиент в этой области. Эффект неуловимый, но тем не менее это предубеждение. (Верность Rbf еще лучше с выбором по умолчанию радиальных функций, дублированных multiquadratic.)

злые функции и апсэмплинг

немного сложнее задача выполнить апсэмплинг на нашей злой функции:

четкие различия начинают проявляться среди трех методов. Глядя на поверхностные графики, есть четкие паразитные экстремумы, появляющиеся на выходе из interp2d (обратите внимание на два горба на правой стороне нанесенной поверхности). В то время как griddata и Rbf производят похожие результаты на первый взгляд, последний, кажется, производит более глубокий минимум вблизи [0.4, -0.4] это отсутствует в основной функции.

однако есть один важный аспект, в котором Rbf значительно превосходит: он уважает симметрию основной функции (что, конечно, также стало возможным благодаря симметрии сетки образца). Выход из griddata нарушает симметрию точек выборки, которая уже слабо видна в гладком случае.

ровное функция и рассеянные данные

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

выход для гладкой функция:

теперь там уже немного Шоу ужасов происходит. Я обрезал выход из interp2d до [-1, 1] исключительно для построения графиков, чтобы сохранить хотя бы минимальный объем информации. Ясно, что в то время как некоторые из основных форм присутствуют, есть огромные шумные области, где метод полностью ломается. Второй случай griddata воспроизводит форму довольно красиво, но обратите внимание на белый области на границе контурного участка. Это связано с тем, что griddata работает только внутри выпуклой оболочки исходных точек данных (другими словами, он не выполняет каких-либо экстраполяция). Я сохранил значение NaN по умолчанию для выходных точек, лежащих вне выпуклой оболочки.² учитывая эти особенности, Rbf кажется, работает лучше всего.

злые функции и рассеянные данные

и момент, которого мы все ждали ибо:

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

---

Итак, позвольте мне сказать несколько слов о трех методах, в убывающий порядок предпочтений (так что худшее, скорее всего, будет прочитано кем-либо).

scipy.interpolate.Rbf

The Rbf класс означает "радиальные базисные функции". Честно говоря, я никогда не рассматривал этот подход, пока не начал исследовать этот пост, но я уверен, что буду использовать их в будущем.

так же, как и методы на основе сплайнов (см. ниже), использование происходит в два этапа: первый создает вызываемый Rbf экземпляр класса на входных данных, а затем вызывает этот объект для данного вывода сетки для получения интерполированного результата. Пример из теста smooth upsampling:

import scipy.interpolate as interp
zfun_smooth_rbf = interp.Rbf(x_sparse, y_sparse, z_sparse_smooth, function='cubic', smooth=0)  # default smooth=0 for interpolation
z_dense_smooth_rbf = zfun_smooth_rbf(x_dense, y_dense)  # not really a function, but a callable class instance

обратите внимание, что как входные, так и выходные точки были 2d-массивами в этом случае, а выход z_dense_smooth_rbf имеет ту же форму, что и x_dense и y_dense без каких-либо усилий. Также обратите внимание, что Rbf поддерживает произвольные размеры для интерполяции.

и scipy.interpolate.Rbf

• производит хорошее поведение выход даже для сумасшедших входных данных
• поддерживает интерполяцию в более высоких измерениях
• экстраполирует вне выпуклой оболочки входных точек (конечно, экстраполяция всегда азартная игра, и вы вообще не должны на нее полагаться)
• создает интерполятор в качестве первого шага, поэтому оценка его в различных выходных точках меньше дополнительных усилий
• может иметь выходные точки произвольной формы (в отличие от ограничения прямоугольные сетки, см. Позже)
• склонны к сохранению симметрии входных данных
• поддерживает несколько видов радиальных функций для сайта function:multiquadric,inverse,gaussian,linear,cubic,quintic,thin_plate и определяемый пользователем произвольный

scipy.interpolate.griddata

мой бывший любимый, griddata, является общей рабочей лошадкой для интерполяции в произвольных измерениях. Он не выполняет экстраполяцию за пределы настройки одного заданного значения для точек вне выпуклой оболочки узловых точек, но так как экстраполяция-это очень непостоянная и опасная вещь, это не обязательно обман. Пример использования:

z_dense_smooth_griddata = interp.griddata(np.array([x_sparse.ravel(),y_sparse.ravel()]).T,
                                          z_sparse_smooth.ravel(),
                                          (x_dense,y_dense), method='cubic')   # default method is linear

обратите внимание на слегка kludgy синтаксис. Входные точки должны быть указаны в массиве формы [N, D] in D габариты. Для этого мы сначала должны сгладить наши 2d массивы координат (используя ravel), затем объединить массивы и транспонировать результат. Существует несколько способов чтобы сделать это, но все они кажутся громоздкими. Вход z данные также должны быть сглажены. У нас есть немного больше свободы, когда дело доходит до выходных точек: по какой-то причине они также могут быть указаны как кортеж многомерных массивов. Обратите внимание, что help на griddata вводит в заблуждение, так как предполагает, что то же самое верно для input точки (по крайней мере для версии 0.17.0):

griddata(points, values, xi, method='linear', fill_value=nan, rescale=False)
    Interpolate unstructured D-dimensional data.

    Parameters
    ----------
    points : ndarray of floats, shape (n, D)
        Data point coordinates. Can either be an array of
        shape (n, D), or a tuple of `ndim` arrays.
    values : ndarray of float or complex, shape (n,)
        Data values.
    xi : ndarray of float, shape (M, D)
        Points at which to interpolate data.

в двух словах scipy.interpolate.griddata

• производит хорошо себя ведет вывод даже для сумасшедших входных данных
• поддерживает интерполяцию в более высоких измерениях
• не выполнять экстраполяцию, одно значение может быть установлено для выхода за пределы выпуклой оболочки исходных точек (см. fill_value)
• вычисляет интерполированные значения в одном вызове, поэтому зондирование нескольких наборов выходных точек начинается с нуля
• может иметь выходные точки произвольной формы
• поддерживает ближайшего соседа и линейная интерполяция в произвольной размеры, кубические в 1D и 2D. Ближайшего соседа и линейной интерполяции использовать NearestNDInterpolator и LinearNDInterpolator под капотом, соответственно. 1D кубическая интерполяция использует сплайн, 2d кубическая интерполяция использует CloughTocher2DInterpolator построить непрерывно дифференцируемый кусочно-кубический интерполятор.
• может нарушить симметрию входных данных

scipy.interpolate.interp2d/scipy.interpolate.bisplrep

единственная причина, по которой я обсуждаем interp2d а его родственники считают, что у него обманчивое название, и люди, скорее всего, попытаются его использовать. Предупреждение о спойлере: не используйте его (начиная с версии scipy 0.17.0). Он уже более особенный, чем предыдущие темы, поскольку он специально используется для двумерной интерполяции, но я подозреваю, что это, безусловно, самый распространенный случай для многомерной интерполяции.

что касается синтаксиса,interp2d похож на Rbf в том, что сначала нужно построить интерполяцию экземпляр, который может быть вызван для предоставления фактических интерполированных значений. Однако есть уловка: выходные точки должны быть расположены на прямоугольной сетке, поэтому входы, входящие в вызов интерполятора, должны быть 1D векторами, которые охватывают выходную сетку, как будто из numpy.meshgrid:

# reminder: x_sparse and y_sparse are of shape [6, 7] from numpy.meshgrid
zfun_smooth_interp2d = interp.interp2d(x_sparse, y_sparse, z_sparse_smooth, kind='cubic')   # default kind is 'linear'
# reminder: x_dense and y_dense are of shape [20, 21] from numpy.meshgrid
xvec = x_dense[0,:] # 1d array of unique x values, 20 elements
yvec = y_dense[:,0] # 1d array of unique y values, 21 elements
z_dense_smooth_interp2d = zfun_smooth_interp2d(xvec,yvec)   # output is [20, 21]-shaped array

одна из самых распространенных ошибок при использовании interp2d помещает ваши полные 2d сетки в вызов интерполяции, что приводит к взрывному потреблению памяти и, надеюсь, к поспешности MemoryError.

теперь, самая большая проблема с interp2d это то, что он часто не работает. Чтобы понять это, надо заглянуть под капот. Оказывается, что interp2d это обертка для функций нижнего уровня bisplrep+bisplev, которые в свою очередь являются обертками для подпрограмм FITPACK (написанных на Fortran). Эквивалентный вызов в предыдущем примере будет

kind = 'cubic'
if kind=='linear':
    kx=ky=1
elif kind=='cubic':
    kx=ky=3
elif kind=='quintic':
    kx=ky=5
# bisplrep constructs a spline representation, bisplev evaluates the spline at given points
bisp_smooth = interp.bisplrep(x_sparse.ravel(),y_sparse.ravel(),z_sparse_smooth.ravel(),kx=kx,ky=ky,s=0)
z_dense_smooth_bisplrep = interp.bisplev(xvec,yvec,bisp_smooth).T  # note the transpose

теперь, вот в чем дело interp2d: (в scipy версии 0.17.0) есть хороший комментарий interpolate/interpolate.py на interp2d:

if not rectangular_grid:
    # TODO: surfit is really not meant for interpolation!
    self.tck = fitpack.bisplrep(x, y, z, kx=kx, ky=ky, s=0.0)

и вообще в interpolate/fitpack.py, в bisplrep есть некоторые настройки и в конечном итоге

tx, ty, c, o = _fitpack._surfit(x, y, z, w, xb, xe, yb, ye, kx, ky,
                                task, s, eps, tx, ty, nxest, nyest,
                                wrk, lwrk1, lwrk2)                 

и это все. Процедуры, лежащие в основе interp2d на самом деле не предназначены для выполнения интерполяции. Их может быть достаточно для достаточно корректных данных, но при реальных обстоятельствах вы, вероятно, захотите использовать что-то еще.

просто в заключение, interpolate.interp2d

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

---

¹я совершенно уверен, что cubic и linear вид базисных функций Rbf точно не соответствуют другим интерполяторам с тем же именем.
²эти NaNs также являются причиной того, почему поверхностный участок кажется таким странным: matplotlib исторически имеет трудности с построением сложных 3d-объектов с надлежащей информацией о глубине. Значения NaN в данных путают рендерер, поэтому части поверхности, которые должны быть сзади, отображаются спереди. Это проблема визуализации, а не интерполяции.