#!/usr/bin/env python
# coding: utf-8

# # Czytanie danych zawierających ciągi znaków NaN
# 
# NaN to oznaczenie ,,Not a Number'' które może pojawić się (na wydrukach z programu komputerowego) gdy procesor wykona jakąś operację ,,niedozwoloną''.
# 
# Dzielenie przez zero wykonywane jest przez współczesne bez problemów, generując jako wynik wartość ,,nieskończoność'' (która podczas wyprowdzania oznaczona będzie `Inf` albo `-Inf`). Procesor poprawnie wykonuje proste operacje na nieskończonościach, na przykłąd `Inf + 1` to będzie `Inf`. Natomiast próba dodawania $-\infty + \infty$ wygeneruje wynik `NaN`.
# ### Dygresja: przykład NaNów i nieskończoności

# In[1]:


duze = 1e300
male = 1e-300
dzielenie = duze/male


# In[2]:


dzielenie


# In[3]:


dzielenie +5


# In[4]:


dzielenie - dzielenie


# ### Koniec dygresji
# Dane generowane przez RRDB będą zawierały takie znaczniki wszędzie tam, gdzie z jakichś powodów, nie udało się dokonać pomiarów. Na przykład plik AVERAGE86400.dat zawiera takie wartości 

# In[5]:


import numpy as np
import numpy.ma as ma


# In[6]:


import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (12,9) # Tu definiujemy rozmiary wykresu. 
# Standardowo wydaje się być zbyt mały.


# Funkcja [`genfromtext()`](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.genfromtxt.html#numpy.genfromtxt) pozwala takie dane importować (ale zachowując wartości `nan` w wybranych punktach. Sprawdźmy co się dzieje w 19 wierszu pliku `AVERAGE86400.dat`?

# In[7]:


np.genfromtxt('AVERAGE86400.dat')[19,:]


# In[8]:


np.isnan(np.genfromtxt('AVERAGE86400.dat')[19,1])


# In[9]:


np.isnan(np.genfromtxt('AVERAGE86400.dat')[19,0])


# Funkcja [`np.isnan()`](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.isnan.html) pozwala sprawdzić czy zmienna przyjmuje wartośc NaN.
# 
# Zatem, wydaje się, że można tu użyć metody maskowania tablic opisanej [gdzie indziej]() używając funkcji [`masked_invalid`](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.ma.masked_invalid.html#numpy.ma.masked_invalid) lub [`masked_where`](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.ma.masked_where.html#numpy.ma.masked_where).
# ### Odrzucenie NaNów

# In[10]:


dane = np.genfromtxt('AVERAGE1800.dat')


# In[11]:


temperatura = ma.masked_invalid(dane)


# In[12]:


plt.plot(dane[:,0],dane[:,1],'.');
plt.show()


# ### Odrzucenie wartości błędnych
# Okazuje się, że dane ciągle są nieprawdopodobne (zawierają błędne wartości). Zatem trzeba je filtrować ze względu na te wartości i NaN-y. Przyjmijmy, że temperatury większe od $55^o$ i mniejsze od $-25^o$ są nieprawdopodobne.
# 
# Można to zapewne zrobić na wiele sposobów: albo stworzyć dwie maski i połączyć je operatorem OR, albo spróbować od razu stworzyć maskę na podstawie złożonego warunku. Skorzystamy z funkcji `np.logical_or`, która pozwala wykonywać operacje logiczne na elemantach całych tablic.

# In[13]:


temp = dane[:,1]
temp = ma.masked_where(np.logical_or(np.isnan(temp), np.logical_or(temp < -25, temp > 55)), temp)
czas = ma.masked_array(dane[:,0], mask=temp.mask)


# In[14]:


temp[~temp.mask].data


# In[15]:


czas = czas[~czas.mask].data
temp = temp[~temp.mask].data


# In[16]:


plt.plot(czas, temp, '.')
plt.show()


# Teraz powinniśmy jeszcze zastąpić czymś odrzucone wartości. Można, w tym celu, posłużyć sę interpolacją jak to opisano w innym notatniku. Nie trzeba — choć można — robić resamplingu, wystarczy użyć interpolacji aby zastąpić „puste miejsca” na wykresie danymi z interpolacji.
# 
# ### Opis osi X
# 
# Unix timestamp może jest i wygodny, ale trudno zorientować się jakiego okresu czasu dotyczy powyższy wykres. Chcielibyśmy mieś os X opisaną jako daty.
# 
# W tym celu należy (niestety) skonwertować timestamp do postaci zrozumiałej dla oprogramowania. Służy do tego funkcja [`fromtimestamp`](https://docs.python.org/3/library/datetime.html#datetime-objects)

# In[17]:


import matplotlib.dates as md
import datetime as dt


# In[18]:


dates=[dt.datetime.fromtimestamp(ts) for ts in czas]


# In[19]:


plt.xticks( rotation=25 )
ax = plt.gca()
xfmt = md.DateFormatter('%Y-%m-%d')
ax.xaxis.set_major_formatter(xfmt)
plt.plot(dates, temp, '.')
plt.show()


# Czasami będziemy potrzebowali aby zastąpić czymś wartości w odrzuconych punktach. Można wówczas posłużyć się interpolacją w sposób opisany [gdzie indziej]().