In [2]:
import pandas as pd
import numpy as np
import scipy.io as sio
In [3]:
from pykalman import KalmanFilter
from scipy.signal import lfilter, cheby1, cheby2, butter
In [4]:
import matplotlib.pyplot as plt
import seaborn as sns
from os import path
import glob
%matplotlib inline
sns.set('notebook')
current_palette = sns.color_palette()
In [5]:
grimm = pd.read_csv('data/alignedGrimmData.csv')
grimm.index = grimm['Unix Time'].values/3600.0
del grimm['Unix Time']
In [6]:
grimm.index
Out[6]:
In [7]:
plt.plot(grimm.index, grimm['PM 2.5'])
Out[7]:
In [8]:
def load_specks(filename, cov=1075):
speckdata = sio.loadmat(filename)
specknames = ['s1']#,'s2','s3']
specks = []
for s in specknames:
speck = pd.DataFrame(speckdata[s], columns=['time','humidity','conc','count','raw','temp'])
speck.index = pd.to_datetime(speck['time'], unit='s')
speck.time = speck.time/3600.0
kf = KalmanFilter(initial_state_mean=0, n_dim_obs=1, transition_matrices=[[1]], observation_matrices=[[1]], observation_covariance=[[cov]], transition_covariance=[[1]])
speck['kfiltered'] = kf.filter(speck.raw.values)[0]
b3, a3 = cheby1(2, 0.5, 0.06)
b4, a4 = butter(1, [0.03, 0.95], btype='bandstop')
speck['filtered'] = lfilter(b4,a4, lfilter(b3, a3, speck.kfiltered))
specks.append(speck)
return specks
In [9]:
specks = load_specks('data/alignedGrimmSpeckData.mat')
In [10]:
import plotly.plotly as pypp
import plotly.graph_objs as pygo
In [14]:
cov = 1075
kf = KalmanFilter(initial_state_mean=0, n_dim_obs=1, transition_matrices=[[1]], observation_matrices=[[1]], observation_covariance=[[cov]], transition_covariance=[[1]])
speck = specks[0]
speck['ksmooth'] = kf.smooth(speck.raw.values)[0]
In [15]:
subs=20
speckf = pygo.Scatter(x=speck.time[::subs], y=speck.ksmooth[::subs]/10, opacity=0.5, name='Filtered')
speckc = pygo.Scatter(x=speck.time[::subs], y=speck.conc[::subs], opacity=0.5, name='Concentration')
grimmd = pygo.Scatter(x=grimm.index[::10], y=grimm['PM 2.5'][::10], opacity=0.5, name='GRIMM2.5')
pypp.iplot([grimmd, speckc, speckf])
Out[15]:
In [16]:
plt.plot(grimm.index/3600.0, grimm['PM 2.5'], alpha=0.5)
for speck in specks:
plt.plot(speck.time/3600.0, speck.conc, alpha=0.5)
plt.plot(speck.time/3600.0, speck.filtered/10.0, alpha=0.5)
In [17]:
grimm.columns
Out[17]:
In [29]:
particle_volumes = {x:(float(x)**3)/6 for x in grimm.columns[4:]}
In [53]:
grimm.columns[4:]
Out[53]:
In [52]:
print particle_volumes
In [45]:
grimm_volumes['0.8']/=2
In [46]:
grimm_volumes = grimm.copy()
for x in particle_volumes:
grimm_volumes[x] *= particle_volumes[x]
In [47]:
grimm_volumes.head()
Out[47]:
In [31]:
grimm.head()
Out[31]:
In [48]:
import itertools
def group(n, iterable):
args = [iter(iterable)] * n
return ([e for e in t if e != None] for t in itertools.izip_longest(*args))
In [51]:
for columns in group(3, grimm.columns):
grimm[columns].plot(alpha=0.5)
In [50]:
for columns in group(3, grimm.columns):
grimm_volumes[columns].plot(alpha=0.5)
In [19]:
specks[0].humidity.plot()
Out[19]:
In [ ]: