In [1]:
matplotlib inline
In [2]:
from __future__ import print_function
%matplotlib inline
import matplotlib.pyplot as plt
#plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 11, 4
plt.rcParams['figure.figsize'] = 19, 6
plt.rcParams['figure.figsize'] = 22, 6
In [3]:
SMALL_SIZE = 22
MEDIUM_SIZE = 24
BIGGER_SIZE = 32
plt.rc('font', size=SMALL_SIZE) # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
In [4]:
import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy import interpolate
from scipy.signal import decimate
from obspy import UTCDateTime
In [5]:
import scipy as sp
import obspy as ob
print("# numpy version = ",np.__version__)
print("# scipy version = ",sp.__version__)
print("# pandas version = ",pd.__version__)
print("# obspy version = ",ob.__version__)
In [6]:
def heaviside(x):
return 0.5*(np.sign(x) + 1)
In [7]:
import math
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
In [8]:
def decimal_year(time):
"""
Return (floating point) decimal year representation of UTCDateTime
input value
"""
start_of_year = UTCDateTime(time.year, 1, 1).timestamp
end_of_year = UTCDateTime(time.year + 1, 1, 1).timestamp
timestamp = time.timestamp
year_fraction = ((timestamp - start_of_year) /
(end_of_year - start_of_year))
return time.year + year_fraction
In [9]:
def year_decimal(start):
year = int(start)
rem = start - year
base = datetime(year, 1, 1)
result = base + timedelta(seconds=(base.replace(year=base.year + 1) - base).total_seconds() * rem)
#print(result)
return result
In [10]:
## inc long term trend
def func(x, e1, e2, f, A, B, e3, e4):
return A + B*x+e1*np.sin(2*math.pi*f*x) + e2*np.cos(2*math.pi*f*x) + \
e3*np.sin(4*math.pi*f*x) + e4*np.cos(4*math.pi*f*x)
In [11]:
## inc long term trend
def func_eq(x, e1, e2, f, A, B, e3, e4, C, D, E, eqT):
return A + B*x+e1*np.sin(2*math.pi*f*x) + e2*np.cos(2*math.pi*f*x) + \
e3*np.sin(4*math.pi*f*x) + e4*np.cos(4*math.pi*f*x) + \
(C+ (D*np.exp(-(x-eqT)/E) ) )*heaviside(x-eqT)
In [12]:
## inc long term trend
def func_twoeq(x, e1, e2, f, A, B, e3, e4, C, D, E, eqT, C2, D2, E2, eqT2):
return A + B*x+e1*np.sin(2*math.pi*f*x) + e2*np.cos(2*math.pi*f*x) + \
e3*np.sin(4*math.pi*f*x) + e4*np.cos(4*math.pi*f*x) + \
(C+ (D*np.exp(-(x-eqT)/E) ) )*heaviside(x-eqT) + \
(C2+ (D2*np.exp(-(x-eqT2)/E2) ) )*heaviside(x-eqT2)
In [13]:
kOPT = 1 # Kik-net
In [14]:
smooth_tw = "5D" # 5days
smooth_tw = "30D" # 30 days
In [15]:
sta = "IBUH03"# Hokkaido velocoty degrease, change direction, no change anino coeff. no increase aniso rms_coeff
In [16]:
if kOPT:
aniso_fi = "/Users/taira/Desktop/Hokkaido/aniso_out/"+sta+".out3"
#rmedian_center.10.IBUH03.out3
#aniso_fi = "http://ncedc.org/ftp/outgoing/taira/Hokkaido/aniso_out/rmedian_center.20."+sta+".out3"
print(aniso_fi)
In [17]:
if kOPT:
aniso_data = pd.read_csv(aniso_fi,
sep=" ",names=["sta", "timeUTC", "lat", "long", "depth", "mag", "elapse_diff", "evid", "viso", "vfast", "vslow", "azfast", "azslow", "azcoeff", "rms_coeff", "leng", "ddeg", "ns", "ne", "f1", "f2", "viso_pvel","vpvs","elapse_days"],header=None)
#aniso_data = pd.read_csv(aniso_fi,
# sep=" ",names=["timeUTC", "viso"],header=None)
aniso_data['time'] = pd.to_datetime(aniso_data['timeUTC'])
data_select = aniso_data.copy()
data_select.set_index('time', inplace=True)
#print(data_select.index)
#ch01_data_select = ch01_data_select[start_timeUTC:end_timeUTC]
aniso_median = data_select.resample(smooth_tw, label='right').median()
aniso_median = aniso_median.dropna()
#print(aniso_median)
In [18]:
if kOPT:
print(aniso_fi)
print(aniso_data.time[1])
In [19]:
if kOPT:
decimal_time_aniso = np.zeros(len(aniso_median.index))
decimal_time_lapse_aniso = np.zeros(len(aniso_median.index))
decimal_ref_time_aniso = 2002.00000
for i in range(len(aniso_median.index)):
#print(i)
dt = UTCDateTime(aniso_median.index[i])
decimal_time_aniso[i] = decimal_year(dt)
decimal_time_lapse_aniso = decimal_time_aniso - decimal_ref_time_aniso
In [20]:
len(decimal_time_lapse_aniso)
Out[20]:
In [21]:
if kOPT:
iburiEQ = pd.DataFrame({'time': ["2018-09-05 18:07:58", "2018-09-05 18:07:58"],
'val': [-9999999999, 9999999999999]})
tokachiEQ = pd.DataFrame({'time': ["2003-09-25 19:50:06", "2003-09-25 19:50:06"],
'val': [-9999999999, 9999999999999]})
#2010-04-04 22:40:42 UTC 32.286°N 115.295°W 10.0 km depth
IBeq = UTCDateTime("2018-09-05T18:07:58")
print(IBeq)
decimal_time_IBeq = decimal_year(IBeq)
decimal_ref_time = 2002.00000
decimal_time_lapse_IBeq = decimal_time_IBeq - decimal_ref_time_aniso
print (decimal_time_lapse_IBeq)
TKeq = UTCDateTime("2003-09-25T19:50:06")
print(TKeq)
decimal_time_TKeq = decimal_year(TKeq)
decimal_ref_time = 2002.00000
decimal_time_lapse_TKeq = decimal_time_TKeq - decimal_ref_time_aniso
print (decimal_time_lapse_TKeq)
In [22]:
if kOPT:
#plt.plot(decimal_time_lapse, data.dvv)
#plt.plot(decimal_time_lapse_aniso, aniso_data.viso*1)
plt.plot(decimal_time_lapse_aniso, aniso_median.viso*1)
plt.xlabel('Year from 2002')
plt.ylabel('Viso (m/s)')
plt.title('Viso Station: '+sta)
plt.grid(True)
In [23]:
print(decimal_time_lapse_aniso[0])
print(decimal_time_lapse_aniso[-1])
In [24]:
st_time = decimal_time_lapse_aniso[0] + 0.01
et_time = decimal_time_lapse_aniso[-1] - 0.01
#st_time = 0.65
#et_time = 5.0
dttime = 0.01
In [25]:
if kOPT:
#f = interpolate.interp1d(decimal_time_lapse_aniso, aniso_data.viso*1, kind="cubic")
f = interpolate.interp1d(decimal_time_lapse_aniso, aniso_median.viso*1, kind="cubic")
#f = interpolate.interp1d(decimal_time_lapse_aniso, aniso_data.viso*1, kind="slinear")
xnew_viso = np.arange(st_time, et_time, dttime)
ynew_viso = f(xnew_viso) # use interpolation function returned by `interp1d
In [26]:
xnew_viso
Out[26]:
In [27]:
len(ynew_viso)
Out[27]:
In [28]:
if kOPT:
#plt.plot(decimal_time_lapse_aniso, aniso_data.viso*1, 'o', xnew_viso, ynew_viso, '-')
plt.plot(decimal_time_lapse_aniso, aniso_median.viso*1, 'o' , label='observed data after medinan filter with '+smooth_tw)
plt.plot(xnew_viso, ynew_viso, '-' , label='interpolated data with cubic spline')
# plt.plot(decimal_time_lapse_aniso, aniso_median.viso*1, 'o' label='data', xnew_viso, ynew_viso, '-')
#plt.plot(xnew_viso, ynew_viso, '-')
plt.ylim((210,240))
plt.xlabel('Year from 2002')
plt.ylabel('Viso (m/s)')
plt.title('Viso Station: '+sta)
plt.grid(True)
plt.legend(loc='upper left')
In [29]:
if kOPT:
xnew_fit = xnew_viso
ynew_fit = ynew_viso
In [30]:
from datetime import datetime, timedelta
In [31]:
def year_decimal(start):
year = int(start)
rem = start - year
base = datetime(year, 1, 1)
result = base + timedelta(seconds=(base.replace(year=base.year + 1) - base).total_seconds() * rem)
#print(result)
return result
In [32]:
#xnew_fitOUT = xnew_fit + 2008
xnew_fitOUT = xnew_fit + 2002
In [33]:
print(xnew_fitOUT[0])
In [34]:
#year_time = np.zeros(len(decimal_time))
#year_time = (len(decimal_time))
year_time = []
#decimal_time_lapse = np.zeros(len(time))
for i in range(len(xnew_fitOUT)):
#print(i)
yt_tmp = year_decimal(xnew_fitOUT[i])
yt_tmp2 = yt_tmp.strftime("%Y-%m-%dT%H:%M:%S.%f")
#decimal_time[i] = decimal_year(dt)
#year_time.append(yt_tmp)
year_time.append(yt_tmp2)
#decimal_time_lapse = decimal_time - decimal_ref_time
In [35]:
print(year_time[0])
In [36]:
year_timedf = pd.DataFrame(year_time)
In [37]:
if kOPT:
eqT_low = decimal_time_lapse_TKeq -0.00001
eqT_high = decimal_time_lapse_TKeq +0.00001
eqT_low2 = decimal_time_lapse_IBeq -0.00001
eqT_high2 = decimal_time_lapse_IBeq +0.00001
In [38]:
print(decimal_time_lapse_TKeq)
print(decimal_time_lapse_IBeq)
In [39]:
twoEQOPT = 1 # if you wante to model two events
#twoEQOPT = 0 # only Tokachi-Oki event
if kOPT:
# D = -0.251 B = 0.04
#def func_eq(x, e1, e2, f, A, B, e3, e4, C, D, E, eqT):
## inc long term trend
#def func_twoeq(x, e1, e2, f, A, B, e3, e4, C, D, E, eqT, C2, D2, E2, eqT2):
# return A + B*x+e1*np.sin(2*math.pi*f*x) + e2*np.cos(2*math.pi*f*x) + \
# e3*np.sin(4*math.pi*f*x) + e4*np.cos(4*math.pi*f*x) + \
# (C+ (D*np.exp(-(x-eqT)/E) ) )*heaviside(x-eqT) + \
# (C2+ (D2*np.exp(-(x-eqT2)/E2) ) )*heaviside(x-eqT2)
if twoEQOPT:
param_bounds_eq=([-np.inf, -np.inf, 0.999999, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf, 0, eqT_low, -np.inf, -np.inf, 0, eqT_low2],
[ np.inf, np.inf, 1.000001, np.inf, np.inf, np.inf, np.inf, np.inf, 0, np.inf, eqT_high, np.inf, 0, np.inf, eqT_high2])
#param_bounds_eq=([-np.inf,-np.inf,0.999999,-np.inf, 0.03997,-np.inf,-np.inf, -0.0000001, -0.26, 0, eqT_low, -0.0000001, -0.120, 0, eqT_low2 ],
# [np.inf, np.inf,1.000001, np.inf, 0.042, np.inf, np.inf, 0.0000001, -0.251, np.inf,eqT_high, 0.0000001, -0.101, np.inf,eqT_high2 ])
popt, pcov = curve_fit(func_twoeq, xnew_fit, ynew_fit, p0=(1.0,0.2,1.0,0.1, 0.1,0.1, 0.1,0.0, -1.0, 1.0, decimal_time_lapse_TKeq,0.0,-1.0,1.0,decimal_time_lapse_IBeq), bounds=param_bounds_eq)
else:
param_bounds_eq=([ -np.inf, -np.inf, 0.999999, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf, 0, eqT_low],
[ np.inf, np.inf, 1.000001, np.inf, np.inf, np.inf, np.inf, np.inf, 0, np.inf, eqT_high])
popt, pcov = curve_fit(func_eq, xnew_fit, ynew_fit, p0=(1.0,0.2,1.0,0.1, 0.1,0.1, 0.1,0.0, -1.0, 1.0,decimal_time_lapse_TKeq), bounds=param_bounds_eq)
In [ ]:
In [ ]:
In [ ]:
In [40]:
print (popt)
print (pcov)
#Print results
print("e1 = %.5f +/- %.5f" % (popt[0], math.sqrt(pcov[0, 0])))
print("e2 = %.5f +/- %.5f" % (popt[1], math.sqrt(pcov[1, 1])))
print("f = %.5f +/- %.5f" % (popt[2], math.sqrt(pcov[2, 2])))
print("A = %.5f +/- %.5f" % (popt[3], math.sqrt(pcov[3, 3])))
print("B = %.5f +/- %.5f" % (popt[4], math.sqrt(pcov[4, 4])))
print("e3 = %.5f +/- %.5f" % (popt[5], math.sqrt(pcov[5, 5])))
print("e4 = %.5f +/- %.5f" % (popt[6], math.sqrt(pcov[6, 6])))
In [41]:
print("C = %.5f +/- %.5f" % (popt[7], math.sqrt(pcov[7, 7])))
print("D = %.5f +/- %.5f" % (popt[8], math.sqrt(pcov[8, 8])))
print("E = %.5f +/- %.5f" % (popt[9], math.sqrt(pcov[9, 9])))
print("eqT = %.5f +/- %.5f" % (popt[10], math.sqrt(pcov[10, 10])))
In [42]:
if twoEQOPT == 1:
#if kOPT == 1:
print("C2 = %.5f +/- %.5f" % (popt[11], math.sqrt(pcov[11, 11])))
print("D2 = %.5f +/- %.5f" % (popt[12], math.sqrt(pcov[12, 12])))
print("E2 = %.5f +/- %.5f" % (popt[13], math.sqrt(pcov[13, 13])))
print("eqT2 = %.5f +/- %.5f" % (popt[14], math.sqrt(pcov[14, 14])))
In [43]:
perr = np.sqrt(np.diag(pcov))
e1=popt[0]
e2=popt[1]
f=popt[2]
A=popt[3]
B=popt[4]
e3=popt[5]
e4=popt[6]
e1sd = math.sqrt(pcov[0, 0])
e2sd = math.sqrt(pcov[1, 1])
fsd = math.sqrt(pcov[2, 2])
Asd = math.sqrt(pcov[3, 3])
Bsd = math.sqrt(pcov[4, 4])
e3sd = math.sqrt(pcov[5, 5])
e4sd = math.sqrt(pcov[6, 6])
#e1 = 0.5
#e2 = 0.5
phi=math.atan2(1*e1,e2)
phi_deg=(phi*180.0)/math.pi
phi_deg2=360+phi_deg
phi_day = (phi_deg / 360.0) *366.0
e1_2=e1*e1
e2_2=e2*e2
etmp=e1_2+e2_2
e=math.sqrt(etmp)
phi4p=math.atan2(1*e3,e4)
phi4p_deg=(phi4p*180.0)/math.pi
phi4p_deg2=360+phi4p_deg
phi4p_day = (phi4p_deg / 360.0) *366.0
e3_2=e3*e3
e4_2=e4*e4
etmp4p=e3_2+e4_2
e4p=math.sqrt(etmp4p)
e_ratio = e4p / e
print("e1 = %.5f "% (e1))
print("e2 = %.5f "% (e2))
print("e = %.5f "% (e))
print("f(1/year) = %.5f "% (f))
print("A = %.5f "% (A))
print("B = %.5f "% (B))
print("e1sd = %.5f "% (e1sd))
print("e2sd = %.5f "% (e2sd))
print("f(1/year)sd = %.5f "% (fsd))
print("Asd = %.5f "% (Asd))
print("Bsd = %.5f "% (Bsd))
print("phi = %.5f "% (phi))
print("phi_deg = %.5f "% (phi_deg))
print("phi_deg2 = %.5f "% (phi_deg2))
print("e3 = %.5f "% (e3))
print("e4 = %.5f "% (e4))
print("e3sd = %.5f "% (e3sd))
print("e4sd = %.5f "% (e4sd))
print("e4p = %.5f "% (e4p))
print("phi4p = %.5f "% (phi4p))
print("phi4p_deg = %.5f "% (phi4p_deg))
print("phi4p_deg2 = %.5f "% (phi4p_deg2))
print("phi_day = %.5f "% (phi_day))
print("phi4p_day = %.5f "% (phi4p_day))
print("e_ratio = %.5f "% (e_ratio))
# cos (w - phi)
In [44]:
C=popt[7]
D=popt[8]
E=popt[9]
eqT=popt[10]
Csd = math.sqrt(pcov[7, 7])
Dsd = math.sqrt(pcov[8, 8])
Esd = math.sqrt(pcov[9, 9])
eqTsd = math.sqrt(pcov[10, 10])
print("C = %.5f "% (C))
print("D = %.5f "% (D))
print("E = %.5f "% (E))
print("eqT = %.5f "% (eqT))
print("Csd = %.5f "% (Csd))
print("Dsd = %.5f "% (Dsd))
print("Esd = %.5f "% (Esd))
print("eqTsd = %.5f "% (eqTsd))
In [45]:
if twoEQOPT == 1:
#if kOPT == 1:
C2=popt[11]
D2=popt[12]
E2=popt[13]
eqT2=popt[14]
Csd2 = math.sqrt(pcov[7, 7])
Dsd2 = math.sqrt(pcov[8, 8])
Esd2 = math.sqrt(pcov[9, 9])
eqTsd2 = math.sqrt(pcov[10, 10])
print("C2 = %.5f "% (C2))
print("D2 = %.5f "% (D2))
print("E2 = %.5f "% (E2))
print("eqT2 = %.5f "% (eqT2))
print("Csd2 = %.5f "% (Csd2))
print("Dsd2 = %.5f "% (Dsd2))
print("Esd2 = %.5f "% (Esd2))
print("eqTsd2 = %.5f "% (eqTsd2))
In [46]:
#print (perr) # one sigma for each parameter
In [47]:
#def func_eq(x, e1, e2, f, A, B, e3, e4, C, D, E, eqT):
## inc long term trend
#def func_eq(x, e1, e2, f, A, B, e3, e4, C, D, E, eqT):
# return A + B*x+e1*np.sin(2*math.pi*f*x) + e2*np.cos(2*math.pi*f*x) + \
# e3*np.sin(4*math.pi*f*x) + e4*np.cos(4*math.pi*f*x) + \
# (C+D*np.exp(-(x-eqT)/E))*heaviside(x-eqT)
In [48]:
resi_fit = np.zeros(len(xnew_fit))
yOUT_fit = np.zeros(len(ynew_fit))
yOUT_fit_1st = np.zeros(len(ynew_fit))
yOUT_fit_2nd = np.zeros(len(ynew_fit))
for i in range(len(xnew_fit)):
# print xdata[i]
#yOUT_fit[i]= A + B*xnew_fit[i]+e1*np.sin(2*math.pi*f*xnew_fit[i]) + e2*np.cos(2*math.pi*f*xnew_fit[i])
yOUT_fit[i]= A + B*xnew_fit[i]+e1*np.sin(2*math.pi*f*xnew_fit[i]) +\
e2*np.cos(2*math.pi*f*xnew_fit[i])+\
e3*np.sin(4*math.pi*f*xnew_fit[i]) + \
e4*np.cos(4*math.pi*f*xnew_fit[i])+ \
(C+ (D*np.exp(-(xnew_fit[i]-eqT)/E) ) ) *heaviside(xnew_fit[i]-eqT)
#print yOUT[i]
resi_fit[i]=ynew_fit[i]-yOUT_fit[i]
#print resi[i]
yOUT_fit_1st[i]= A + B*xnew_fit[i]+e1*np.sin(2*math.pi*f*xnew_fit[i]) +\
e2*np.cos(2*math.pi*f*xnew_fit[i])
yOUT_fit_2nd[i]= A + B*xnew_fit[i]+e3*np.sin(4*math.pi*f*xnew_fit[i]) +\
e4*np.cos(4*math.pi*f*xnew_fit[i])
In [49]:
if twoEQOPT == 1:
#if kOPT == 1:
resi_fit = np.zeros(len(xnew_fit))
yOUT_fit = np.zeros(len(ynew_fit))
for i in range(len(xnew_fit)):
# print xdata[i]
#yOUT_fit[i]= A + B*xnew_fit[i]+e1*np.sin(2*math.pi*f*xnew_fit[i]) + e2*np.cos(2*math.pi*f*xnew_fit[i])
yOUT_fit[i]= A + B*xnew_fit[i]+e1*np.sin(2*math.pi*f*xnew_fit[i]) +\
e2*np.cos(2*math.pi*f*xnew_fit[i])+\
e3*np.sin(4*math.pi*f*xnew_fit[i]) + \
e4*np.cos(4*math.pi*f*xnew_fit[i])+ \
(C+ (D*np.exp(-(xnew_fit[i]-eqT)/E) ) ) *heaviside(xnew_fit[i]-eqT) + \
(C2+ (D2*np.exp(-(xnew_fit[i]-eqT2)/E2) ) ) *heaviside(xnew_fit[i]-eqT2)
#print yOUT[i]
resi_fit[i]=ynew_fit[i]-yOUT_fit[i]
#print resi[i]
In [50]:
resi_sq= np.sum(resi_fit**2)
print(resi_sq)
data_sq= np.sum(ynew_fit**2)
print(data_sq)
VR = (1 - (resi_sq/data_sq)) * 100
print(VR)
In [51]:
type(f)
Out[51]:
In [52]:
s = "e1 = %.5f "% (e1)
print(s)
In [53]:
VROUTOPT = 0
if VROUTOPT:
vr_fi = "VR.out_rmedian_center."+str(tw)+netIn+"_"+staIn+"_"+netIn+"_"+staIn+"."+xccomIn+"."+freqIn+".030.dttwithoffset.static_1.0_5.0_10.0_both_0.85_0.1_0.5.out"
if dvv_2sdOPT:
#dvv2sd_"+str(dvv_2sd_max)+"."
vr_fi = "VR.out_rmedian_center."+str(tw)+"_dvv2sd_"+str(dvv_2sd_max)+"."+netIn+"_"+staIn+"_"+netIn+"_"+staIn+"."+xccomIn+"."+freqIn+".030.dttwithoffset.static_1.0_5.0_10.0_both_0.85_0.1_0.5.out"
fi = open(vr_fi, 'w')
#s = 'The value of x is ' + repr(x) + ', and y is ' + repr(y) + '...'
s = "# VR = "+str(VR)+"\n"
fi.write(s)
s = "# st_time = "+str(st_time)+"\n"
fi.write(s)
s = "# et_time = "+str(et_time)+"\n"
fi.write(s)
s = "# dttime = "+str(dttime)+"\n"
fi.write(s)
s = "# e1 = %.5f "% (e1)+"\n"
fi.write(s)
s = "# e2 = %.5f "% (e2)+"\n"
fi.write(s)
s = "# e = %.5f "% (e)+"\n"
fi.write(s)
s = "# f = %.5f "% (f)+"\n"
fi.write(s)
s = "# A = %.5f "% (A)+"\n"
fi.write(s)
s = "# B = %.5f "% (B)+"\n"
fi.write(s)
s = "# Bsd = %.5f "% (Bsd)+"\n"
fi.write(s)
s = "# C = %.5f "% (C)+"\n"
fi.write(s)
s = "# D = %.5f "% (D)+"\n"
fi.write(s)
s = "# E = %.5f "% (E)+"\n"
fi.write(s)
s = "# eqT = %.5f "% (eqT)+"\n"
fi.write(s)
if twoEQOPT == 1:
#if kOPT == 1:
s = "# C2 = %.5f "% (C2)+"\n"
fi.write(s)
s = "# D2 = %.5f "% (D2)+"\n"
fi.write(s)
s = "# E2 = %.5f "% (E2)+"\n"
fi.write(s)
s = "# eqT2 = %.5f "% (eqT2)+"\n"
fi.write(s)
s = "# e1sd = %.5f "% (e1sd)+"\n"
fi.write(s)
s = "# e2sd = %.5f "% (e2sd)+"\n"
fi.write(s)
s = "# fsd = %.5f "% (fsd)+"\n"
fi.write(s)
s = "# Asd = %.5f "% (Asd)+"\n"
fi.write(s)
s = "# Csd = %.5f "% (Csd)+"\n"
fi.write(s)
s = "# Dsd = %.5f "% (Dsd)+"\n"
fi.write(s)
s = "# Esd = %.5f "% (Esd)+"\n"
fi.write(s)
s = "# eqTsd = %.5f "% (eqTsd)+"\n"
fi.write(s)
if twoEQOPT == 1:
#if kOPT == 1:
s = "# Csd2 = %.5f "% (Csd2)+"\n"
fi.write(s)
s = "# Dsd2 = %.5f "% (Dsd2)+"\n"
fi.write(s)
s = "# Esd2 = %.5f "% (Esd2)+"\n"
fi.write(s)
s = "# eqTsd2 = %.5f "% (eqTsd2)+"\n"
fi.write(s)
s = "# phi = %.5f "% (phi)+"\n"
fi.write(s)
s = "# phi_deg = %.5f "% (phi_deg)+"\n"
fi.write(s)
s = "# phi_deg2 = %.5f "% (phi_deg2)+"\n"
fi.write(s)
s = "# e3 = %.5f "% (e3)+"\n"
fi.write(s)
s = "# e4 = %.5f "% (e4)+"\n"
fi.write(s)
s = "# e3sd = %.5f "% (e3sd)+"\n"
fi.write(s)
s = "# e4sd = %.5f "% (e4sd)+"\n"
fi.write(s)
s = "# e4p = %.5f "% (e4p)+"\n"
fi.write(s)
s = "# phi4p = %.5f "% (phi4p)+"\n"
fi.write(s)
s = "# phi4p_deg = %.5f "% (phi4p_deg)+"\n"
fi.write(s)
s = "# phi4p_deg2 = %.5f "% (phi4p_deg2)+"\n"
fi.write(s)
s = "# phi_day = %.5f "% (phi_day)+"\n"
fi.write(s)
s = "# phi4p_day = %.5f "% (phi4p_day)+"\n"
fi.write(s)
s = "# e_ratio = %.5f "% (e_ratio)+"\n"
fi.write(s)
#f.write('2: 2nd line\n')
#f.write('2: last line\n')
fi.close()
In [54]:
#vr_fi = "VR.out_"+netIn+"_"+staIn+"_"+netIn+"_"+staIn+"."+xccomIn+"."+freqIn+".030.dttwithoffset.static_1.0_5.0_10.0_both_0.85_0.1_0.5.out"
#print(vr_fi)
# store output in pandas dataframe
#df = pd.DataFrame({ 'VR' : VRout
# })
# save
#df.to_csv(vr_fi,index=False)
In [55]:
year_timedf = pd.DataFrame(year_time)
In [56]:
dvvOUTOPT = 0
if dvvOUTOPT:
dvvOUT_fi = "dvvOUT_rmedian_center."+str(tw)+netIn+"_"+staIn+"_"+netIn+"_"+staIn+"."+xccomIn+"."+freqIn+"."+paraIn+".out"
if dvv_2sdOPT:
#dvv2sd_"+str(dvv_2sd_max)+"."
dvvOUT_fi = "dvvOUT_rmedian_center."+str(tw)+"_dvv2sd_"+str(dvv_2sd_max)+"."+netIn+"_"+staIn+"_"+netIn+"_"+staIn+"."+xccomIn+"."+freqIn+"."+paraIn+".out"
print(dvvOUT_fi)
df = pd.DataFrame({ 'ynew_fit' : ynew_fit,
'xnew_fit': xnew_fit,
'yOUT_fit' : yOUT_fit,
'year_time': year_time})
df.to_csv(dvvOUT_fi,index=False)
In [57]:
pngOUTOPT = 0
if pngOUTOPT:
pngOUT_fi = "fitOUT_rmedian_center."+str(tw)+netIn+"_"+staIn+"_"+netIn+"_"+staIn+"."+xccomIn+"."+freqIn+"."+paraIn+".png"
if dvv_2sdOPT:
#dvv2sd_"+str(dvv_2sd_max)+"."
pngOUT_fi = "fitOUT_rmedian_center."+str(tw)+"_dvv2sd_"+str(dvv_2sd_max)+"."+netIn+"_"+staIn+"_"+netIn+"_"+staIn+"."+xccomIn+"."+freqIn+"."+paraIn+".png"
In [58]:
#plt.plot(xnew_fit, yOUT_fit_2nd)
#plt.plot(xnew_fit, yOUT_fit_1st)
plt.plot(xnew_fit, yOUT_fit,label='fitting')
plt.plot(xnew_fit, ynew_fit,label='interpolated data with cubic spline')
plt.xlabel('Year from 2002')
plt.ylabel('Viso (m/s)')
plt.title('Viso Station: '+sta)
plt.grid(True)
plt.legend(loc='upper left')
plt.ylim((210,240))
if pngOUTOPT:
plt.savefig(pngOUT_fi)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: