-
Notifications
You must be signed in to change notification settings - Fork 1
/
pltspecYEAR.py
executable file
·124 lines (91 loc) · 3.51 KB
/
pltspecYEAR.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
#!/usr/bin/env python
import glob
import numpy as np
import sys
from obspy.signal.spectral_estimation import get_NLNM
import matplotlib.pyplot as plt
from multiprocessing import Pool
path = '/TEMP_TEST/PSDS/'
# Grab the NLNM and get power and per in the microseism band
minper = 5.
maxper = 10.
per, NLNM = get_NLNM()
micNLNM = NLNM[(minper <= per) & (per <= maxper)]
micper = per[(minper <= per) & (per <= maxper)]
def checkAlive(curfile):
# Boolean version of dead channel metric
try:
with open(curfile,'r') as f:
staper=[]
stapow=[]
for line in f:
line = (line.strip()).split(',')
if (1./float(line[1]) >= minper) and (1./float(line[1]) <= maxper):
staper.append(1./float(line[1]))
stapow.append(float(line[0]))
staper = np.asarray(staper)
stapow = np.asarray(stapow)
NLNMinterp = np.interp(1./staper,1./micper,micNLNM)
dbdiff = 1./float(len(NLNMinterp))*sum(stapow-NLNMinterp)
if dbdiff <= 0.:
result = False
else:
result = True
except:
result = False
return result
def getPercentile(snst):
year = snst.split('.')[3]
chan = snst.split('.')[2]
loc = snst.split('.')[1]
sta = snst.split('.')[0]
files = glob.glob(path + sta + '/' + year + '/*' + loc + '_' + chan + '*')
for idx, curfile in enumerate(files):
print 'On ' + sta + ' ' + loc + ' ' + chan + ' ' + year + ' ' + str(100*float(idx)/float(len(files))) + '%'
if checkAlive(curfile):
power = []
period =[]
with open(curfile,'r') as f:
for line in f:
power.append(float(line.split(",")[0]))
period.append(1./float(line.split(",")[1]))
if not 'powerArray' in vars():
powerArray = np.asarray(power)
else:
powerArray = np.vstack((powerArray,np.asarray(power)))
if 'powerArray' in vars():
# We Now have an array of all the powers and periods
(numfiles, length) = powerArray.shape
# Grab the 1st, 5th, 25th, and 50th
#percentiles = np.empty([4,length])
if '[12NE]' in chan:
chanName = chan[:2] + 'H'
else:
chanName = chan
fw = open('results/' + sta + loc + chanName + year + 'percentile','w')
fw.write('period, 1st, 5th, 25th, 50th \n')
for idx in range(0,length):
fw.write(str(period[idx]) + ', ' + str(np.percentile(powerArray[:,idx],1.)) + ', ' + \
str(np.percentile(powerArray[:,idx],5.)) + ', ' + str(np.percentile(powerArray[:,idx],25.)) + \
', ' + str(np.percentile(powerArray[:,idx],50.)) + '\n')
fw.close()
return
if __name__ == "__main__":
# Used for testing
#year = '1989'
#chan = 'LH[12NE]'
#sta = 'HRV'
#loc = ''
#getPercentile(sta + '.' + loc + '.' + chan + '.' + year)
#sys.exit()
for year in range(1989,2015):
stas = glob.glob(path + '*/' + str(year) + '/PSD*')
goodsncl = []
for val in stas:
chan = val.split('_')[5]
if chan[-1] in ['1','2','N','E']:
chan = chan[:-1] + '[12NE]'
goodsncl.append(val.split('_')[3] + '.' + val.split('_')[4] + '.' + chan + '.' + str(year))
goodsncl = list(set(goodsncl))
pool = Pool(6)
pool.map(getPercentile,goodsncl)