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).
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (12,9) # Tu definiujemy rozmiary wykresu. Standardowo wydaje się być zbyt mały.
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
N = 21 # liczba losowanych wartości
M = 5 * N # liczba wartości do wyliczenia funkcji interpolowanej
p = 3 # stopień wielomianu interpolacyjnego
x0 = np.linspace(-1, 1, N) # Oś x zawiera się między -1 a 1
y0 = np.random.random(N)
Przygotowujemy wykres
plt.plot(x0, y0, 'o', label='Data');
x = np.linspace(-1, 1, M)
tablica options
zawiara możliwe typy interpoacji. W pętli będziemy je przebiegać.
options = ('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic', p)
for o in options:
f = interp1d(x0, y0, kind=o) # interpolacja
plt.plot(x, f(x), label=o) # wykres funkcji interpolującej
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.
data = np.loadtxt("AVERAGE300.dat")
plt.plot(data[:,0],data[:,1],label='.')
[<matplotlib.lines.Line2D at 0x7efe275822b0>]
plt.show()
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.
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.
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
Chyba lepszym pomysłem będzie stworzenie macierzy „z maską” (masked array). Warto nauczyć się tej techniki.
import numpy.ma as ma
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.
ma.getmask(temperatura)
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])
np.shape(temperatura[~temperatura.mask].data)[0]
263
Spróbujmy obejrzeć „zamaskowane” wartości:
temperatura[temperatura.mask].data
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ą:
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
.
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
plt.plot(CZAS[0:30],TEMPERATURA[0:30],label='.')
[<matplotlib.lines.Line2D at 0x7efe2755d8b0>]
f = lagrange(CZAS[0:30], TEMPERATURA[0:30])
x = np.linspace(0,9000, 200)
plt.plot(x, f(x))
[<matplotlib.lines.Line2D at 0x7efe274b0be0>]
plt.show()
CZAS[0:30]
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)
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()
from math import pi
from math import sqrt
from math import sin
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])
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])
my_sin=BarycentricInterpolator(X,Y)
x=np.linspace(0, 2*pi, 100)
plt.plot(X,Y)
plt.plot(x,my_sin(x))
plt.show()
s = np.sin(x)
plt.plot(x, s-my_sin(x))
[<matplotlib.lines.Line2D at 0x7efe259f2940>]
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}$.