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 [ ]:
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 [ ]:
In [7]:
# a -> eqT
# c == 1000000 -> steep step
# c = 1 smooth step
def step(x, a, c):
return 1/(1+(np.exp(-2*c *(x-a))) )
In [8]:
import math
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
In [9]:
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 [10]:
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 [11]:
## 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 [12]:
## 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 [13]:
## 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 [14]:
def func_eq_step(x, e1, e2, f, A, B, e3, e4, C, D, E, eqT, sc):
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) ) )*step(x,eqT, sc)
In [15]:
stepOPT = 0 # use step func
In [16]:
dvv_2sdOPT = 0 # no info. error
dvv_2sd_max = 0.1 # %
dvv_2sd_max = 0.2 # %
In [17]:
# use 1D
smooth_tw = "1D" # 1
#smooth_tw = "11D" # 1
#smooth_tw = "30D" # 1
In [18]:
sta = "IBU03"
In [19]:
xc_para = "all_event_120-day.info"
xc_para = "all_event_60-day.info"
xc_para = "all_event.info"
xc_para = "all_event_180-day.info"
In [20]:
#dvv_fi = "./"+sta+".ZRT.01.030.dttwithoffset.dynamic_2.0_20.0_30.0_both_0.85_0.1_0.5.out" #modified by an 20200409 lagtop
data_dir = "/Users/taira/work/python_work/Gaussian-stack"
#data_dir = "/Users/taira/work/python_work/event_stacking/event_5"
#data_dir = "/Users/taira/work/python_work/event_stacking/event_11"
#dvv_fi = "./"+sta+"."+xc_para
# test no sta included
dvv_fi = data_dir+"/"+sta+"."+xc_para
dvv_fi = data_dir+"/"+xc_para
print(dvv_fi)
In [21]:
#20020617153100 2002-06-17 15:31:00 216.32 1091.3059 5.0449 0.88 174.80
dvv_data = pd.read_csv(dvv_fi,
sep=" ",names=["eqid", "time1", "time2", "viso", "vp", "vpvs", "aniso_st", "fast_az"],header=None)
#dvv_data['time'] = pd.to_datetime(dvv_data['timeUTC'])
In [22]:
dvv_data['timeUTC'] = dvv_data['time1'].astype(str) + 'T' + dvv_data['time2'].astype(str)
dvv_data['time'] = pd.to_datetime(dvv_data['timeUTC'])
In [23]:
print(dvv_data)
In [24]:
if dvv_2sdOPT == 1:
data_select = dvv_data[ dvv_data['dvv_2sd'] <= dvv_2sd_max ]
else:
data_select = dvv_data.copy()
# old just copy
#
#catalog_data_eq_normalF = catalog_data_eq[catalog_data_eq['fault_type'] == 0]
In [25]:
data_select.set_index('time', inplace=True)
In [26]:
dvv_median = data_select.resample(smooth_tw, label='right').median()
dvv_median = dvv_median.dropna()
In [27]:
print(dvv_fi)
print(dvv_data.time[1])
In [28]:
decimal_time_dvv = np.zeros(len(dvv_median.index))
decimal_time_lapse_dvv = np.zeros(len(dvv_median.index))
decimal_ref_time_dvv = 2002.00000
for i in range(len(dvv_median.index)):
#print(i)
dt = UTCDateTime(dvv_median.index[i])
decimal_time_dvv[i] = decimal_year(dt)
decimal_time_lapse_dvv = decimal_time_dvv - decimal_ref_time_dvv
In [29]:
len(decimal_time_lapse_dvv)
Out[29]:
In [30]:
iburiEQ = pd.DataFrame({'time': ["2018-09-05 18:07:58", "2018-09-05 18:07:58"],
'val': [-9999999999, 9999999999999]})
iburiEQ['time'] = pd.to_datetime(iburiEQ['time'])
IBeq = UTCDateTime("2018-09-05 18:07:58")
tokachiEQ = pd.DataFrame({'time': ["2003-09-25 19:50:06", "2003-09-25 19:50:06"],
'val': [-9999999999, 9999999999999]})
tokachiEQ['time'] = pd.to_datetime(tokachiEQ['time'])
TKeq = UTCDateTime("2003-09-25 19:50:06")
decimal_time_TKeq = decimal_year(TKeq)
decimal_time_lapse_TKeq = decimal_time_TKeq - decimal_ref_time_dvv
print (decimal_time_lapse_TKeq )
decimal_time_IBeq = decimal_year(IBeq)
decimal_time_lapse_IBeq = decimal_time_IBeq - decimal_ref_time_dvv
print (decimal_time_lapse_IBeq )
In [31]:
print(dvv_median)
In [32]:
plt.plot(decimal_time_lapse_dvv, dvv_median.viso*1)
#plt.fill_between(decimal_time_lapse_dvv, (dvv_median.dvv-dvv_median.dvv_2sd),
# (dvv_median.dvv+dvv_median.dvv_2sd), alpha=0.25,linewidth = 2.5, color="orange")
plt.xlabel('Year from 2002')
#plt.ylim((-0.2,0.2))
plt.ylabel('Viso (m/s)') #modified by an 20200407
plt.title('Viso Station: '+sta) #modified by an 20200409
plt.grid(True)
In [ ]:
In [33]:
print(decimal_time_lapse_dvv[0])
print(decimal_time_lapse_dvv[-1])
In [34]:
st_time = decimal_time_lapse_dvv[0] + 0.01
et_time = decimal_time_lapse_dvv[-1] - 0.01
print(st_time)
print(et_time)
#st_time = 0.65
#et_time = 5.0
In [35]:
### adds
dttime = 0.01
#dttime = 0.05
#dttime = 0.1
### adds
In [36]:
f = interpolate.interp1d(decimal_time_lapse_dvv, dvv_median.viso*1, kind="cubic")
xnew_dvv = np.arange(st_time, et_time, dttime)
ynew_dvv = f(xnew_dvv) # use interpolation function returned by `interp1d
In [37]:
len(xnew_dvv)
Out[37]:
In [38]:
len(ynew_dvv)
Out[38]:
In [39]:
plt.plot(decimal_time_lapse_dvv, dvv_median.viso*1, 'o' , label='observed data after medinan filter with '+smooth_tw)
plt.plot(xnew_dvv, ynew_dvv, '-' , 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.ylim((200,230))
plt.legend(loc='upper right')
Out[39]:
In [40]:
xnew_fit = xnew_dvv
ynew_fit = ynew_dvv
In [41]:
print(xnew_fit)
print(ynew_fit)
In [42]:
from datetime import datetime, timedelta
In [43]:
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 [44]:
xnew_fitOUT = xnew_fit + 2002
In [45]:
print(xnew_fitOUT[0])
In [46]:
#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 [47]:
print(year_time[0])
In [48]:
year_timedf = pd.DataFrame(year_time)
In [49]:
print(decimal_time_lapse_IBeq)
In [50]:
# seasonality
e1_min = -np.inf
e1_max = np.inf
e2_min = -np.inf
e2_max = np.inf
# seasonality
e3_min = -np.inf
e3_max = np.inf
e4_min = -np.inf
e4_max = np.inf
# freq. seasonality
f_min = 0.999999
f_max = 1.000001
# offset
A_min = -np.inf
A_max = np.inf
# liner trend
B_min = -np.inf
B_max = np.inf
# no longer used should be 0
C_min = -0.0000001
C_max = 0.0000001
# coseismic change. currently assuming negative value
D_min = -np.inf
D_max = 0
#D_max = np.inf
# recover time -> should be positive
E_min = 0
E_max = np.inf
# eqT
eqT_low = decimal_time_lapse_TKeq - 0.000001
eqT_low = decimal_time_lapse_TKeq - 0.2
eqT_high = decimal_time_lapse_TKeq + 0.2
eqT_high = decimal_time_lapse_TKeq + 0.00001
In [ ]:
In [51]:
# # e1 e2 f A B e3 e4 C D E eqT C2
# ([-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=([e1_min, e2_min, f_min, A_min, B_min, e3_min, e4_min, C_min, D_min, E_min, eqT_low],
[e1_max, e2_max, f_max, A_max, B_max, e3_max, e4_max, C_max, D_max, E_max, eqT_high])
In [52]:
e1_init = 1.0
e2_init = 1.0
e3_init = 1.0
e4_init = 1.0
f_init = 1.0
A_init = 0.1
B_init = 0.1
C_init = 0.0
D_init = -1.0
E_init = 1.0
popt, pcov = curve_fit(func_eq, xnew_fit, ynew_fit, p0=(e1_init, e2_init,f_init, A_init, B_init, e3_init, e4_init, C_init, D_init, E_init, decimal_time_lapse_TKeq), bounds=param_bounds_eq)
In [ ]:
In [53]:
if stepOPT:
sc_init = 5.0
sc_min = 5.0
sc_max = 100
#sc_init = 1.0
#sc_min = 1.0
#sc_max = 100
sc_init = 10.0
sc_min = 10.0
sc_max = 1000
param_bounds_eq=([e1_min, e2_min, f_min, A_min, B_min, e3_min, e4_min, C_min, D_min, E_min, eqT_low, sc_min],
[e1_max, e2_max, f_max, A_max, B_max, e3_max, e4_max, C_max, D_max, E_max, eqT_high, sc_max])
popt, pcov = curve_fit(func_eq_step, xnew_fit, ynew_fit, p0=(e1_init, e2_init,f_init, A_init, B_init, e3_init, e4_init, C_init, D_init, E_init, decimal_time_lapse_TKeq, sc_init), bounds=param_bounds_eq)
In [54]:
## 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)
C2_init = 0.0
D2_init = -1.0
E2_init = 1.0
# no longer used should be 0
C2_min = -0.0000001
C2_max = 0.0000001
# coseismic change. currently assuming negative value
D2_min = -np.inf
D2_max = 0
#D_max = np.inf
# recover time -> should be positive
E2_min = 0
E2_max = np.inf
eqT2_min = decimal_time_lapse_IBeq -0.2
eqT2_max = decimal_time_lapse_IBeq + 0.00001
twoEQOPT = 1
if twoEQOPT == 1:
# D = -0.251 B = 0.04
#def func_eq(x, e1, e2, f, A, B, e3, e4, C, D, E, eqT, C2, D2, E2, eqT2 ):
param_bounds_eq=([e1_min, e2_min, f_min, A_min, B_min, e3_min, e4_min, C_min, D_min, E_min, eqT_low, C2_min, D2_min, E2_min, eqT2_min],
[e1_max, e2_max, f_max, A_max, B_max, e3_max, e4_max, C_max, D_max, E_max, eqT_high, C2_max, D2_max, E2_max, eqT2_max])
popt, pcov = curve_fit(func_twoeq, xnew_fit, ynew_fit, p0=(e1_init, e2_init,f_init, A_init, B_init, e3_init, e4_init, C_init, D_init, E_init, decimal_time_lapse_TKeq, C2_init, D2_init, E2_init, decimal_time_lapse_IBeq), bounds=param_bounds_eq)
In [55]:
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 [56]:
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 [57]:
if twoEQOPT == 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 [58]:
if stepOPT:
#print("sb = %.5f +/- %.5f" % (popt[11], math.sqrt(pcov[11, 11])))
print("sc = %.5f +/- %.5f" % (popt[11], math.sqrt(pcov[11, 11])))
In [59]:
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 [60]:
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 [61]:
if twoEQOPT == 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 [62]:
if stepOPT:
#sb=popt[11]
sc=popt[11]
scsd = math.sqrt(pcov[11, 11])
#print("sb = %.5f "% (sb))
print("sc = %.5f "% (sc))
print("scsd = %.5f "% (scsd))
In [63]:
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) #modified by an 20200421
# 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) ) ) *h1*(heaviside(xnew_fit[i] -eqT)+ h2) #modified by an 20200421
#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 [64]:
if twoEQOPT == 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 [65]:
#def func_eq_step(x, e1, e2, f, A, B, e3, e4, C, D, E, eqT, sb, sc):
# 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) ) )*step(x,eqT,sb,sc)
if stepOPT:
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)):
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) ) ) * (step(xnew_fit[i],eqT, sc)) #modified by an 20200421
#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 [66]:
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 [67]:
type(f)
Out[67]:
In [68]:
s = "e1 = %.5f "% (e1)
print(s)
In [69]:
VROUTOPT = 1
if VROUTOPT:
#vr_fi = "VR.out_rmedian_center."+str(smooth_tw )+netIn+"_"+staIn+"_"+netIn+"_"+staIn+"."+xccomIn+"."+freqIn+".030.dttwithoffset.static_1.0_5.0_10.0_both_0.85_0.1_0.5.out"
vr_fi = "VR.out_"+sta+"_smtw"+str(smooth_tw )+"_"+xc_para
if dvv_2sdOPT:
#dvv2sd_"+str(dvv_2sd_max)+"."
#vr_fi = "VR.out_rmedian_center."+str(smooth_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"
#vr_fi = "VR.out_"+sta+"."+xc_para+str(smooth_tw )+"_dvv2sd_"+str(dvv_2sd_max)
vr_fi = "VR.out_"+sta+"_smtw"+str(smooth_tw )+"_dvv2sd_"+str(dvv_2sd_max)+"_"+xc_para
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 = "# smooth_tw = "+str( smooth_tw )+"\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 = "# 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)
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 = "# Bsd = %.5f "% (Bsd)+"\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:
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)
if stepOPT == 1:
#if kOPT == 1:
s = "# sc = %.5f "% (sc)+"\n"
fi.write(s)
s = "# scsd = %.5f "% (scsd)+"\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 [70]:
#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 [71]:
year_timedf = pd.DataFrame(year_time)
In [72]:
dvvOUTOPT = 0
if dvvOUTOPT:
dvvOUT_fi = "dvvOUT_rmedian_center."+str(smooth_tw )+netIn+"_"+staIn+"_"+netIn+"_"+staIn+"."+xccomIn+"."+freqIn+"."+paraIn+".out"
if dvv_2sdOPT:
#dvv2sd_"+str(dvv_2sd_max)+"."
dvvOUT_fi = "dvvOUT_rmedian_center."+str(smooth_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 [73]:
pngOUTOPT = 1
if pngOUTOPT:
#pngOUT_fi = "fitOUT_rmedian_center."+str(tw)+netIn+"_"+staIn+"_"+netIn+"_"+staIn+"."+xccomIn+"."+freqIn+"."+paraIn+".png"
pngOUT_fi = "fitOUT_"+sta+"_smtw"+str(smooth_tw )+"_"+xc_para+".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"
pngOUT_fi = "fitOUT_"+sta+"_smtw"+str(smooth_tw )+"_dvv2sd_"+str(dvv_2sd_max)+"_"+xc_para+".png"
print(pngOUT_fi)
In [74]:
#plt.plot(xnew_fit, yOUT_fit_2nd)
#plt.plot(xnew_fit, yOUT_fit_1st)
plt.plot(decimal_time_lapse_dvv, dvv_median.viso*1, 'o' , label='observed data after medinan filter with '+smooth_tw)
#plt.plot(decimal_time_lapse_dvv, dvv_median.dvv*1, 'o' , label='observed data after medinan filter with '+smooth_tw+'-'+(str)dvv_2sd_max)
plt.plot(xnew_fit, ynew_fit,label='interpolated data with cubic spline', color="black")
plt.plot(xnew_fit, yOUT_fit,label='fitting', color="red", linewidth = 5)
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((205,230))
if pngOUTOPT:
plt.savefig(pngOUT_fi)
In [75]:
print("VR = %.5f" %VR)
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])))
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])))
if twoEQOPT == 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 [76]:
if stepOPT:
#print("sb = %.5f +/- %.5f" % (popt[11], math.sqrt(pcov[11, 11])))
print("sc = %.5f +/- %.5f" % (popt[11], math.sqrt(pcov[11, 11])))
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: