In [1]:
import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

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__)

# numpy version =  1.19.1
# scipy version =  1.5.2
# pandas version =  1.0.4
# obspy version =  1.2.2
In [2]:
SMALL_SIZE = 16
MEDIUM_SIZE = 18
BIGGER_SIZE = 22

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 [3]:
to_rad=math.pi/180.0
def aniso_fit(x, Viso, Ac, As):
  return   Viso + Ac*np.cos(2*x*to_rad) + As*np.sin(2*x*to_rad)
In [4]:
aniso_dir = "http://ncedc.org/ftp/outgoing/taira/Hokkaido"
aniso_para = "7_5_-1_0_1_15"
In [5]:
#sta = "IBUH01"
sta = "IBUH03"
In [6]:
aniso_fi = aniso_dir+"/VEL_neg_"+sta+"_"+aniso_para+".csv"
aniso_df = pd.read_csv(aniso_fi,
                       sep=",",names=["az", "vel"], header=None)
print("# aniso_fi = ", aniso_fi)

# aniso_fi =  http://ncedc.org/ftp/outgoing/taira/Hokkaido/VEL_neg_IBUH03_7_5_-1_0_1_15.csv
In [7]:
aniso_fi = "./mccc_lag.out"
In [8]:
aniso_fi = aniso_dir+"/"+sta+".mccc_lag.out"
aniso_fi = aniso_dir+"/"+sta+".mccc_lag.out.enve"

aniso_df = pd.read_csv(aniso_fi,
                       sep=" ",names=["az", "vel"], header=None)
print("# aniso_fi = ", aniso_fi)

# aniso_fi =  http://ncedc.org/ftp/outgoing/taira/Hokkaido/IBUH03.mccc_lag.out.enve
In [9]:
#print(aniso_df.head())
In [10]:
#statOUT = aniso_df['vel'].describe()
#print(statOUT)
#print(statOUT['mean'])
#aniso_df['vel'] = aniso_df['vel']/statOUT['mean']
In [11]:
daz = aniso_df['az'][1] - aniso_df['az'][0]
print("# daz = ", daz)

# daz =  5.0
In [12]:
Viso_init = 200
Ac_init = 1
As_init = 1

Viso_min = 0
Viso_max = 30000

Ac_min = -np.inf
Ac_max =  np.inf

As_min = -np.inf
As_max =  np.inf
In [13]:
param_bounds = ([Viso_min, Ac_min, As_min],
                [Viso_max, Ac_max, As_max])
In [14]:
popt, pcov = curve_fit(aniso_fit, aniso_df['az'], aniso_df['vel'], p0=(Viso_init, Ac_init, As_init), bounds=param_bounds)     #modified by an 20200421
In [15]:
#print (popt)
#print (pcov)
#Print results

Viso = popt[0]
Ac = popt[1]
As = popt[2]
print("# Viso =  %.5f +/- %.5f" % (Viso, math.sqrt(pcov[0, 0])))
print("# Ac = %.5f +/- %.5f" % (Ac, math.sqrt(pcov[1, 1])))
print("# As =  %.5f +/- %.5f" % (As, math.sqrt(pcov[2, 2])))

Viso_out = "{:.5f}".format(Viso)

# Viso =  198.81137 +/- 0.03245
# Ac = 3.15560 +/- 0.04590
# As =  -0.27780 +/- 0.04590
In [16]:
az_syn = np.arange(0,180,daz)
vel_syn = np.zeros(len(az_syn))
for i in range(len(az_syn)):
    vel_syn[i] = aniso_fit(az_syn[i], Viso, Ac, As)

resi_fit = np.zeros(len(az_syn))
for i in range(len(az_syn)):
    resi_fit[i]=vel_syn[i]-aniso_df['vel'][i]
    #print("# resi_fit = ", resi_fit[i])
    
resi_sq= np.sum(resi_fit**2)
#print(resi_sq)
data_sq= np.sum(aniso_df['vel']**2)
#print(data_sq)

VR = (1 - (resi_sq/data_sq)) * 100
#print("# VR (%) =  ", VR)

VR_out = "{:.5f}".format(VR)
print("# VR (%) =  ", VR_out)


#Root-mean-square deviation
rms_coeff = math.sqrt(resi_sq/len(az_syn))
rms_coeff_out = "{:.5f}".format(rms_coeff)

print("# rms_coeff = ", rms_coeff_out)

# VR (%) =   99.99991
# rms_coeff =  0.18643
In [17]:
x = aniso_df['vel']
y = vel_syn 
slope, intercept, r_value, p_value, std_err = sp.stats.linregress(x, y)
In [18]:
print("# r_value = ", r_value)

# r_value =  0.996554483065
In [19]:
R_squared = r_value**2
R_squared_out = "{:.5f}".format(R_squared)

print("# R_squared = ", R_squared_out)

# R_squared =  0.99312
In [20]:
syn_data = pd.DataFrame({ 'az' : az_syn, 'vel' : vel_syn })
#print(syn_data)
In [21]:
syn_output_fi = "fitting.out"
syn_data.to_csv(syn_output_fi, columns=['az','vel'], header=False, index=False, sep=' ')
In [22]:
#daz2 = 1.0 # every one degree
daz2 = 0.01 # every one degree

az_syn2 = np.arange(0,180,daz2)
vel_syn2 = np.zeros(len(az_syn2))
for i in range(len(az_syn2)):
    vel_syn2[i] = aniso_fit(az_syn2[i], Viso, Ac, As)
In [23]:
v_slow =  np.amin(vel_syn2)
v_fast =  np.amax(vel_syn2)

v_slow_index = np.where(  (vel_syn2 <= v_slow ) )
v_fast_index = np.where(  (v_fast <= vel_syn2) )

print("# v_slow_az_index = ", v_slow_index )
print("# v_fast_az_index = ", v_fast_index )

az_slow = az_syn2[v_slow_index]
az_fast = az_syn2[v_fast_index]


print("# az_slow = ", az_slow[0])
print("# az_fast = ", az_fast[0])


aniso_strength = ((v_fast - v_slow)/v_fast) * 100




az_fast_out = "{:.5f}".format(az_fast[0])

v_slow_out = "{:.5f}".format(v_slow)
v_fast_out = "{:.5f}".format(v_fast)
aniso_strength_out = "{:.5f}".format(aniso_strength)

print("# v_slow (m/s) = ", v_slow_out)
print("# v_fast (m/s) = ", v_fast_out)
print("# aniso_strength (%) =", aniso_strength_out)

# v_slow_az_index =  (array([8748]),)
# v_fast_az_index =  (array([17748]),)
# az_slow =  87.48
# az_fast =  177.48
# v_slow (m/s) =  195.64356
# v_fast (m/s) =  201.97917
# aniso_strength (%) = 3.13677
In [24]:
plotOPT = 0
plotOPT = 1
In [25]:
if plotOPT:
    #plt.rcParams['figure.figsize'] = 14, 5
    png_fi= "aniso_fit.png"
    fig = plt.figure(figsize=(14, 5))

    x_pos = 0.68
    y_pos = 0.80
    dy = 0.055
    t=fig.text(x_pos, y_pos, "VR (%): "+str(VR_out), ha='left')
    t.set_bbox(dict(facecolor='white', alpha=0.7, edgecolor='white'))


    t=fig.text(x_pos, y_pos-dy, "Viso (m/s): "+str(Viso_out), ha='left')
    t.set_bbox(dict(facecolor='white', alpha=0.7, edgecolor='white'))


    t=fig.text(x_pos, y_pos-(dy*2), "aniso_strength (%): "+str(aniso_strength_out), ha='left')
    t.set_bbox(dict(facecolor='white', alpha=0.7, edgecolor='white'))


    t=fig.text(x_pos, y_pos-(dy*3), "az_fast (deg.): "+str(az_fast_out), ha='left')
    t.set_bbox(dict(facecolor='white', alpha=0.7, edgecolor='white'))

    t=fig.text(x_pos, y_pos-(dy*4), "rms_coeff: "+str(rms_coeff_out), ha='left')
    t.set_bbox(dict(facecolor='white', alpha=0.7, edgecolor='white'))

    t=fig.text(x_pos, y_pos-(dy*5), "R_squared : "+str(R_squared_out), ha='left')
    t.set_bbox(dict(facecolor='white', alpha=0.7, edgecolor='white'))
    

    plt.scatter(aniso_df['az'], aniso_df['vel'], label = "data")
    plt.plot(az_syn, vel_syn, label = 'fitting')

    plt.xlabel('Azimuth (deg.)') 
    plt.ylabel('Vs (m/s)') 

    plt.grid(True)
    plt.legend(loc='upper left')
    plt.savefig(png_fi)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: