In [127]:

    

%pylab inline
import pandas as pd
from pykalman import KalmanFilter

df = pd.read_csv("../data/ChungCheonDC/CompositeETCdata.csv")
df_DC = pd.read_csv("../data/ChungCheonDC/CompositeDCdata.csv")
df_DCprc = pd.read_csv("../data/ChungCheonDC/CompositeDCdata_processed.csv")
df_DCstd = pd.read_csv("../data/ChungCheonDC/CompositeDCstddata.csv")


    

    




Populating the interactive namespace from numpy and matplotlib

In [128]:

    

# missininds = np.arange(df_DC[electrodeID[elecind]].values.size)[np.isnan(df_DC[electrodeID[elecind]].values)]
electrodeID = df_DC.keys()[1:-1]


    

In [129]:

    

from scipy import interpolate
sys.path.append("../codes/")
from DCdata import readReservoirDC_all
directory = "../data/ChungCheonDC/"
dat_temp,height_temp, ID = readReservoirDC_all(directory+"20151231180000.apr")
locs = dat_temp[:,:4]
mida = locs[:,:2].sum(axis=1)
midb = locs[:,2:].sum(axis=1)
mid = (mida + midb)*0.5
dz = mida-midb
x = np.linspace(mid.min(), mid.max(), 100)
z = np.linspace(dz.min(), dz.max(), 100)
grid_x, grid_z = np.meshgrid(x,z)

def vizDCtimeSeries(idatum, itime, itime_ref, colors, flag, df_DC):
    fig = plt.figure(figsize = (12, 12))
    ax1 = plt.subplot(411)
    ax2 = plt.subplot(412)
    
    valsratio = df_DC[electrodeID].values[itime,:].flatten() / df_DC[electrodeID].values[itime_ref,:].flatten()
    valsDC = np.log10(df_DC[electrodeID].values[itime,:].flatten())
    valsDCstd = df_DCstd[electrodeID].values[itime,:].flatten()
    grid_rho_ratio = griddata(mid, dz, valsratio, grid_x, grid_z, interp='linear')
    grid_rho_ratio = grid_rho_ratio.reshape(grid_x.shape)
    if flag =="std":
        vmin, vmax = 0, 10
        grid_rho = griddata(mid, dz, valsDCstd, grid_x, grid_z, interp='linear')        
    elif flag =="rho":
        vmin, vmax = np.log10(20), np.log10(200)
        grid_rho = griddata(mid, dz, valsDC, grid_x, grid_z, interp='linear')
    grid_rho = grid_rho.reshape(grid_x.shape)
        
    
    ax1.contourf(grid_x, grid_z, grid_rho, 200, vmin =vmin, vmax = vmax, clim=(vmin, vmax), cmap="jet")    
    vmin, vmax = 0.9, 1.1
    ax2.contourf(grid_x, grid_z, grid_rho_ratio, 200, vmin =vmin, vmax = vmax, clim=(vmin, vmax), cmap="jet")        
    ax1.scatter(mid, dz, s=20, c = valsDC, edgecolor="None", vmin =vmin, vmax = vmax, clim=(vmin, vmax))
    ax1.plot(mid, dz, 'k.')
    ax2.scatter(mid, dz, s=20, c = valsratio, edgecolor="None", vmin =vmin, vmax = vmax, clim=(vmin, vmax))
    ax2.plot(mid, dz, 'k.')
    
    for i in range(len(colors)):
        ax1.plot(mid[idatum[i]], dz[idatum[i]], 'o', color=colors[i])    
        ax2.plot(mid[idatum[i]], dz[idatum[i]], 'o', color=colors[i])    
        

    ax3 = plt.subplot(413)
    ax3_1 = ax3.twinx()
    df.plot(x='date', y='reservoirH', ax=ax3_1, color='k', linestyle='-', lw=2)
    df.plot(x='date', y='upperH_med', ax=ax3_1, color='b', linestyle='-', lw=2)
    df.plot(x='date', y='Temp (degree)', ax=ax3, color='r', linestyle='-', lw=2)
    df.plot(x='date', y='Rainfall (mm)', ax=ax3, color='b', linestyle='-', marker="o", lw=2)   # ms=4)
    ax3.legend(loc=3, bbox_to_anchor=(1.05, 0.7))
    ax3_1.legend(loc=3, bbox_to_anchor=(1.05, 0.4))
    itime_ref0 = itime_ref
    itime_ref1 = itime
    ax3.plot(np.r_[itime_ref0, itime_ref0], np.r_[-5, 40], 'k--', lw=2)
    ax3.plot(np.r_[itime_ref1, itime_ref1], np.r_[-5, 40], 'k--', lw=2)

    ax4 = plt.subplot(414)
    ax4_1 = ax4.twinx()
    ax4.legend(loc=3, bbox_to_anchor=(1.05, 0.7))
    ax4.set_yscale('log')
    temp = df_DC[electrodeID[elecind]].values
    vmax = np.median(temp[~np.isnan(temp)]) + np.std(temp[~np.isnan(temp)])*20
    vmin = np.median(temp[~np.isnan(temp)]) - np.std(temp[~np.isnan(temp)])*20
    ax4.plot(np.r_[itime_ref1, itime_ref1], np.r_[vmin, vmax], 'k--', lw=2)
    ax4.plot(np.r_[itime_ref0, itime_ref0], np.r_[vmin, vmax], 'k--', lw=2)
    
    df_DC.plot(x='date', y=electrodeID[idatum], ax=ax4,color='r')
    df.plot(x='date', y='Rainfall (mm)', ax=ax4_1, color='b', linestyle='-', marker="o", ms=4)
  
    ax4.set_ylim(vmin, vmax)


    

In [130]:

    

ax1 = plt.subplot(111)
ax1_1 = ax1.twinx()
df.plot(figsize=(12,3), x='date', y='reservoirH', ax=ax1_1, color='k', linestyle='-', lw=2)
df.plot(figsize=(12,3), x='date', y='upperH_med', ax=ax1_1, color='b', linestyle='-', lw=2)
df.plot(figsize=(12,3), x='date', y='Temp (degree)', ax=ax1, color='r', linestyle='-', lw=2)
ax1.legend(loc=3, bbox_to_anchor=(1.05, 0.7))
ax1_1.legend(loc=3, bbox_to_anchor=(1.05, 0.4))
itime_ref0 = 255
itime_ref1 = 115
ax1.plot(np.r_[itime_ref0, itime_ref0], np.r_[-5, 35], 'k-')
ax1.plot(np.r_[itime_ref1, itime_ref1], np.r_[-5, 35], 'k-')
# print df['date'].values[itime_ref]


    

    
Out[130]:





[<matplotlib.lines.Line2D at 0x11892a210>]






    

In [131]:

    

# ax1 = plt.subplot(111)
# ax1_1 = ax1.twinx()
# df_DC.plot(figsize=(12,3), x='date', y=electrodeID[elecind], ax=ax1, colors=['k', 'b', 'r'])
# df.plot(figsize=(12,3), x='date', y='reservoirH', ax=ax1_1, color='k', linestyle='-', lw=2)
# ax1.legend(loc=3, bbox_to_anchor=(1.05, 0.7))
# ax1_1.legend(loc=3, bbox_to_anchor=(1.05, 0.4))
# ax1.set_yscale('linear')


    

In [132]:

    

# ax1 = plt.subplot(111)
# df_DCstd.plot(figsize=(12,3), x='date', y=electrodeID[elecind], ax=ax1, colors=['k', 'b', 'r'], linestyle="-", marker='.', lw=1)
# ax1.set_yscale('log')
# ax1.legend(loc=3, bbox_to_anchor=(1.05, 0.7))


    

In [133]:

    

txrxID =  df_DC.keys()[1:-1]
xmasking = lambda x: np.ma.masked_where(np.isnan(x.values), x.values)


    

In [134]:

    

#x= electrodeID[elecind] 
x= df_DC[txrxID]
max3 = pd.rolling_max(x, 3)


    

    




/Users/sklim/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:3: FutureWarning: pd.rolling_max is deprecated for DataFrame and will be removed in a future version, replace with 
	DataFrame.rolling(window=3,center=False).max()
  app.launch_new_instance()

In [135]:

    

from ipywidgets import interact


    

In [136]:

    

# making matrix like max3 (but with zeros)
newdata = np.zeros_like(max3)


    

In [137]:

    

newdata.shape


    

    
Out[137]:





(365, 380)

In [138]:

    

df_DC_new = df_DC.copy()  
for i,index in enumerate(txrxID):
    df_DC_new.loc[:,index] = newdata[:,i].flatten()
# df_DC_new.to_csv("../data/ChungCheonDC/CompositeDCdata_processed.csv")


    

In [146]:

    

from ipywidgets import interact, IntSlider, ToggleButtons
itime = 201
itime_ref = 201
print df['date'].values[itime]
elecind = [5, 150,200]
 

# vizDCtimeSeries(elecind, itime, itime_ref, ['k','b','r'])
viz = lambda idatum, itime, flag: vizDCtimeSeries([idatum], itime, itime_ref, ['r'], flag, df_DC_new)
interact(viz, idatum=IntSlider(min=0, max=379, step=1, value=144)\
         ,itime=IntSlider(min=0, max=360, step=1, value=201)\
         ,flag=ToggleButtons(options=["std", "rho"]))


    

    




---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-146-87dace07c2d1> in <lambda>(idatum, itime, flag)
      7 
      8 # vizDCtimeSeries(elecind, itime, itime_ref, ['k','b','r'])
----> 9 viz = lambda idatum, itime, flag: vizDCtimeSeries([idatum], itime, itime_ref, ['r'], flag, df_DC_new)
     10 interact(viz, idatum=IntSlider(min=0, max=379, step=1, value=144)         ,itime=IntSlider(min=0, max=360, step=1, value=201)         ,flag=ToggleButtons(options=["std", "rho"]))

<ipython-input-129-e5cc0812882e> in vizDCtimeSeries(idatum, itime, itime_ref, colors, flag, df_DC)
     65     vmax = np.median(temp[~np.isnan(temp)]) + np.std(temp[~np.isnan(temp)])*20
     66     vmin = np.median(temp[~np.isnan(temp)]) - np.std(temp[~np.isnan(temp)])*20
---> 67     ax4.plot(np.r_[itime_ref1, itime_ref1], np.r_[vmin, vmax], 'k--', lw=2)
     68     ax4.plot(np.r_[itime_ref0, itime_ref0], np.r_[vmin, vmax], 'k--', lw=2)
     69 

/Users/sklim/anaconda/lib/python2.7/site-packages/matplotlib/__init__.pyc in inner(ax, *args, **kwargs)
   1817                     warnings.warn(msg % (label_namer, func.__name__),
   1818                                   RuntimeWarning, stacklevel=2)
-> 1819             return func(ax, *args, **kwargs)
   1820         pre_doc = inner.__doc__
   1821         if pre_doc is None:

/Users/sklim/anaconda/lib/python2.7/site-packages/matplotlib/axes/_axes.pyc in plot(self, *args, **kwargs)
   1384             lines.append(line)
   1385 
-> 1386         self.autoscale_view(scalex=scalex, scaley=scaley)
   1387         return lines
   1388 

/Users/sklim/anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.pyc in autoscale_view(self, tight, scalex, scaley)
   2190                 y1 += delta
   2191             if not _tight:
-> 2192                 y0, y1 = ylocator.view_limits(y0, y1)
   2193             self.set_ybound(y0, y1)
   2194 

/Users/sklim/anaconda/lib/python2.7/site-packages/matplotlib/ticker.pyc in view_limits(self, vmin, vmax)
   1615         if minpos <= 0 or not np.isfinite(minpos):
   1616             raise ValueError(
-> 1617                 "Data has no positive values, and therefore can not be "
   1618                 "log-scaled.")
   1619 

ValueError: Data has no positive values, and therefore can not be log-scaled.





    
Out[146]:





<function __main__.<lambda>>






    

In [140]:

    

print df['date'].values[itime_ref]


    

    




2015-07-21

In [141]:

    

for i in range(0,379,100):
    x= df_DC[txrxID[i]]
    x1 = df_DC_new[txrxID[i]]
    plt.plot(newdata[:,i], 'k')
    plt.plot(x1, 'ro')    
    plt.plot(x, 'k.', ms=2)


    

In [142]:

    

i = 245
x= df_DC[txrxID[i]]
#median10 = pd.rolling_median(x, 6)
mean10 = pd.rolling_max(x, 3)
#x1 = median10
#x2 = mean10
# Masking array having NaN
xm = xmasking(mean10)
# Construct a Kalman filter
kf = KalmanFilter(transition_matrices = [1],
                  observation_matrices = [1],
                  initial_state_mean = 67.6,
                  initial_state_covariance = 1,
                  observation_covariance=1,
                  transition_covariance=1)
# Use the observed values of the price to get a rolling mean
state_means, _ = kf.filter(xm)

#plt.plot(x1)
plt.plot(x)
#plt.plot(x1)
#plt.plot(x2)
plt.plot(state_means)
plt.legend([  'origin x','Kalman'])


    

    




/Users/sklim/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:4: FutureWarning: pd.rolling_max is deprecated for Series and will be removed in a future version, replace with 
	Series.rolling(window=3,center=False).max()






    
Out[142]:





<matplotlib.legend.Legend at 0x112eb6e50>






    

In [143]:

    

y='Rainfall (mm)'


    

In [144]:

    

print y


    

    




Rainfall (mm)

In [145]:

    

i = 144


x= df_DC[txrxID[i]]
#median10 = pd.rolling_median(x, 6)
mean10 = pd.rolling_max(x, 3)
#x1 = median10
#x2 = mean10
# Masking array having NaN
xm = xmasking(mean10)
# Construct a Kalman filter
kf = KalmanFilter(transition_matrices = [1],
                  observation_matrices = [1],
                  initial_state_mean = 67.6,
                  initial_state_covariance = 1,
                  observation_covariance=1,
                  transition_covariance=1)
# Use the observed values of the price to get a rolling mean
state_means, _ = kf.filter(xm)

#plt.plot(x1)
plt.plot(x)
#plt.plot(x1)
#plt.plot(x2)
plt.plot(state_means)
plt.legend([  'origin x','Kalman'])


    

    




/Users/sklim/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:6: FutureWarning: pd.rolling_max is deprecated for Series and will be removed in a future version, replace with 
	Series.rolling(window=3,center=False).max()






    
Out[145]:





<matplotlib.legend.Legend at 0x11dd76950>






    

In [ ]:

    




    

In [ ]: