array_processing/example.py at master · uafgeotools/array_processing · GitHub

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
#%% User-defined parameters

fromwaveform_collectionimportgather_waveforms
fromobspy.coreimportUTCDateTime
importnumpyasnp

# Data collection
SOURCE='IRIS'
NETWORK='IM'
STATION='I53H?'
LOCATION='*'
CHANNEL='BDF'
START=UTCDateTime('2018-12-19T01:45:00')
END=START+20*60

# Filtering
FMIN=0.1# [Hz]
FMAX=1.0# [Hz]

# Array processing
WINLEN=50# [s]
WINOVER=0.5

#%% Grab and filter waveforms

st=gather_waveforms(SOURCE, NETWORK, STATION, LOCATION, CHANNEL, START, END,
remove_response=True)

st.filter('bandpass', freqmin=FMIN, freqmax=FMAX, corners=2, zerophase=True)
st.taper(max_percentage=0.01)

#%% Array processing and plotting using least squares

fromarray_processing.tools.plottingimportarray_plot
fromlts_arrayimportltsva

latlist= [tr.stats.latitudefortrinst]
lonlist= [tr.stats.longitudefortrinst]

#%% Array processing. ALPHA = 1.0: least squares processing.
ALPHA=1.0
vel, baz, t, mdccm, stdict, sigma_tau, conf_int_vel, conf_int_baz=ltsva(st, latlist, lonlist, WINLEN, WINOVER, ALPHA)

fig1, axs1=array_plot(st, t, mdccm, vel, baz, ccmplot=True, mcthresh=0.6,
sigma_tau=sigma_tau)

#%% Array processing. 0.5 <= ALPHA < 1.0: least trimmed squares processing.
ALPHA_LTS=0.50
vel_lts, baz_lts, t_lts, mdccm_lts, stdict_lts, sigma_tau, conf_int_vel, conf_int_baz=ltsva(st, latlist, lonlist, WINLEN, WINOVER, ALPHA_LTS)

fig1_lts, axs1_lts=array_plot(st, t_lts, mdccm_lts, vel_lts, baz_lts, ccmplot=True, mcthresh=0.6, sigma_tau=None, stdict=stdict_lts)

#%% Array uncertainty
fromarray_processing.algorithms.helpersimportgetrij
fromarray_processing.toolsimportarraySig
fromarray_processing.tools.plottingimportarraySigPlt, arraySigContourPlt

SIGLEVEL=1/st[0].stats.sampling_rate
KMAX=400
TRACE_VELOCITY=0.33
rij=getrij(latlist, lonlist)


sigV, sigTh, impResp, vel, th, kvec=arraySig(rij, kmax=KMAX,
sigLevel=SIGLEVEL)

fig2=arraySigPlt(rij, SIGLEVEL, sigV, sigTh, impResp, vel, th, kvec)

fig3=arraySigContourPlt(sigV, sigTh, vel, th, trace_v=TRACE_VELOCITY)

#%% Delay and sum beam

fromarray_processing.toolsimportbeamForm

data=np.array([tr.datafortrinst]).transpose()
beam=beamForm(data, rij, st[0].stats.sampling_rate, 50)

#%% Pure state filter

fromarray_processing.toolsimportpsf

x_psf, P=psf(data, p=2, w=3, n=3, window=None)