%matplotlib inline
import pandas as pd
from math import pi
from bokeh.models import ColumnDataSource
from bokeh.models import HoverTool
from bokeh.plotting import figure, output_file, show
from bokeh.models import LinearAxis, Range1d
from bokeh.layouts import gridplot
from scipy.signal import detrend
print("# pandas version = ",pd.__version__)
import bokeh as bk
print("# bokeh version = ",bk.__version__)
import scipy as sp
print("# scipy version = ",sp.__version__)
# ftp directory
pfo_dir = "http://ncedc.org/ftp/outgoing/taira/PFO"
# sts1v data
pfo_sts1v_fi = pfo_dir + "/pfo_sts1v_0.00001-0.01hz.300sps.acc.ts.wh.out"
print(pfo_sts1v_fi)
# pfo volumetric strain data. assuming poisson's ratio = 0.25
pfo_evolume_fi = pfo_dir + "/ez_volume_nu0.25_.ts.wh.out"
print(pfo_evolume_fi)
# pfo gravitatity data
pfo_grav_fi = pfo_dir + "/PFO.grav.out.ts.wh.out"
# CIB data
cib_fi = pfo_dir + "/CIB.WTHT.ts.wh.out"
# CIC data
cic_fi = pfo_dir + "/CIC.WTHT.ts.wh.out"
# read file
# format
# time acc (nano m/s*2)
#2009-7-14T0:0:00 -53298
pfo_sts1v =pd.read_csv(pfo_sts1v_fi,
sep=" ",names=["time", "acc_nano"],
header=None, skiprows=1)
#print(hhz_psd_ts)
# convert "time" to datetime. need for plot
pfo_sts1v['time'] = pd.to_datetime(pfo_sts1v['time'])
# check
pfo_sts1v.head()
# check
pfo_sts1v.tail()
# PFO strain data
# time vol ez0 (az=0deg) ez90 (az=90deg)
pfo_evolume =pd.read_csv(pfo_evolume_fi,
sep=",",names=["time", "vol","ez0","ez90"],
header=None, skiprows=1)
pfo_evolume.head()
# time format conversion
pfo_evolume['time'] = pd.to_datetime(pfo_evolume['time'])
# PFO gravity data
# time gravity microgal, positive for decreasing
# microgal = microgal2. same data
pfo_grav =pd.read_csv(pfo_grav_fi,
sep=" ",names=["time", "microgal", "microgal2"],
header=None, skiprows=1)
# microgal2 is Nan.... Not sure why but microgal2 is not used
pfo_grav.head()
# time format conversion
pfo_grav['time'] = pd.to_datetime(pfo_grav['time'])
# check
pfo_grav['time'].head()
# CIB data
# time counts meters
cib =pd.read_csv(cib_fi,
sep=" ",names=["time", "counts","meters"],
header=None, skiprows=1)
# check
cib.head()
# time format conversion
cib['time'] = pd.to_datetime(cib['time'])
# check
cib['time'].head()
# trim data. same as the sts1 data
cib_select = cib[ ("2009-07-14 00:00:00" <= cib['time']) & (cib['time'] <= "2009-09-21 23:55:00") ]
# check
cib_select.head()
# CIC data
# time counts meters
cic =pd.read_csv(cic_fi,
sep=" ",names=["time", "counts","meters"],
header=None, skiprows=1)
# check
cic.head()
# time format conversion
cic['time'] = pd.to_datetime(cic['time'])
# check
cic['time'].head()
# trim data
cic_select = cic[ ("2009-07-14 00:00:00" <= cic['time']) & (cic['time'] <= "2009-09-21 23:55:00") ]
# check
cic_select.head()
# trim data
pfo_grav_select = pfo_grav[ ("2009-07-14 00:00:00" <= pfo_grav['time']) & (pfo_grav['time'] <= "2009-09-21 23:55:00") ]
# check
pfo_grav_select.head()
# trim daat
pfo_evolume_select = pfo_evolume[ ("2009-07-14 00:00:00" <= pfo_evolume['time']) & (pfo_evolume['time'] <= "2009-09-21 23:55:00") ]
# check
pfo_evolume_select.head()
# dataframe to numpy matrix
evolume_matrix = pfo_evolume_select.as_matrix()
# check to see evol
evolume_matrix[:,1] # 0 time, 1 evol # 2 ez0 #3 ez90
# now add evol data as pfo_sts1v['vol']
pfo_sts1v['vol'] =evolume_matrix[:,1]
# check
pfo_sts1v['vol'].head()
# check
pfo_sts1v['time'].head()
# dataframe to numpy matrix
grav_matrix = pfo_grav_select.as_matrix()
# check
grav_matrix[:,1] # microgal
# now add microgal data as pfo_sts1v['microgal']
pfo_sts1v['microgal'] = grav_matrix[:,1]
# check
pfo_sts1v['microgal'].head()
# dataframe to numpy matrix
cib_matrix = cib_select.as_matrix()
# check
cib_matrix[:,2] # 0 time, 1 counts, 2 meters
# now add cib data as pfo_sts1v['cib']
pfo_sts1v['cib'] = detrend(cib_matrix[:,2])
# dataframe to numpy matrix
cic_matrix = cic_select.as_matrix()
# check
cic_matrix[:,2] # 0 time, 1 counts, 2 meters
# now add cib data as pfo_sts1v['cic']
pfo_sts1v['cic'] = detrend(cic_matrix[:,2])
#TOOLS = "pan,wheel_zoom,box_zoom,undo,reset,save"
TOOLS = "pan,box_zoom,undo,reset,save"
#data = {'x_values': [1, 2, 3, 4, 5],
# 'y_values': [6, 7, 2, 3, 6]}
# make "data". all parameters need to be added
data = {
'time': pfo_sts1v['time'],
'acc_mili': pfo_sts1v['acc_nano']/1000/1000,
'microgal': pfo_sts1v['microgal'],
'vol': pfo_sts1v['vol'],
'cib': pfo_sts1v['cib'],
'cic': pfo_sts1v['cic'],
#'time': pfo_evolume['time'],
}
# make source. needed
source = ColumnDataSource(data=data)
# x_axis is datetime. add tools. define size and add title
p = figure(x_axis_type="datetime", tools=TOOLS, plot_width=600, height=600, title = "PFO")
# labels
p.xaxis.axis_label = 'Time (UTC)'
p.yaxis.axis_label = 'Acceleration (mm/s*2) '
# label orientation
p.xaxis.major_label_orientation = pi/4
# grid line
p.grid.grid_line_alpha=1.0
# define Y-range
bottom, top = -5, 5
p.y_range=Range1d(bottom, top)
# plot STS1-v data
p.line(x='time', y='acc_mili', source=source,legend="STS1v")
# another Y-axis (right side)
p.extra_y_ranges = {"y2": Range1d(start=-160, end=60)}
# another Y-axis label
p.add_layout(LinearAxis(y_range_name="y2", axis_label='Nanostrain (positive extension)'), 'right')
# plot volumetric strain data with y2 axis
p.line(x='time', y='vol', y_range_name="y2", source=source,legend="Volumetric strain",line_color="black")
# font and etc
p.title.text_font="calibri"
p.title.text_font_size="20pt"
p.title.align="center"
p.xaxis.axis_label_text_font_size = "16pt"
p.xaxis.axis_label_text_font="calibri"
p.xaxis.major_label_text_font_size = "14pt"
p.xaxis.major_label_text_font="calibri"
p.yaxis.axis_label_text_font_size = "16pt"
p.yaxis.axis_label_text_font="calibri"
p.yaxis.major_label_text_font_size = "14pt"
p.yaxis.major_label_text_font="calibri"
p.xaxis.axis_label_text_font_style = 'normal'
p.yaxis.axis_label_text_font_style = 'normal'
p.xaxis[0].formatter.days = '%m/%d/%Y'
# p2 is the second figure
# x_axis is datetime. add tools. define size and add title
p2 = figure(x_axis_type="datetime", tools=TOOLS, plot_width=600, height=600, title = "PFO")
p2.xaxis.axis_label = 'Time (UTC)'
p2.yaxis.axis_label = 'Water level (m) '
p2.xaxis.major_label_orientation = pi/4
p2.grid.grid_line_alpha=1.0
bottom, top = -0.2, 0.2
p2.y_range=Range1d(bottom, top)
p2.line(x='time', y='cib', source=source,legend="CIB")
p2.line(x='time', y='cic', source=source,legend="CIC",line_color="red")
p2.extra_y_ranges = {"y2": Range1d(start=-160, end=60)}
p2.add_layout(LinearAxis(y_range_name="y2", axis_label='Nanostrain (positive extension)'), 'right')
p2.line(x='time', y='vol', y_range_name="y2", source=source,legend="Volumetric strain",line_color="black")
p2.title.text_font="calibri"
p2.title.text_font_size="20pt"
p2.title.align="center"
p2.xaxis.axis_label_text_font_size = "16pt"
p2.xaxis.axis_label_text_font="calibri"
p2.xaxis.major_label_text_font_size = "14pt"
p2.xaxis.major_label_text_font="calibri"
p2.yaxis.axis_label_text_font_size = "16pt"
p2.yaxis.axis_label_text_font="calibri"
p2.yaxis.major_label_text_font_size = "14pt"
p2.yaxis.major_label_text_font="calibri"
p2.xaxis.axis_label_text_font_style = 'normal'
p2.yaxis.axis_label_text_font_style = 'normal'
p2.xaxis[0].formatter.days = '%m/%d/%Y'
# save .html file
output_file("pfo_ts.html", title="PFO")
# if you want to plot only "p"
#show(p) # open a browser
# two figure with "vertical"
#grid = gridplot([[p], [p2]])
# two figure with "horizontal"
grid = gridplot([[p,p2]])
# will open a browser
show(grid)