forked from ROBelgium/MSNoise
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy paths04stack.py
242 lines (205 loc) · 10.8 KB
/
s04stack.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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
"""MSNoise is capable of using a reference function defined by absolute or
relative dates span. For example, an absolute range could be "from 1 January
2010 to 31 December 2011" and a relative range could be "the last 200 days".
In the latter case, the REF will need to be exported at every run, meaning the
following steps (MWCS and DTT) will be executed on the whole configured period.
If the REF is defined between absolute dates, excluding "today", the MWCS and
DTT will only be calculated for new data (e.g. "yesterday" and "today").
The corresponding configuration bits are ``ref_begin`` and ``ref_end``. In the
future, we plan on allowing multiple references to be defined.
Only data for new/modified dates need to be exported. If any CC-job has been
marked "Done" within the last day, the stacks will be calculated and a new DTT
job will be inserted in the database. For dates in the period of interest, the
moving-window stack will only be exported if new/modified CCF is available.
The export directory are "REF/" and "DAY%03i/" where %03i will be replaced by
the number of days stacked together (DAYS_005 for a 5-days stack, e.g.).
Please note that within MSNoise, stacks are always *inclusive* of the time/day
mentionned. For example, a 5-days stack on January 10, will contain
cross-correlation functions computed for January 6, 7, 8, 9 AND 10!
The graphical representation centered on a "January 10" tick might then display
changes in the CCF that occurred *on* the 10th !
Moving-window stack length(s) are configured using the ``mov_stack`` bit.
Usage:
~~~~~~
The best way to call this code is to start it from the console (-h shows the
help)
.. code-block:: sh
$ s04stack.py -h
usage: s04stack.py [-h] [-r] [-m] [-i INTERVAL]
Compute [REF,MOV] stacks if jobs have been modified in the last i days.
optional arguments:
-h, --help show this help message and exit
-r, --ref Triggers the computation of REF stacks
-m, --mov Triggers the computation of MOV stacks
-i INTERVAL, --interval INTERVAL
Number of days before now to search for modified CC
jobs [default:1]
On a routine basis, one should thus run the following to compute REF *and* MOV
stacks:
.. code-block:: sh
python s04stack.py -r -m
While, when playing around with data, and surely on the first run, one should
define the *-i INTERVAL*, as jobs might have been marked "Done" more than 24
hours before running the stack. This, for example, will tell the code to search
for jobs marked in the last 10 days:
.. code-block:: sh
python s04stack.py -r -m -i 10
"""
import os
import argparse
import numpy as np
import scipy.signal
from scipy.stats.stats import nanmean
from database_tools import *
import logging
logging.basicConfig(level=logging.DEBUG,
filename="./stack.log",
format='%(asctime)s [%(levelname)s] %(message)s',
filemode='w')
console = logging.StreamHandler()
console.setLevel(logging.DEBUG)
formatter = logging.Formatter('%(asctime)s [%(levelname)s] %(message)s')
console.setFormatter(formatter)
logging.getLogger('').addHandler(console)
def stack(stype, interval=1):
"""Computes the REF/MOV stacks.
Parameters
----------
stype : {'mov', 'ref'}
Defines which of the REF or Moving-window stacks must be exported
interval : int, optional
Number of days before now to search for modified CC jobs
"""
db = connect()
components_to_compute = get_components_to_compute(db)
export_format = get_config(db, 'export_format')
if export_format == "BOTH":
mseed = True
sac = True
elif export_format == "SAC":
mseed = False
sac = True
elif export_format == "MSEED":
mseed = True
sac = False
maxlag = float(get_config(db, "maxlag"))
cc_sampling_rate = float(get_config(db, "cc_sampling_rate"))
if stype == "day":
start, end, datelist = build_daystack_datelist(db)
format = "stack"
elif stype == "mov":
start, end, datelist = build_movstack_datelist(db)
format = "matrix"
mov_stack = get_config(db, "mov_stack")
if mov_stack.count(',') == 0:
mov_stacks = [int(mov_stack), ]
else:
mov_stacks = [int(mi) for mi in mov_stack.split(',')]
if 1 in mov_stacks:
mov_stacks.remove(1) # remove 1 day stack, it should exist already
elif stype == "ref":
start, end, datelist = build_ref_datelist(db)
format = "stack"
for f in get_filters(db, all=False):
filterid = int(f.ref)
for components in components_to_compute:
for station1, station2 in get_station_pairs(db, used=True):
sta1 = "%s_%s" % (station1.net, station1.sta)
sta2 = "%s_%s" % (station2.net, station2.sta)
pair = "%s:%s" % (sta1, sta2)
logging.debug('Processing %s-%s-%i' %
(pair, components, filterid))
updated_days = updated_days_for_dates(db, start, end, pair.replace('_', '.'), type='CC', interval=datetime.timedelta(days=interval),returndays=True)
if len(updated_days) != 0:
logging.debug("New Data for %s-%s-%i" %
(pair, components, filterid))
nstack, stack_total = get_results(
db, sta1, sta2, filterid, components, datelist, format=format)
if nstack > 0:
if stype == "mov":
for i, date in enumerate(datelist):
jobadded = False
for mov_stack in mov_stacks:
if i < mov_stack:
low = 0
high = mov_stack
else:
low = i - mov_stack + 1
high = i + 1
newdata = False
for uday in datelist[low:high]:
if uday in updated_days:
newdata = True
break
if newdata:
corr = stack_total[low:high]
if not np.all(np.isnan(corr)):
day_name = "%s_%s" % (
sta1, sta2)
logging.debug("%s %s [%s - %s] (%i day stack)" % (
day_name, date, datelist[low], datelist[i], mov_stack))
corr = nanmean(corr, axis=0)
corr = scipy.signal.detrend(
corr)
stack_path = os.path.join(
"STACKS", "%02i" % filterid, "%03i_DAYS" % mov_stack, components, day_name)
filename = os.path.join(
stack_path, str(date))
if mseed:
export_mseed(
db, filename, pair, components, filterid, corr, maxlag=maxlag, cc_sampling_rate=cc_sampling_rate)
if sac:
export_sac(
db, filename, pair, components, filterid, corr, maxlag=maxlag, cc_sampling_rate=cc_sampling_rate)
day_name = "%s:%s" % (
sta1, sta2)
if not jobadded:
update_job(
db, date, day_name.replace('_', '.'), 'DTT', 'T')
jobadded = True
del corr
elif stype == "ref":
stack_path = os.path.join(
"STACKS", "%02i" % filterid, "REF", components)
ref_name = "%s_%s" % (sta1, sta2)
filename = os.path.join(stack_path, ref_name)
stack_total = scipy.signal.detrend(stack_total)
if mseed:
export_mseed(
db, filename, pair, components, filterid, stack_total)
if sac:
export_sac(
db, filename, pair, components, filterid, stack_total)
ref_name = "%s:%s" % (sta1, sta2)
update_job(
db, "REF", ref_name.replace('_', '.'), 'DTT', 'T')
del stack_total
def refstack(interval):
stack('ref', interval)
def movstack(interval):
stack('mov', interval)
if __name__ == "__main__":
parser = argparse.ArgumentParser(description='Compute [REF,MOV] stacks if\
jobs have been modified in the last i days.',
epilog=__doc__)
parser.add_argument('-r', '--ref', action="store_true",
help='Triggers the computation of REF stacks',
default=False)
parser.add_argument('-m', '--mov', action="store_true",
help='Triggers the computation of MOV stacks',
default=False)
parser.add_argument('-i', '--interval',
help='Number of days before now to search for\
modified CC jobs [default:1]', default=1, type=int)
args = parser.parse_args()
logging.info('Starting this program with: ref=%s, mov=%s, interval=%i'
% (args.ref, args.mov, args.interval))
db = connect()
if args.ref:
logging.info("*** Starting: REF Stack ***")
refstack(args.interval)
logging.info("*** Finished: REF Stack ***")
if args.mov:
logging.info("*** Starting: MOV Stack ***")
movstack(args.interval)
logging.info("*** Finished: MOV Stack ***")