To pytanie i odpowiedź ma charakter kanoniczny (-ish) dotyczący dwuwymiarowej (i wielowymiarowej) interpolacji przy użyciu scipy. Często pojawiają się pytania dotyczące podstawowej składni różnych wielowymiarowych metod interpolacji, mam nadzieję, że również je wyjaśnię.

Mam zestaw rozproszonych dwuwymiarowych punktów danych i chciałbym wykreślić je jako ładną powierzchnię, najlepiej używając czegoś takiego jak contourflub plot_surfacew matplotlib.pyplot. Jak mogę interpolować moje dwuwymiarowe lub wielowymiarowe dane do siatki za pomocą Scipy?

Znalazłem scipy.interpolatepakiet podrzędny, ale nadal pojawiają się błędy podczas używania interp2dlub bisplreplub griddatalub rbf. Jaka jest prawidłowa składnia tych metod?

Zastrzeżenie: piszę ten post głównie z myślą o składni i ogólnym zachowaniu. Nie jestem zaznajomiony z aspektem pamięci i procesora opisanych metod i kieruję tę odpowiedź do tych, którzy mają dość małe zestawy danych, tak że jakość interpolacji może być głównym aspektem do rozważenia. Zdaję sobie sprawę, że pracując z bardzo dużymi zbiorami danych, skuteczniejsze metody (a mianowicie griddatai Rbf) mogą nie być wykonalne.

Porównam trzy rodzaje wielowymiarowych metod interpolacji ( interp2d/ splajny griddatai 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ą dwuwymiarową interpolację, ale realne metody można zastosować w dowolnych wymiarach. Każda metoda zapewnia różne rodzaje interpolacji; we wszystkich przypadkach użyję interpolacji sześciennej (lub czegoś zbliżonego do ¹ ). Ważne jest, aby pamiętać, że za każdym razem, gdy używasz interpolacji, wprowadzasz odchylenie w porównaniu z surowymi danymi, a określone metody mają wpływ na artefakty, które uzyskasz. Zawsze bądź tego świadomy i interpoluj odpowiedzialnie.

Będą dwa zadania interpolacji

1. upsampling (dane wejściowe na siatce prostokątnej, dane wyjściowe na siatce gęstszej)
2. interpolacja rozproszonych danych do zwykłej siatki

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

1. płynna i przyjazna funkcja cos(pi*x)*sin(pi*y):; zakres w[-1, 1]
2. funkcja zła (a zwłaszcza nieciągła): x*y/(x^2+y^2)o wartości 0,5 w pobliżu źródła; zakres w[-0.5, 0.5]

Oto jak wyglądają:

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

Dane testowe

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

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)

Płynna funkcja i upsampling

Zacznijmy od najłatwiejszego zadania. Oto jak działa upsampling z siatki kształtu [6,7]do jednej z [20,21]płynnych funkcji testu:

Mimo że jest to proste zadanie, istnieją już subtelne różnice między wynikami. Na pierwszy rzut oka wszystkie trzy wyniki są rozsądne. Warto zwrócić uwagę na dwie cechy, opierając się na naszej wcześniejszej wiedzy na temat podstawowej funkcji: środkowy przypadek griddatanajbardziej zniekształca dane. Zwróć uwagę na y==-1granicę wykresu (najbliżej xetykiety): funkcja powinna być ściśle zerowa (ponieważ y==-1jest linią węzłową dla funkcji gładkiej), ale tak nie jest griddata. Zwróć także uwagę na x==-1granicę wykresów (z tyłu, po lewej): podstawowa funkcja ma lokalne maksimum (implikujące zerowy gradient w pobliżu granicy) przy [-1, -0.5], ale griddatawynik pokazuje wyraźnie niezerowy gradient w tym regionie. Efekt jest subtelny, ale mimo wszystko jest to błąd. (WiernośćRbfjest jeszcze lepszy z domyślnym wyborem funkcji radialnych, dubbingowanych multiquadratic).

Zła funkcja i upsampling

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

Między tymi trzema metodami zaczynają się pojawiać wyraźne różnice. Patrząc na wykresy powierzchni, na wyjściu z interp2d(zwróć uwagę na dwa garby po prawej stronie wykreślonej powierzchni) widać wyraźne fałszywe ekstrema . Choć griddatai Rbfwydaje się produkować podobne wyniki na pierwszy rzut oka, ten ostatni wydaje się produkować głębsze minimum blisko [0.4, -0.4], że jest nieobecny z funkcją podstawową.

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

Gładka funkcja i rozproszone dane

Najczęściej chce się dokonać interpolacji na danych rozproszonych. Z tego powodu oczekuję, że te testy będą ważniejsze. Jak pokazano powyżej, punkty próbkowania wybrano pseudojednorodnie w interesującej nas dziedzinie. W realistycznych scenariuszach możesz mieć dodatkowy szum przy każdym pomiarze i powinieneś rozważyć, czy na początku ma sens interpolacja surowych danych.

Wyjście dla funkcji płynnej:

Teraz jest już trochę horroru. Przyciąłem dane wyjściowe od interp2ddo między [-1, 1]wyłącznie na potrzeby kreślenia, aby zachować przynajmniej minimalną ilość informacji. Oczywiste jest, że chociaż niektóre z podstawowych kształtów są obecne, istnieją ogromne hałaśliwe obszary, w których metoda całkowicie się psuje. Drugi przypadek griddatacałkiem ładnie odtwarza kształt, ale zwróć uwagę na białe obszary na granicy działki konturowej. Wynika to z faktu, że griddatadziała tylko wewnątrz wypukłego kadłuba punktów danych wejściowych (innymi słowy, nie dokonuje ż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, Rbfwydaje się działać najlepiej.

Zła funkcja i rozproszone dane

I moment, na który wszyscy czekaliśmy:

To żadna niespodzianka, że ​​się interp2dpoddaje. W rzeczywistości podczas rozmowy interp2dpowinieneś spodziewać się przyjaznych RuntimeWarningnarzekań na niemożność zbudowania splajnu. Jeśli chodzi o pozostałe dwie metody, Rbfwydaje się , że daje najlepsze wyniki, nawet w pobliżu granic dziedziny, w której wynik jest ekstrapolowany.

Pozwólcie więc, że powiem kilka słów na temat tych trzech metod, w kolejności malejącej preferencji (tak, aby najgorsze było najmniej prawdopodobne).

scipy.interpolate.Rbf

RbfKlasa oznacza „radialnych funkcji bazowych”. Szczerze mówiąc, nigdy nie rozważałem tego podejścia, dopóki nie zacząłem szukać informacji na temat tego postu, ale jestem prawie pewien, że będę ich używać w przyszłości.

Podobnie jak w przypadku metod opartych na splajnach (patrz dalej), użycie odbywa się w dwóch krokach: najpierw tworzy się wywoływalną Rbfinstancję klasy na podstawie danych wejściowych, a następnie wywołuje ten obiekt dla danej siatki wyjściowej, aby uzyskać interpolowany wynik. Przykład z testu płynnego 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 dane wyjściowe z_dense_smooth_rbfmają taki sam kształt jak x_densei y_densebez żadnego wysiłku. Należy również pamiętać, że Rbfobsługuje dowolne wymiary interpolacji.

Więc, scipy.interpolate.Rbf

• generuje grzeczne dane wyjściowe nawet dla szalonych danych wejściowych
• obsługuje interpolację w wyższych wymiarach
• ekstrapoluje poza wypukłą powłokę punktów wejściowych (oczywiście ekstrapolacja jest zawsze ryzykowna i generalnie nie należy w ogóle na niej polegać)
• tworzy interpolator jako pierwszy krok, więc ocena go w różnych punktach wyjściowych wymaga mniej dodatkowego wysiłku
• może mieć punkty wyjściowe o dowolnym kształcie (w przeciwieństwie do ograniczania do prostokątnych siatek, patrz dalej)
• skłonne do zachowania symetrii danych wejściowych
• obsługuje wiele rodzajów funkcji radialnych dla słowa kluczowego function: multiquadric, inverse, gaussian, linear, cubic, quintic, thin_platei zdefiniowany przez użytkownika dowolna

scipy.interpolate.griddata

Mój poprzedni faworyt griddatato ogólny konik roboczy do interpolacji w dowolnych wymiarach. Nie wykonuje ekstrapolacji poza ustaleniem jednej zadanej wartości dla punktów znajdujących się poza wypukłym kadłubem punktów węzłowych, ale ponieważ ekstrapolacja jest bardzo zmienną i niebezpieczną rzeczą, niekoniecznie jest to oszustwo. 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 niezdarną składnię. Punkty wejściowe muszą być określone w tablicy kształtu [N, D]w Dwymiarach. W tym celu musimy najpierw spłaszczyć nasze tablice współrzędnych 2D (używając ravel), a następnie połączyć tablice i transponować wynik. Można to zrobić na wiele sposobów, ale wszystkie wydają się być nieporęczne. Dane wejściowe zrównież muszą zostać spłaszczone. Mamy trochę więcej swobody, jeśli chodzi o punkty wyjściowe: z jakiegoś powodu można je również określić jako krotkę wielowymiarowych tablic. Zwróć uwagę, że helpof griddatajest mylący, 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

• generuje grzeczne dane wyjściowe nawet dla szalonych danych wejściowych
• obsługuje interpolację w wyższych wymiarach
• nie dokonuje ekstrapolacji, można ustawić pojedynczą wartość dla wyjścia poza wypukłą otoczką punktów wejściowych (patrz fill_value)
• oblicza interpolowane wartości 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 interpolację liniową w dowolnych wymiarach sześciennych w 1d i 2d. Stosuje się odpowiednio interpolację najbliższego sąsiada i interpolację liniową NearestNDInterpolatororaz LinearNDInterpolatorpod maską. Interpolacja sześcienna 1d wykorzystuje splajn, interpolacja sześcienna 2d wykorzystuje CloughTocher2DInterpolatordo konstruowania ciągłego różniczkowalnego interpolatora odcinkowo-sześciennego.
• może naruszyć symetrię danych wejściowych

scipy.interpolate.interp2d/scipy.interpolate.bisplrep

Jedynym powodem, dla którego omawiam interp2di jego krewnych, jest to, że ma zwodniczą nazwę i ludzie prawdopodobnie będą próbować go używać. Alert spoilera: nie używaj go (od wersji Scipy 0.17.0). Jest już bardziej wyjątkowy niż poprzednie przedmioty, ponieważ jest specjalnie używany do interpolacji dwuwymiarowej, ale podejrzewam, że jest to zdecydowanie najczęstszy przypadek interpolacji wielowymiarowej.

Jeśli chodzi o składnię, interp2djest podobna do Rbftego, że najpierw wymaga skonstruowania instancji interpolacji, którą można wywołać w celu podania rzeczywistych wartości interpolowanych. Jest jednak pewien haczyk: punkty wyjściowe muszą znajdować się na prostokątnej siatce, więc dane wejściowe do wywołania interpolatora muszą być wektorami 1d, które obejmują siatkę wyjściową, 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 interp2djest umieszczanie pełnych siatek 2d w wywołaniu interpolacyjnym, co prowadzi do gwałtownego zużycia pamięci i miejmy nadzieję, do pośpiechu MemoryError.

Największym problemem interp2djest to, że często nie działa. Aby to zrozumieć, musimy zajrzeć pod maskę. Okazuje się, że interp2djest to opakowanie dla funkcji niższego poziomu bisplrep+ bisplev, które z kolei są opakowaniami dla procedur FITPACK (napisanych w języku Fortran). Równoważne wywołanie z poprzedniego przykładu byłoby

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

A teraz sprawa interp2d: (w wersji scipy 0.17.0) jest miły komentarzinterpolate/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, bisplrepjest trochę konfiguracji 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. Podstawowe procedury interp2dnie są tak naprawdę przeznaczone do wykonywania interpolacji. Mogą wystarczyć do wystarczająco dobrze wychowanych danych, ale w realnych okolicznościach prawdopodobnie będziesz chciał użyć czegoś innego.

Podsumowując, interpolate.interp2d

• może prowadzić do artefaktów nawet w przypadku dobrze temperowanych danych
• jest przeznaczony specjalnie do problemów z dwiema zmiennymi (chociaż istnieje ograniczenie interpndla punktów wejściowych zdefiniowanych w siatce)
• dokonuje ekstrapolacji
• tworzy interpolator jako pierwszy krok, więc ocena go w różnych punktach wyjściowych wymaga mniej dodatkowego wysiłku
• może generować wyjście tylko przez prostokątną siatkę, dla wyjścia rozproszonego musiałbyś wywołać interpolator w pętli
• obsługuje interpolację liniową, sześcienną i kwintyczną
• może naruszyć symetrię danych wejściowych

¹ Jestem prawie pewien, że cubici linearrodzaj funkcji bazowych Rbfnie odpowiadają dokładnie innym interpolatorom o tej samej nazwie.
² Te NaN są również powodem, dla którego wykres powierzchni wydaje się tak dziwny: matplotlib historycznie miał trudności z wykreślaniem złożonych obiektów 3D z odpowiednimi informacjami o głębokości. Wartości NaN w danych dezorientują mechanizm renderujący, więc części powierzchni, które powinny znajdować się z tyłu, są wykreślane z przodu. Jest to problem dotyczący wizualizacji, a nie interpolacji.