Zastrzeżenie: piszę ten post głównie z rozważań składniowych i ogólnego zachowania w umyśle. Nie jestem zaznajomiony z aspektem pamięci i procesora opisywanych metod i kieruję tę odpowiedź do tych, którzy mają dość małe zbiory danych, tak że jakość interpolacji może być głównym aspektem do rozważenia. Zdaję sobie sprawę, że podczas pracy z bardzo dużymi zbiorami danych, lepiej działające metody (mianowicie griddata i Rbf) mogą nie być wykonalne.

Zamierzam porównaj trzy rodzaje metod interpolacji wielowymiarowej(interp2d/spliny, griddata oraz Rbf). Poddam je dwóm rodzajom zadań interpolacyjnych i dwóm rodzajom funkcji bazowych (punkty, z których mają być interpolowane). Konkretne przykłady zademonstrują interpolację dwuwymiarową, ale realne metody mają zastosowanie w dowolnych wymiarach. Każda metoda zapewnia różne rodzaje interpolacji; we wszystkich przypadkach użyję sześciennych interpolacja (lub coś podobnego¹). Ważne jest, aby pamiętać, że za każdym razem, gdy korzystasz z interpolacji, wprowadzasz błąd w porównaniu z surowymi danymi, a konkretne metody wpływają na Artefakty, z którymi skończysz. Zawsze bądź tego świadomy i Interpoluj odpowiedzialnie.

Dwa zadania interpolacji będą

Upsampling (dane wejściowe znajdują się na prostokątnej siatce, dane wyjściowe na gęstszej siatce)]}
1. interpolacja rozproszonych danych na siatka regularna

Dwie funkcje (nad domeną [x,y] in [-1,1]x[-1,1]) będą

1. funkcja gładka i przyjazna: cos(pi*x)*sin(pi*y); zakres w [-1, 1]
2. funkcja zła (a w szczególności nieciągła): x*y/(x^2+y^2) o wartości 0,5 blisko początku; zakres w [-0.5, 0.5]

Oto jak wyglądają:

Najpierw zademonstruję, jak te trzy metody zachowują się w tych czterech testach, a następnie szczegółowo przedstawię składnię wszystkich trzy. Jeśli wiesz, czego powinieneś oczekiwać od metody, możesz nie chcieć tracić czasu na naukę jej składni (patrząc na Ciebie, interp2d).

Dane testowe

Dla jasności, oto kod, za pomocą którego wygenerowałem dane wejściowe. Podczas gdy w tym konkretnym przypadku jestem oczywiście świadomy funkcji leżącej u podstaw danych, będę używać tego tylko do generowania danych wejściowych dla metod interpolacji. Używam numpy dla wygody (i głównie do generowania danych) , ale scipy sama też by wystarczyła.

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)

Funkcja gładka i upsampling

Zacznijmy od najprostszego zadania. Oto jak upsampling z siatki kształtu [6,7] do jednej z [20,21] Działa dla funkcji gładkiego testu:

Choć jest to proste zadanie, istnieją już subtelne różnice między wyjściami. Na pierwszy rzut oka wszystkie trzy wyjścia są rozsądne. Istnieją dwie cechy, które należy zauważyć, w oparciu o naszą wcześniejszą wiedzę na temat funkcja podstawowa: środkowy przypadek griddata najbardziej zniekształca dane. Zwróć uwagę na y==-1 granicę wykresu (najbliższą etykiecie x): funkcja powinna być ściśle zerowa (ponieważ y==-1 jest linią węzłową dla funkcji gładkiej), jednak tak nie jest dla griddata. Zwróć również uwagę na x==-1 granice Wykresów (z tyłu, po lewej): funkcja bazowa ma lokalne maksimum (sugerujące zerowy gradient w pobliżu granicy) w [-1, -0.5], jednak wyjście griddata pokazuje wyraźnie niezerowy gradient w tym regionie. Efekt jest subtelny, ale mimo to jest to stronniczość. (Wierność Rbf jest jeszcze lepsza przy domyślnym wyborze funkcji promieniowych, nazwanych multiquadratic.)

Zła funkcja i upsampling

Nieco trudniejszym zadaniem jest wykonanie upsamplingu na naszej złej funkcji:]}

Widoczne są wyraźne różnice pomiędzy tymi trzema metodami. Patrząc na wykresy powierzchniowe, w wyjściu z interp2d pojawiają się wyraźne skrajności fałszywe (zwróć uwagę na dwa Garby po prawej stronie wykreślonej powierzchni). Podczas gdy griddata i Rbf wydają się dawać podobne wyniki na pierwszy rzut oka, te ostatnie wydają się produkować głębsze minimum w pobliżu [0.4, -0.4], które jest nieobecne w podstawowej funkcji.

Jest jednak jeden istotny aspekt, w którym Rbf jest o wiele lepszy: szanuje symetrię funkcji bazowej (co jest oczywiście możliwe dzięki symetrii siatki próbki). Wyjście z griddata łamie symetrię punktów próbki, co jest już słabo widoczne w gładkiej obudowie.

Funkcja gładka i dane rozproszone

NajczÄ ™ Ĺ " ciej chce siÄ ™ wykonaÄ ‡ interpolacjä ™ na danych rozproszonych. Z tego powodu oczekuję, że te testy będą ważniejsze. Jak pokazano powyżej, punkty próbne zostały wybrane pseudo jednolicie w dziedzinie zainteresowania. W realistycznych scenariuszach przy każdym pomiarze może wystąpić dodatkowy szum i należy rozważyć, czy warto interpolować surowe dane, aby rozpocząć z.

Wyjście dla funkcji gładkiej:

Teraz jest już trochę horroru dzieje. Przyciąłem wyjście z interp2d do between [-1, 1] wyłącznie do kreślenia, aby zachować przynajmniej minimalną ilość informacji. Jest jasne, że podczas gdy niektóre z podstawowych kształt jest obecny, istnieją ogromne hałaśliwe regiony, w których metoda całkowicie rozpada się. Drugi przypadek griddata odtwarza kształt dość ładnie, ale zwróć uwagę na białe regiony na granicy działki konturowej. Wynika to z faktu, że griddata działa tylko wewnątrz wypukłego kadłuba wejściowych punktów danych (innymi słowy, nie wykonuje żadnej ekstrapolacji). Zachowałem domyślną wartość NaN dla punktów wyjściowych leżących poza wypukłym kadłubem.² biorąc pod uwagę te cechy, Rbf wydaje się działać najlepiej.

Zła funkcja i rozproszone dane

And the moment we ' ve all been waiting za:

Nie jest wielką niespodzianką, że interp2d poddaje się. W rzeczywistości, podczas wywołania do interp2d należy spodziewać się przyjaznych RuntimeWarnings narzekających na niemożność skonstruowania splajnu. Jeśli chodzi o pozostałe dwie metody, {[9] } wydaje się produkować najlepsze wyniki, nawet w pobliżu granic domeny, w której wynik jest ekstrapolowany.

---

Więc pozwól mi powiedzieć kilka słów o trzech metodach, w kolejności malejącej preferencji (tak, że najgorsza jest najmniej prawdopodobna do przeczytania przez kogokolwiek).

scipy.interpolate.Rbf

Klasa Rbf oznacza "radial basis functions". Szczerze mówiąc, nigdy nie brałem pod uwagę tego podejścia, dopóki nie zacząłem badać tego postu, ale jestem prawie pewien, że będę z nich korzystać w przyszłości.

Podobnie jak metody oparte na spline( patrz później), użycie odbywa się w dwóch krokach: pierwszy tworzy wywoływalną instancję klasy Rbf na podstawie danych wejściowych, a następnie wywołuje ten obiekt dla danego wyjścia siatka w celu uzyskania interpolowanego wyniku. Przykład z testu gładkiego upsamplingu:

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

Zauważ, że zarówno punkty wejściowe, jak i wyjściowe były w tym przypadku tablicami 2d, A wyjście z_dense_smooth_rbf ma taki sam kształt jak x_dense i y_dense bez żadnego wysiłku. Zauważ również, że Rbf obsługuje dowolne wymiary dla interpolacji.

Więc, scipy.interpolate.Rbf

• tworzy dobrze zachowujące się wyjście nawet dla szalonych danych wejściowych
• wspiera interpolację w wyższych wymiary
[113]}ekstrapolacja poza wypukłym kadłubem punktów wejściowych (oczywiście ekstrapolacja jest zawsze hazardem i generalnie nie należy na niej w ogóle polegać) {114]} [113]}tworzy interpolator jako pierwszy krok, więc ocena go w różnych punktach wyjściowych jest mniejszym dodatkowym wysiłkiem Może mieć punkty wyjściowe o dowolnym kształcie (w przeciwieństwie do ograniczonych do prostokątnych siatek, patrz później).]}
• podatne na zachowanie symetrii wejścia dane
• obsługuje wiele rodzajów funkcji promieniowych dla słów kluczowych function: multiquadric, inverse, gaussian, linear, cubic, quintic, thin_plate i definiowane przez użytkownika dowolne

scipy.interpolate.griddata

Mój poprzedni ulubieniec, griddata, jest ogólnym koniem roboczym do interpolacji w dowolnych wymiarach. Nie wykonuje ekstrapolacji poza ustawieniem pojedynczej wartości domyślnej dla punktów poza wypukłym kadłubem punktów węzłowych, ale ponieważ ekstrapolacja jest bardzo zmienna i niebezpieczna, jest to nie koniecznie przekręt. przykład użycia:

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

Zwróć uwagę na nieco kludyczną składnię. Punkty wejściowe muszą być określone w tablicy kształtu [N, D] w wymiarach D. W tym celu najpierw musimy spłaszczyć nasze tablice współrzędnych 2d (używając ravel), a następnie połączyć tablice i transponować wynik. Istnieje wiele sposobów, aby to zrobić, ale wszystkie z nich wydają się być nieporęczne. Dane wejściowe z również muszą zostać spłaszczone. Mamy nieco większą swobodę jeśli chodzi o punkty wyjściowe: dla niektórych powód mogą być również określone jako krotka wielowymiarowych tablic. Należy zauważyć, że help Z griddata jest mylące, ponieważ sugeruje, że to samo dotyczy punktów wejściowych (przynajmniej dla wersji 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.

W skrócie, scipy.interpolate.griddata

• tworzy dobrze zachowujące się wyjście nawet dla szalonych danych wejściowych
• obsługuje interpolację w wyższych wymiarach
• nie wykonuje ekstrapolacji, można ustawić pojedynczą wartość dla wyjścia poza wypukły kadłub punktów wejściowych (patrz fill_value)
• oblicza wartości interpolowane w jednym wywołaniu, więc sondowanie wielu zestawów punktów wyjściowych zaczyna się od zera
• może mieć punkty wyjściowe o dowolnym kształcie
• obsługuje interpolację najbliższego sąsiada i liniową w dowolnych wymiarach, sześciennych w 1d i 2D. Interpolacja najbliższego sąsiada i liniowa używa odpowiednio NearestNDInterpolator i LinearNDInterpolator pod maską. Interpolacja sześcienna 1D wykorzystuje splajn, interpolacja sześcienna 2D używa CloughTocher2DInterpolator do skonstruowania interpolatora różniczkowalnego w sposób ciągły.
• może naruszać symetrię danych wejściowych

scipy.interpolate.interp2d/scipy.interpolate.bisplrep

Jedynym powodem, dla którego omawiam interp2di jego krewnych jest to, że ma zwodnicze imię, a ludzie prawdopodobnie spróbują go użyć. Uwaga Spoiler: nie używaj go (od wersji scipy 0.17.0). Jest już bardziej wyjątkowy niż poprzednie przedmioty, ponieważ jest specjalnie używany do dwuwymiarowych interpolacji, ale podejrzewam, że jest to zdecydowanie najczęstszy przypadek interpolacji wielowymiarowej.

Jeśli chodzi o składnię, {[10] } jest podobna do Rbf w tym, że najpierw trzeba zbudować instancję interpolacji, którą można wywołać, aby dostarczyć rzeczywiste wartości interpolowane. Jest jednak pewien haczyk: punkty wyjściowe muszą być umieszczone na prostokątnej siatce, więc wejścia do wywołania interpolatora muszą być wektorami 1d, które rozciągają się na siatkę wyjściową, tak jakby z 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

Jednym z najczęstszych błędów podczas używania interp2d jest umieszczenie pełnych siatek 2d w wywołaniu interpolacyjnym, co prowadzi do wybuchowego zużycia pamięci i miejmy nadzieję, że pochopnego MemoryError.

Największym problemem z

Jest to, że często nie działa. Aby to zrozumieć, musimy zajrzeć pod maskę. Okazuje się, że interp2d jest opakowaniem dla funkcji niższego poziomubisplrep+bisplev, które z kolei są owijarkami dla procedur FITPACK (zapisanych w Fortran). Analogiczne wywołanie do poprzedniego przykładu to

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

Oto rzecz o interp2d: (w wersji 0.17.0) jest fajny komentarz w interpolate/interpolate.py dla 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)

I rzeczywiście w interpolate/fitpack.py, w bisplrep jest jakaś konfiguracja i ostatecznie

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)                 

I to wszystko. Procedury leżące u podstaw interp2d nie są tak naprawdę przeznaczone do wykonywania interpolacji. Mogą wystarczyć do dostatecznie dobrze zachowanych danych, ale w realistycznych okolicznościach prawdopodobnie będziesz chciał użyć czegoś innego.

Tak na zakończenie, interpolate.interp2d

Może prowadzić do artefaktów nawet przy dobrze hartowanych danych.]} Jest to rozwiązanie przeznaczone dla problemów dwuargumentowych (choć istnieje ograniczenie interpn dla punktów wejściowych zdefiniowanych na siatce)
• wykonuje ekstrapolację
[113]}tworzy interpolator jako pierwszy krok, więc ocena go w różnych punktach wyjściowych jest mniejszym dodatkowym wysiłkiem
• może tylko aby uzyskać wyjście na prostokątną siatkę, trzeba będzie wywołać interpolator w pętli
• obsługuje interpolację liniową, sześcienną i kwintyczną
• może naruszać symetrię danych wejściowych

---

¹jestem całkiem pewien, że cubic i linear funkcje bazowe Rbf nie odpowiadają dokładnie innym interpolatorom o tej samej nazwie.
²te Nany są również powodem, dla którego powierzchnia fabuła wydaje się taka dziwna: matplotlib historycznie ma trudności z kreśleniem złożonych obiektów 3d z odpowiednią informacją o głębi. Wartości NaN w danych mylą renderer, więc części powierzchni, które powinny znajdować się z tyłu, są wykreślane tak, aby znajdowały się z przodu. Jest to problem z wizualizacją, a nie interpolacją.