-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathstatistics.py
131 lines (106 loc) · 4.62 KB
/
statistics.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
#!/usr/bin/env python
import numpy as np
from matplotlib import pyplot as plt
import plot
import obspy
import multiprocessing as mp
ppt_font = {'family' : 'sans-serif',
'style' : 'normal',
'variant' : 'normal',
'weight' : 'bold',
'size' : 'xx-large'}
paper_font = {'family' : 'serif',
'style' : 'normal',
'variant' : 'normal',
'weight' : 'medium',
'size' : 'large'}
def bootstrap_compute(st,**kwargs):
repeat = kwargs.get('repeat',100)
i_resamp = kwargs.get('resamp',len(st))
p_lim = kwargs.get('p_lim',(-1.5,1.5))
x_lim = kwargs.get('x_lim',(-10,160))
p_tick = kwargs.get('p_tick',-0.1)
vesp_row = kwargs.get('vesp_row',1)
def random_gen(i_resamp):
rand_list = []
for ii in range(i_resamp):
rand_list.append(np.random.randint(0,i_resamp))
return rand_list
def resample_stream(st,rand_list):
st_resamp = obspy.core.Stream()
for ii in rand_list:
st_resamp.append(st[ii])
return st_resamp
vesp_list = []
for ii in range(repeat):
rand_list = random_gen(len(st))
st_resamp = resample_stream(st,rand_list)
vesp = plot.vespagram(st_resamp,p_lim=p_lim,
x_lim=x_lim,p_tick=p_tick,plot=False)
vesp_list.append(vesp)
std_array = np.std(vesp_list,axis=0)
return std_array
def bootstrap_color(std_array,**kwargs):
window_tuple = kwargs.get('window_tuple',(-10,230))
window_slowness = kwargs.get('window_slowness',(-1.5,1.5))
slowness_tick = kwargs.get('slowness_tick',-0.1)
plot_line = kwargs.get('plot_line',False)
save = kwargs.get('save',False)
clim = kwargs.get('clim',(np.log10(std_array).min(),
np.log10(std_array).max()))
cmap = kwargs.get('cmap','gnuplot')
font = kwargs.get('font',paper_font)
plot = kwargs.get('plot',True)
fig, ax = plt.subplots(figsize=(15,5))
image = ax.imshow(np.log10(std_array),aspect='auto',interpolation='lanczos',
extent=[window_tuple[0],window_tuple[1],
window_slowness[0],window_slowness[1]],cmap=cmap,
vmin=clim[0],vmax=clim[1])
cbar = fig.colorbar(image,ax=ax)
cbar.set_label('Log(Standard Deviation)',fontdict=font)
ax.set_xlim(window_tuple)
ax.xaxis.set(ticks=range(window_tuple[0],window_tuple[1],10))
ax.set_ylim(window_slowness)
ax.grid(color='w',lw=2,alpha=0.6)
ax.set_xlabel('Seconds after P',fontdict=font)
ax.set_ylabel('Slowness (s/deg)',fontdict=font)
plt.show()
def bootstrap_wave(vesp,std_array,**kwargs):
window_tuple = kwargs.get('window_tuple',(-10,230))
window_slowness = kwargs.get('window_slowness',(-1.5,1.5))
slowness_tick = kwargs.get('slowness_tick',-0.1)
plot_line = kwargs.get('plot_line',False)
save = kwargs.get('save',False)
clim = kwargs.get('clim',(np.log10(std_array).min(),
np.log10(std_array).max()))
cmap = kwargs.get('cmap','gnuplot')
font = kwargs.get('font',paper_font)
plot = kwargs.get('plot',True)
vesp_wave = vesp.copy()
fig, ax = plt.subplots(figsize=(15,5))
ax.set_xlim(window_tuple)
ax.set_ylim((window_slowness[0],window_slowness[1]+0.1))
ax.set_xlabel('Seconds after P',fontdict=font)
ax.set_ylabel('Slowness (s/deg)',fontdict=font)
ax.xaxis.set(ticks=range(window_tuple[0],window_tuple[1],10))
time_vec = np.linspace(window_tuple[0],window_tuple[1],
num=vesp_wave.shape[1])
for idx, ii in enumerate(np.arange(window_slowness[1],window_slowness[0],
slowness_tick)):
vesp_wave[idx,:] +=ii
std_upper = vesp_wave[idx,:]+std_array[idx,:]
std_lower = vesp_wave[idx,:]-std_array[idx,:]
ax.plot(time_vec,vesp_wave[idx,:],color='k',lw=0.5)
ax.fill_between(time_vec,std_lower,std_upper,alpha=0.5,color='r')
fig2, ax2 = plt.subplots(figsize=(15,5))
ax2.set_xlim(window_tuple)
ax2.set_ylim((window_slowness[0],window_slowness[1]+0.1))
ax2.set_xlabel('Seconds after P',fontdict=font)
ax2.set_ylabel('Stacked Amplitude',fontdict=font)
ax2.xaxis.set(ticks=range(window_tuple[0],window_tuple[1],10))
vesp_single = vesp_wave[15,:]/np.abs(vesp_wave[15.:]).max()
std_upper = vesp_single+std_array[15,:]/np.abs(vesp_wave[15,:]).max()
std_lower = vesp_single-std_array[15,:]/np.abs(vesp_wave[15,:]).max()
ax2.plot(time_vec,vesp_single,color='k',lw=0.5)
ax2.fill_between(time_vec,std_lower,std_upper,alpha=0.5,color='r')
plt.show()