Interpolacja

Do interpolacji użyję funkcji interp1d (interpolacja jednowymiarowa) z pakietu scipy. Żeby móc z niej korzystać trzeba ją „zaimportować” (włączyć do projektu z biblioteki).

Funkcja interp1d jako argumenty przyjmuje dwa dektory danych x i y. Można jej również określić rodzaj interpolacji. Może to być: 'linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic' lub stopień wielomianu interpolacyjnego.

Biblioteka numpy to wygodna biblioteka do operacji na tablicach z danymi.

matplotlib.pyplot to biblioteka oferująca podstawowe funkcje rysowania wykresów (dosyć podobna do tej z Matlaba).

In [1]:
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (12,9) # Tu definiujemy rozmiary wykresu. Standardowo wydaje się być zbyt mały.
In [2]:
import numpy as np
from scipy.interpolate import interp1d
from scipy.interpolate import lagrange
from scipy.interpolate import BarycentricInterpolator

Nasze done to wartości losowe z zakresu od 0 do 1

In [3]:
N = 21 # liczba losowanych wartości
M = 5 * N # liczba wartości do wyliczenia funkcji interpolowanej
p = 3 # stopień wielomianu interpolacyjnego
In [4]:
x0 = np.linspace(-1, 1, N) # Oś x zawiera się między -1 a 1
y0 = np.random.random(N)

Przygotowujemy wykres

In [5]:
plt.plot(x0, y0, 'o', label='Data');
In [6]:
x = np.linspace(-1, 1, M)

tablica options zawiara możliwe typy interpoacji. W pętli będziemy je przebiegać.

In [7]:
options = ('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic', p)
In [8]:
for o in options:
    f = interp1d(x0, y0, kind=o)    # interpolacja
    plt.plot(x, f(x), label=o)      # wykres funkcji interpolującej
In [9]:
plt.legend() # opis wykresu
plt.show()   # generacja wykresu
No handles with labels found to put in legend.

Jak widać w każdym z tych przypadków funkcja interpolująca przechodzi przez wygenerowane punkty. między nimi, zaś, może zachowywać się dosyć dziwnie. I odczytywanie danych pomiędzy węzłami interpolacji traktowane jest — na ogół — jako czary albo wróżenie z fusów.

Teraz dane z pomiarów

In [10]:
data = np.loadtxt("AVERAGE300.dat")
In [11]:
plt.plot(data[:,0],data[:,1],label='.')
Out[11]:
[<matplotlib.lines.Line2D at 0x7f7e4fc02eb8>]
In [12]:
plt.show()

Fltracja danych

W tym wypadku filtracja będzie oznaczać odrzucenie danych złych.

Jak widać dane zawierają błędy. Prawdopodobnie wszystkie wartości mniejsze ok -25 są błędne. Trzeba je usunąć. Funkcja min podaje minimalną wartość z tablicy. Funkcja np.argmin numer indeksu na którym ta wartość się realizuje.

In [13]:
print('min =', min(data[:,1]))
print('argmin =', np.argmin(data[:,1]))
min = -11577.819208
argmin = 178

Pojedyńczą wartość można usunąć za pomocą funkcji np.delete, ale nie jest to najwygodniejsza metoda gdy złych wartości jest więcej.

In [14]:
temperatura = np.delete(data, (178), axis=0);
min(temperatura[:,1])
print(data.shape[0]) # to będzie liczba wierszy
print(temperatura.shape[0])
288
287

Masked array

Chyba lepszym pomysłem będzie stworzenie macierzy „z maską” (masked array). Warto nauczyć się tej techniki.

In [15]:
import numpy.ma as ma 
In [16]:
temperatura = ma.masked_less(data[:,1], -25)

Z pierwszej kolumny wybieramy dane które są mniejsze niż -25 (stopni celsjusza). Powstaje obiekt (temperatura) składający się z danych oraz tablicy logicznej zawierających prawdę wszędzie tam gdzie warunke jest spełniony.

ma.getmask(temperatura) (albo temperatura.mask) zwraca tablicę logiczną.

temperatura.data (albo ma.getdata(temperatura) zwraca „dobre” wartości.

In [17]:
ma.getmask(temperatura)
Out[17]:
array([False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False])
In [18]:
np.shape(temperatura[~temperatura.mask].data)[0]
Out[18]:
263

Spróbujmy obejrzeć „zamaskowane” wartości:

In [19]:
temperatura[temperatura.mask].data
Out[19]:
array([-11577.819208 ,  -1069.5577776,  -1069.5577776,  -1069.5577776,
        -1069.5577776,  -1069.5577776,  -1069.5577776,  -1069.5577776,
        -1069.5577776,  -1069.5577776,  -1069.5577776,  -1069.5577776,
        -1069.5577776,  -1069.5577776,  -1069.5577776,  -1069.5577776,
        -1069.5577776,  -1069.5577776,  -1069.5577776,  -1069.5577776,
        -1069.5577776,  -1069.5577776,  -1069.5577776,  -1069.5577776,
        -1069.5577776])

Teraz musimy wybrać współrzędne czasowe poprawnych wartości z pierwszej kolumny macierzy data. Stworzymy kolejną tablicę z maską (używając maski uzyskanej z temperatury). Funkcja ma.array służy do tworzenia tablicy z maską:

In [20]:
czas=ma.array(data[:,0], mask = temperatura.mask)

Skoro dostęp do zamaskowanych wartości daje temperatura.mask to dostęp do warości poprawnych da negacja, czyli ~temparatura.mask.

In [21]:
CZAS = czas[~czas.mask].data # tylko dane
CZAS = CZAS - min(CZAS) # dzięki temu współrzędne x są liczone względem początku doby
TEMPERATURA = temperatura[~temperatura.mask].data # tylko dane
In [22]:
plt.plot(CZAS[0:30],TEMPERATURA[0:30],label='.')
Out[22]:
[<matplotlib.lines.Line2D at 0x7f7e4fb72fd0>]
In [23]:
f = lagrange(CZAS[0:30], TEMPERATURA[0:30])
x = np.linspace(0,9000, 200)
plt.plot(x, f(x))
Out[23]:
[<matplotlib.lines.Line2D at 0x7f7e4fb52940>]
In [24]:
plt.show()
In [25]:
CZAS[0:30]
Out[25]:
array([   0.,  300.,  600.,  900., 1200., 1500., 1800., 2100., 2400.,
       2700., 3000., 3300., 3600., 3900., 4200., 4500., 4800., 5100.,
       5400., 5700., 6000., 6300., 6600., 6900., 7200., 7500., 7800.,
       8100., 8400., 8700.])

Jak widać ten rodzaj interpolacji, w tym przypadku nie sprawdza się (choć wyniki dla 10 węzłów wyglądają interesująco)

In [26]:
f = BarycentricInterpolator(CZAS[0:30], TEMPERATURA[0:30])
x = np.linspace(1000,7800, 200)
plt.plot(CZAS[5:30],TEMPERATURA[5:30],label='.')
plt.plot(x, f(x))
plt.show()

Teraz funkcja trygonometryczna (sin)

In [27]:
from math import pi
from math import sqrt
from math import sin
In [28]:
X=np.array([0,  pi/6,   pi/4,   pi/3,   pi/2,   2*pi/3, 3*pi/4, 5*pi/6, 
            pi, 7*pi/6, 5*pi/4, 4*pi/3, 3*pi/2, 5*pi/3, 7*pi/4, 11*pi/6, 
            2*pi])
In [29]:
Y=np.array([0,  1/2,  sqrt(2)/2,  sqrt(3)/2,   1,  sqrt(3)/2,  sqrt(2)/2,  1/2, 
            0, -1/2, -sqrt(2)/2, -sqrt(3)/2,  -1, -sqrt(3)/2, -sqrt(2)/2, -1/2, 
            0])
In [30]:
my_sin=BarycentricInterpolator(X,Y)
In [31]:
x=np.linspace(0, 2*pi, 100)
In [32]:
plt.plot(X,Y)
plt.plot(x,my_sin(x))
plt.show()
In [33]:
s = np.sin(x)
In [34]:
plt.plot(x, s-my_sin(x))
Out[34]:
[<matplotlib.lines.Line2D at 0x7f7e4f9ec9b0>]
In [35]:
plt.show()

Zwracam uwagę, że interpolaca jest bardzo dobra! „Przeregulowania” na początku i końcu przedziału co do wartości bezwzględnej są mniejsze niż $1^{-9}$.

In [ ]: