In [1]:
%matplotlib inline
In [2]:
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
In [3]:
print("# pandas version = ",pd.__version__)

# pandas version =  0.24.2
In [4]:
import bokeh as bk
print("# bokeh version = ",bk.__version__)

# bokeh version =  1.3.4
In [5]:
import scipy as sp
print("# scipy version = ",sp.__version__)

# scipy version =  1.3.0
In [6]:
# ftp directory
pfo_dir = "http://ncedc.org/ftp/outgoing/taira/PFO"
In [7]:
# sts1v data
pfo_sts1v_fi = pfo_dir + "/pfo_sts1v_0.00001-0.01hz.300sps.acc.ts.wh.out"
In [8]:
print(pfo_sts1v_fi)

http://ncedc.org/ftp/outgoing/taira/PFO/pfo_sts1v_0.00001-0.01hz.300sps.acc.ts.wh.out
In [9]:
# pfo volumetric strain data. assuming poisson's ratio = 0.25
pfo_evolume_fi = pfo_dir + "/ez_volume_nu0.25_.ts.wh.out"
In [10]:
print(pfo_evolume_fi)

http://ncedc.org/ftp/outgoing/taira/PFO/ez_volume_nu0.25_.ts.wh.out
In [11]:
# pfo gravitatity data
pfo_grav_fi = pfo_dir + "/PFO.grav.out.ts.wh.out"
In [12]:
# CIB data
cib_fi = pfo_dir + "/CIB.WTHT.ts.wh.out"
In [13]:
# CIC data
cic_fi = pfo_dir + "/CIC.WTHT.ts.wh.out"
In [ ]:
In [14]:
# 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)
In [15]:
# convert "time" to datetime. need for plot
pfo_sts1v['time'] = pd.to_datetime(pfo_sts1v['time'])
In [16]:
# check
pfo_sts1v.head()
Out[16]:
time acc_nano
0 2009-07-14 00:00:00 -53298
1 2009-07-14 00:05:00 -52064
2 2009-07-14 00:10:00 -50835
3 2009-07-14 00:15:00 -49611
4 2009-07-14 00:20:00 -48391
In [17]:
# check
pfo_sts1v.tail()
Out[17]:
time acc_nano
20155 2009-09-21 23:35:00 -2343
20156 2009-09-21 23:40:00 -1022
20157 2009-09-21 23:45:00 382
20158 2009-09-21 23:50:00 1873
20159 2009-09-21 23:55:00 3452
In [18]:
# 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)
In [19]:
pfo_evolume.head()
Out[19]:
time vol ez0 ez90
0 2009-01-01T00:00:00.0000 -0.808287 3.45042 -4.66285
1 2009-01-01T00:05:00.0000 -1.087101 3.06842 -4.69907
2 2009-01-01T00:10:00.0000 -1.375547 2.67398 -4.73730
3 2009-01-01T00:15:00.0000 -1.672854 2.26856 -4.77784
4 2009-01-01T00:20:00.0000 -1.979714 1.85230 -4.82187
In [20]:
# time format conversion
pfo_evolume['time'] = pd.to_datetime(pfo_evolume['time'])
In [21]:
# 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)
In [22]:
# microgal2 is Nan.... Not sure why but microgal2 is not used
pfo_grav.head()
Out[22]:
time microgal microgal2
0 2009-01-01T00:00:00.0000 -4.25255 NaN
1 2009-01-01T00:05:00.0000 -5.57433 NaN
2 2009-01-01T00:10:00.0000 -6.94214 NaN
3 2009-01-01T00:15:00.0000 -8.35256 NaN
4 2009-01-01T00:20:00.0000 -9.80871 NaN
In [23]:
# time format conversion
pfo_grav['time'] = pd.to_datetime(pfo_grav['time'])
In [24]:
# check
pfo_grav['time'].head()
Out[24]:
0   2009-01-01 00:00:00
1   2009-01-01 00:05:00
2   2009-01-01 00:10:00
3   2009-01-01 00:15:00
4   2009-01-01 00:20:00
Name: time, dtype: datetime64[ns]
In [25]:
# CIB data
# time counts meters
cib =pd.read_csv(cib_fi,
                          sep=" ",names=["time", "counts","meters"],
                          header=None, skiprows=1)
In [26]:
# check
cib.head()
Out[26]:
time counts meters
0 2009-01-01T00:00:00.0000 -8300 -0.284441
1 2009-01-01T00:05:00.0000 -8292 -0.284167
2 2009-01-01T00:10:00.0000 -8281 -0.283790
3 2009-01-01T00:15:00.0000 -8283 -0.283858
4 2009-01-01T00:20:00.0000 -8275 -0.283584
In [27]:
# time format conversion
cib['time'] = pd.to_datetime(cib['time'])
In [28]:
# check
cib['time'].head()
Out[28]:
0   2009-01-01 00:00:00
1   2009-01-01 00:05:00
2   2009-01-01 00:10:00
3   2009-01-01 00:15:00
4   2009-01-01 00:20:00
Name: time, dtype: datetime64[ns]
In [29]:
# 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")   ]
In [30]:
# check
cib_select.head()
Out[30]:
time counts meters
55872 2009-07-14 00:00:00 -11848 -0.406031
55873 2009-07-14 00:05:00 -11854 -0.406237
55874 2009-07-14 00:10:00 -11841 -0.405791
55875 2009-07-14 00:15:00 -11840 -0.405757
55876 2009-07-14 00:20:00 -11839 -0.405723
In [31]:
# CIC data
# time counts meters
cic =pd.read_csv(cic_fi,
                          sep=" ",names=["time", "counts","meters"],
                          header=None, skiprows=1)
In [32]:
# check
cic.head()
Out[32]:
time counts meters
0 2009-01-01T00:00:00.0000 -21896 -0.764170
1 2009-01-01T00:05:00.0000 -21881 -0.763647
2 2009-01-01T00:10:00.0000 -21862 -0.762984
3 2009-01-01T00:15:00.0000 -21853 -0.762670
4 2009-01-01T00:20:00.0000 -21835 -0.762042
In [33]:
# time format conversion
cic['time'] = pd.to_datetime(cic['time'])
In [34]:
# check
cic['time'].head()
Out[34]:
0   2009-01-01 00:00:00
1   2009-01-01 00:05:00
2   2009-01-01 00:10:00
3   2009-01-01 00:15:00
4   2009-01-01 00:20:00
Name: time, dtype: datetime64[ns]
In [35]:
# trim data
cic_select = cic[          ("2009-07-14 00:00:00" <= cic['time'])  &   (cic['time'] <= "2009-09-21 23:55:00")   ]
In [36]:
# check
cic_select.head()
Out[36]:
time counts meters
55872 2009-07-14 00:00:00 -23635 -0.824862
55873 2009-07-14 00:05:00 -23634 -0.824827
55874 2009-07-14 00:10:00 -23614 -0.824129
55875 2009-07-14 00:15:00 -23610 -0.823989
55876 2009-07-14 00:20:00 -23599 -0.823605
In [37]:
# 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")   ]
In [38]:
# check
pfo_grav_select.head()
Out[38]:
time microgal microgal2
55872 2009-07-14 00:00:00 40.30138 NaN
55873 2009-07-14 00:05:00 39.36400 NaN
55874 2009-07-14 00:10:00 38.35928 NaN
55875 2009-07-14 00:15:00 37.28634 NaN
55876 2009-07-14 00:20:00 36.13591 NaN
In [39]:
# 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")   ]
In [40]:
# check
pfo_evolume_select.head()
Out[40]:
time vol ez0 ez90
55872 2009-07-14 00:00:00 8.590858 9.63600 3.25028
55873 2009-07-14 00:05:00 8.393278 9.48848 3.10143
55874 2009-07-14 00:10:00 8.181044 9.32478 2.94678
55875 2009-07-14 00:15:00 7.953931 9.14435 2.78654
55876 2009-07-14 00:20:00 7.709951 8.94617 2.61875
In [41]:
# dataframe to numpy matrix
evolume_matrix = pfo_evolume_select.as_matrix()

/anaconda3/envs/netops/lib/python3.7/site-packages/ipykernel_launcher.py:2: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
In [42]:
# check to see evol
evolume_matrix[:,1] # 0 time, 1 evol # 2 ez0 #3 ez90
Out[42]:
array([8.590858, 8.393278, 8.181044, ..., -6.06829, -6.609242999999999,
       -7.151617], dtype=object)
In [43]:
# now add evol data as pfo_sts1v['vol']
pfo_sts1v['vol'] =evolume_matrix[:,1]
In [44]:
# check
pfo_sts1v['vol'].head()
Out[44]:
0    8.59086
1    8.39328
2    8.18104
3    7.95393
4    7.70995
Name: vol, dtype: object
In [45]:
# check
pfo_sts1v['time'].head()
Out[45]:
0   2009-07-14 00:00:00
1   2009-07-14 00:05:00
2   2009-07-14 00:10:00
3   2009-07-14 00:15:00
4   2009-07-14 00:20:00
Name: time, dtype: datetime64[ns]
In [46]:
# dataframe to numpy matrix
grav_matrix = pfo_grav_select.as_matrix()

/anaconda3/envs/netops/lib/python3.7/site-packages/ipykernel_launcher.py:2: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
In [47]:
# check 
grav_matrix[:,1] # microgal
Out[47]:
array([40.30138, 39.364, 38.35928, ..., -29.45913, -32.00622, -34.55945],
      dtype=object)
In [48]:
# now add microgal data as pfo_sts1v['microgal']
pfo_sts1v['microgal'] = grav_matrix[:,1]
In [49]:
# check
pfo_sts1v['microgal'].head()
Out[49]:
0    40.3014
1     39.364
2    38.3593
3    37.2863
4    36.1359
Name: microgal, dtype: object
In [50]:
# dataframe to numpy matrix
cib_matrix = cib_select.as_matrix()

/anaconda3/envs/netops/lib/python3.7/site-packages/ipykernel_launcher.py:2: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
In [51]:
# check 
cib_matrix[:,2] # 0 time, 1 counts, 2 meters
Out[51]:
array([-0.406031, -0.406237, -0.405791, ..., -0.45832700000000004,
       -0.45795, -0.45747], dtype=object)
In [52]:
# now add cib data as pfo_sts1v['cib']
pfo_sts1v['cib'] =  detrend(cib_matrix[:,2])
In [53]:
# dataframe to numpy matrix
cic_matrix = cic_select.as_matrix()

/anaconda3/envs/netops/lib/python3.7/site-packages/ipykernel_launcher.py:2: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
In [54]:
# check
cic_matrix[:,2] # 0 time, 1 counts, 2 meters
Out[54]:
array([-0.8248620000000001, -0.8248270000000001, -0.8241290000000001, ...,
       -0.8681030000000001, -0.867195, -0.866323], dtype=object)
In [55]:
# now add cib data as pfo_sts1v['cic']
pfo_sts1v['cic'] =  detrend(cic_matrix[:,2])
In [ ]:
In [ ]:
In [56]:
#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)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: