ifc.py 16.9 KB
Newer Older
Marius Kriegerowski's avatar
wip py3    
Marius Kriegerowski committed
1
from __future__ import print_function
2
import logging
Sebastian Heimann's avatar
Sebastian Heimann committed
3
import numpy as num
Marius Kriegerowski's avatar
Marius Kriegerowski committed
4
from collections import defaultdict
5
from scipy.signal import fftconvolve
Marius Kriegerowski's avatar
Marius Kriegerowski committed
6
from pyrocko.guts import Object, String, Float, Bool, StringChoice, List, Dict
Marius Kriegerowski's avatar
Marius Kriegerowski committed
7
8
from pyrocko import trace, autopick, util, model
from pyrocko.gui import util as gui_util
9
from pyrocko import marker as pmarker
10
from lassie import shifter, common, geo
Sebastian Heimann's avatar
Sebastian Heimann committed
11

12
13
logger = logging.getLogger('lassie.ifc')

Sebastian Heimann's avatar
Sebastian Heimann committed
14
15
16
guts_prefix = 'lassie'


Sebastian Heimann's avatar
Sebastian Heimann committed
17
18
19
20
21
22
23
24
25
26
27
def downsample(tr, deltat):
    try:
        tr.downsample_to(
            deltat, demean=False, snap=True,
            allow_upsample_max=4)

    except util.UnavailableDecimation:
        logger.warn('using resample instead of decimation')
        tr.resample(deltat)


28
29
class TraceSelector(Object):
    '''
30
    Filter traces used in an IFC by NSLC-id lists and/or lists of regular
31
32
33
    expressions.
    '''

34
35
36
37
38
39
40
41
42
43
44
45
46
    white_list = List.T(
        optional=True,
        default=[],
        help='list of NSLC ids'
    )

    white_list_regex = List.T(
        String.T(
            default=[],
            optional=True,
        ),
        help='list of regular expressions'
    )
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

    def __call__(self, trs):
        matched = []
        for tr in trs:
            nslc = '.'.join(tr.nslc_id)

            if self.white_list and nslc in self.white_list:
                matched.append(tr)
                continue
            if self.white_list_regex and util.match_nslc(
                    self.white_list_regex, nslc):
                matched.append(tr)
        return matched

    def __str__(self):
        return '%s:\n wl: %s\n wlre: %s\n' % (
            self.__class__.__name__, self.white_list, self.white_list_regex)


Sebastian Heimann's avatar
Sebastian Heimann committed
66
class IFC(common.HasPaths):
Sebastian Heimann's avatar
Sebastian Heimann committed
67
68
69
    '''Image function contribution.'''

    name = String.T()
70
71
72
    weight = Float.T(
        default=1.0,
        help='global weight for this IFC')
73
74

    weights = Dict.T(
Marius Kriegerowski's avatar
Marius Kriegerowski committed
75
76
                String.T(help='NSL regular expression identifying stations'),
                Float.T(help='weighting factor'),
77
78
79
                optional=True,
                help='weight selected traces')

Sebastian Heimann's avatar
Sebastian Heimann committed
80
81
    fmin = Float.T()
    fmax = Float.T()
82
83
    shifter = shifter.Shifter.T(optional=True)

84
85
86
    trace_selector = TraceSelector.T(
        optional=True, help='select traces to be treated by this IFC')

87
    def __init__(self, *args, **kwargs):
88
        common.HasPaths.__init__(self, *args, **kwargs)
89
        self.shifter.t_tolerance = 1./(self.fmax * 2.)
Sebastian Heimann's avatar
Sebastian Heimann committed
90

Sebastian Heimann's avatar
Sebastian Heimann committed
91
92
93
94
    def setup(self, config):
        if self.shifter:
            self.shifter.setup(config)

95
96
97
    def get_table(self, grid, receivers):
        return self.shifter.get_table(grid, receivers)

Sebastian Heimann's avatar
Sebastian Heimann committed
98
99
100
    def get_tpad(self):
        return 4. / self.fmin

101
102
103
104
105
106
    def get_fsmooth(self):
        return self.fmin

    def prescan(self, p):
        pass

107
108
109
110
    def deltat_cf_is_available(self, deltat_cf):
        return False

    def preprocess(self, trs, wmin, wmax, tpad_new, deltat_cf):
Sebastian Heimann's avatar
Sebastian Heimann committed
111
112
        pass

113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
    def get_weights(self, nsls):
        if self.weights is None:
            return num.ones(len(nsls), dtype=num.float)

        else:
            weights = num.empty(len(nsls))
            selectors = self.weights.keys()
            for insl, nsl in enumerate(nsls):
                weights[insl] = 1.
                for selector in selectors:
                    if util.match_nslc(selector, nsl):
                        weights[insl] = self.weights[selector]
                        break

            return weights


Sebastian Heimann's avatar
Sebastian Heimann committed
130
131
132
133
134
135
136
137
138
139
class WavePacketIFC(IFC):

    '''
    An image function useful as a simple fast general purpose event detector.

    This detector type image function is sensible to transient wave packets
    traveling with a given velocity through the network, e.g. P phase together
    with P coda or S together with S coda or surface waves packets.
    '''

140
    fsmooth = Float.T(optional=True)
Sebastian Heimann's avatar
Sebastian Heimann committed
141
142
143
144
145
146
    fsmooth_factor = Float.T(default=0.1)

    def get_tpad(self):
        return 4. / self.get_fsmooth()

    def get_fsmooth(self):
147
148
149
150
        if self.fsmooth is not None:
            return self.fsmooth
        else:
            return self.fmin * self.fsmooth_factor
Sebastian Heimann's avatar
Sebastian Heimann committed
151

152
153
154
155
156
    def deltat_cf_is_available(self, deltat_cf):
        return deltat_cf < 0.025 / self.get_fsmooth()

    def preprocess(self, trs, wmin, wmax, tpad_new, deltat_cf):
        fsmooth = self.get_fsmooth()
Sebastian Heimann's avatar
Sebastian Heimann committed
157

Marius Kriegerowski's avatar
Marius Kriegerowski committed
158
159
160
        if self.trace_selector:
            trs = self.trace_selector(trs)

Sebastian Heimann's avatar
Sebastian Heimann committed
161
162
163
164
165
166
        if not trs:
            return []

        trs_filt = []
        for orig_tr in trs:
            tr = orig_tr.copy()
Marius Kriegerowski's avatar
Marius Kriegerowski committed
167
            tr.bandpass(4, self.fmin, self.fmax, demean=True)
168

Sebastian Heimann's avatar
Sebastian Heimann committed
169
            tr.ydata = tr.ydata**2
170

Sebastian Heimann's avatar
Sebastian Heimann committed
171
172
            n = int(num.round(1./fsmooth / tr.deltat))
            taper = num.hanning(n)
173
            tr.set_ydata(fftconvolve(tr.get_ydata(), taper))
Sebastian Heimann's avatar
Sebastian Heimann committed
174
175
            tr.set_ydata(num.maximum(tr.ydata, 0.0))
            tr.shift(-(n/2.*tr.deltat))
Sebastian Heimann's avatar
Sebastian Heimann committed
176
            downsample(tr, deltat_cf)
Sebastian Heimann's avatar
Sebastian Heimann committed
177
178
179
180
181
182
183
184
185
186
            trs_filt.append(tr)

        trs_by_nsl = common.grouped_by(trs_filt, lambda tr: tr.nslc_id[:3])

        dataset = []
        for nsl in sorted(trs_by_nsl.keys()):
            sumtr = None
            for tr in trs_by_nsl[nsl]:
                if sumtr is None:
                    sumtr = tr.copy()
187
                    sumtr.set_codes(channel='CF_%s' % self.name)
Sebastian Heimann's avatar
Sebastian Heimann committed
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
                else:
                    sumtr.add(tr)

            navg = 5./fsmooth / deltat_cf

            yavg = trace.moving_avg(sumtr.ydata, navg)

            if num.any(yavg == 0.0):
                continue

            sumtr.ydata /= yavg
            # sumtr.ydata -= 1.0
            sumtr.chop(wmin - tpad_new, wmax + tpad_new)

            dataset.append((nsl, sumtr))

        return dataset


class OnsetIFC(IFC):

    '''
    An image function sensible to sharp onsets in the signal.

    This image function is based on an STA/LTA.
    '''

Sebastian Heimann's avatar
Sebastian Heimann committed
215
    short_window = Float.T()
Sebastian Heimann's avatar
Sebastian Heimann committed
216
    window_ratio = Float.T(default=5.0)
Sebastian Heimann's avatar
Sebastian Heimann committed
217
218
219
220
221
222
223
    fsmooth = Float.T(optional=True)
    fsmooth_factor = Float.T(optional=True)
    fnormalize = Float.T(optional=True)
    fnormalize_factor = Float.T(optional=True)

    def get_tpad(self):
        return 2. / self.get_fnormalize()
Sebastian Heimann's avatar
Sebastian Heimann committed
224
225

    def get_fsmooth(self):
Sebastian Heimann's avatar
Sebastian Heimann committed
226
        if self.fsmooth is None and self.fsmooth_factor is None:
Sebastian Heimann's avatar
flake8    
Sebastian Heimann committed
227
            raise common.LassieError('must use fsmooth or fsmooth_factor')
Sebastian Heimann's avatar
Sebastian Heimann committed
228
229
230
231
232
233
234
235

        if self.fsmooth is not None:
            return self.fsmooth
        else:
            return self.fmin * self.fsmooth_factor

    def get_fnormalize(self):
        if self.fnormalize is None and self.fnormalize_factor is None:
Sebastian Heimann's avatar
flake8    
Sebastian Heimann committed
236
237
            raise common.LassieError(
                'must use fnormalize or fnormalize_factor')
Sebastian Heimann's avatar
Sebastian Heimann committed
238
239
240
241
242

        if self.fnormalize is not None:
            return self.fnormalize
        else:
            return self.fmin * self.fnormalize_factor
Sebastian Heimann's avatar
Sebastian Heimann committed
243

244
245
246
247
248
    def deltat_cf_is_available(self, deltat_cf):
        return deltat_cf < 0.025 / self.get_fsmooth()

    def preprocess(self, trs, wmin, wmax, tpad_new, deltat_cf):
        fsmooth = self.get_fsmooth()
Sebastian Heimann's avatar
Sebastian Heimann committed
249
        fnormalize = self.get_fnormalize()
Sebastian Heimann's avatar
Sebastian Heimann committed
250

251
252
253
        if self.trace_selector:
            trs = self.trace_selector(trs)

Sebastian Heimann's avatar
Sebastian Heimann committed
254
255
256
257
258
259
260
261
262
263
264
265
266
267
        if not trs:
            return []

        trs_filt = []
        for orig_tr in trs:
            tr = orig_tr.copy()
            tr.highpass(4, self.fmin, demean=True)
            tr.lowpass(4, self.fmax, demean=False)
            tr.ydata = tr.ydata**2

            trs_filt.append(tr)

        trs_by_nsl = common.grouped_by(trs_filt, lambda tr: tr.nslc_id[:3])

Sebastian Heimann's avatar
Sebastian Heimann committed
268
        swin = self.short_window
Sebastian Heimann's avatar
Sebastian Heimann committed
269
270
271
272
273
274
275
276
        lwin = self.window_ratio * swin

        dataset = []
        for nsl in sorted(trs_by_nsl.keys()):
            sumtr = None
            for tr in trs_by_nsl[nsl]:
                if sumtr is None:
                    sumtr = tr.copy()
277
                    sumtr.set_codes(channel='CF_%s' % self.name)
Sebastian Heimann's avatar
Sebastian Heimann committed
278
279
280
281
282
                else:
                    sumtr.add(tr)

            sumtr.set_ydata(sumtr.get_ydata().astype(num.float32))
            autopick.recursive_stalta(swin, lwin, 1.0, 4.0, 3.0, sumtr)
283

Sebastian Heimann's avatar
Sebastian Heimann committed
284
            sumtr.shift(-swin/2.)
285

Sebastian Heimann's avatar
Sebastian Heimann committed
286
            normtr = sumtr.copy()
Sebastian Heimann's avatar
Sebastian Heimann committed
287

Sebastian Heimann's avatar
Sebastian Heimann committed
288
289
290
291
            ntap_smooth = int(num.round(1./fsmooth / tr.deltat))
            ntap_normalize = int(num.round(1./fnormalize / tr.deltat))
            taper_smooth = num.hanning(ntap_smooth)
            taper_normalize = num.hanning(ntap_normalize)
292

Sebastian Heimann's avatar
Sebastian Heimann committed
293
294
            sumtr.shift(-(ntap_smooth / 2. * tr.deltat))
            normtr.shift(-(ntap_normalize / 2. * tr.deltat))
Sebastian Heimann's avatar
Sebastian Heimann committed
295
296

            sumtr.set_ydata(
Sebastian Heimann's avatar
Sebastian Heimann committed
297
                fftconvolve(sumtr.get_ydata()/len(trs_by_nsl), taper_smooth))
Sebastian Heimann's avatar
Sebastian Heimann committed
298

Sebastian Heimann's avatar
Sebastian Heimann committed
299
300
301
302
            normtr.set_ydata(
                fftconvolve(
                    normtr.get_ydata()/len(trs_by_nsl),
                    taper_normalize))
303

Sebastian Heimann's avatar
Sebastian Heimann committed
304
            normtr.set_codes(channel='normtr')
305

Sebastian Heimann's avatar
Sebastian Heimann committed
306
307
308
309
            tmin = max(sumtr.tmin, normtr.tmin)
            tmax = min(sumtr.tmax, normtr.tmax)
            sumtr.chop(tmin, tmax)
            normtr.chop(tmin, tmax)
310

Sebastian Heimann's avatar
Sebastian Heimann committed
311
            sumtr.set_ydata(sumtr.get_ydata() / normtr.get_ydata())
Sebastian Heimann's avatar
Sebastian Heimann committed
312

Sebastian Heimann's avatar
Sebastian Heimann committed
313
            downsample(sumtr, deltat_cf)
Sebastian Heimann's avatar
Sebastian Heimann committed
314

315
            sumtr.chop(wmin - tpad_new, wmax + tpad_new)
Sebastian Heimann's avatar
Sebastian Heimann committed
316
317
318
319
320
321

            dataset.append((nsl, sumtr))

        return dataset


322
323
class TemplateMatchingIFC(IFC):

324
    template_event_path = common.Path.T(
325
326
        help='Event parameters of the template')

327
    template_markers_path = common.Path.T(
328
        optional=False,
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
        help='File with markers defining the template')

    sum_square = Bool.T(
        default=False,
        help='Sum square of correlation')

    normalization = StringChoice.T(
        default='gliding',
        choices=['off', 'normal', 'gliding'])

    downsample_rate = Float.T(
        optional=True,
        help='If set, downsample to this sampling rate before processing [Hz]')

    def get_tpad(self):
        tmin_masters = min(tr.tmin for tr in self.masters.values())
        tmax_masters = max(tr.tmax for tr in self.masters.values())
        tmaster = tmax_masters - tmin_masters
        return tmaster

    def get_template_origin(self):
Sebastian Heimann's avatar
Sebastian Heimann committed
350
351
352
353

        event = model.load_one_event(self.expand_path(
            self.template_event_path))

354
355
356
357
358
359
360
361
362
363
364
365
366
        origin = geo.Point(
            lat=event.lat,
            lon=event.lon,
            z=event.depth)

        return origin

    def get_table(self, grid, receivers):
        origin = self.get_template_origin()
        return self.shifter.get_offset_table(grid, receivers, origin)

    def extract_template(self, p):

Sebastian Heimann's avatar
Sebastian Heimann committed
367
368
        markers = gui_util.load_markers(self.expand_path(
            self.template_markers_path))
369

Sebastian Heimann's avatar
flake8    
Sebastian Heimann committed
370
371
        def trace_selector_global(tr):
            return True
372
373
374
375
376
377

        period_highpass = 1./self.fmin
        tpad = 2 * period_highpass

        master_traces = []
        for marker in markers:
378
379
380
            if marker.tmin == marker.tmax:
                logger.warn('tmin == tmax in template marker %s' % marker)

381
382
383
            if not marker.nslc_ids:
                trace_selector = trace_selector_global
            else:
Sebastian Heimann's avatar
flake8    
Sebastian Heimann committed
384
385
386
387
                def trace_selector(tr):
                    return (
                        marker.match_nslc(tr.nslc_id) and
                        trace_selector_global(tr))
388
389
390
391
392
393
394
395
396
397

            master_traces.extend(p.all(
                tmin=marker.tmin,
                tmax=marker.tmax,
                trace_selector=trace_selector,
                tpad=tpad))

        masters = {}
        for xtr in master_traces:
            tr = xtr.copy()
398
399
            if self.downsample_rate is not None:
                downsample(tr, 1./self.downsample_rate)
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417

            tr.highpass(4, self.fmin, demean=False)
            tr.lowpass(4, self.fmax, demean=False)
            smin = round(xtr.wmin / tr.deltat) * tr.deltat
            smax = round(xtr.wmax / tr.deltat) * tr.deltat
            tr.chop(smin, smax)
            if tr.nslc_id in masters:
                raise common.LassieError(
                    'more than one waveform selected on trace with id "%s"'
                    % '.'.join(tr.nslc_id))

            masters[tr.nslc_id] = tr

        return masters

    def prescan(self, p):
        self.masters = self.extract_template(p)

418
419
420
421
    def deltat_cf_is_available(self, deltat_cf):
        return False

    def preprocess(self, trs, wmin, wmax, tpad_new, deltat_cf):
422
423
424
425
426
427
428
429
430
431
432
433

        tmin_masters = min(tr.tmin for tr in self.masters.values())
        tmax_masters = max(tr.tmax for tr in self.masters.values())
        tmaster = tmax_masters - tmin_masters
        tref = tmin_masters

        nsl_to_traces = defaultdict(list)
        for orig_b in trs:

            b = orig_b.copy()

            nslc = b.nslc_id
434
435
436
437
438
439
            a = self.masters.get(nslc, False)
            if not a:
                continue

            if self.downsample_rate is not None:
                downsample(b, 1./self.downsample_rate)
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473

            b.highpass(4, self.fmin, demean=False)
            b.lowpass(4, self.fmax, demean=False)
            smin = round((wmin - tmaster) / b.deltat) * b.deltat
            smax = round((wmax + tmaster) / b.deltat) * b.deltat
            b.chop(smin, smax)

            normalization = self.normalization
            if normalization == 'off':
                normalization = None

            c = trace.correlate(
                a, b, mode='valid', normalization=normalization)

            c.shift(-c.tmin + b.tmin - (a.tmin - tref))
            c.meta = {'tabu': True}
            if self.sum_square:
                c.ydata = c.ydata**2

            c.chop(wmin - tpad_new, wmax + tpad_new)

            nsl_to_traces[nslc[:3]].append(c)

        dataset = []
        for nsl, cs in nsl_to_traces.iteritems():
            csum = cs[0]
            for c in cs[1:]:
                csum.add(c)

            dataset.append((nsl, csum))

        return dataset


474
475
476
477
478
479
480
class ManualPickIFC(IFC):

    '''
    An image function based on manual picks.
    '''

    fsmooth = Float.T(default=0.1)
481
    picks_path = common.Path.T()
482
483
484
485
    picks_phasename = String.T()

    def __init__(self, *args, **kwargs):
        IFC.__init__(self, *args, **kwargs)
Sebastian Heimann's avatar
Sebastian Heimann committed
486
        self._picks_data = None
487

Sebastian Heimann's avatar
Sebastian Heimann committed
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
    def get_picks_data(self):
        if not self._picks_data:
            markers = pmarker.load_markers(self.expand_path(self.picks_path))
            nsl_to_index = {}
            picked_index = []
            picked_times = []
            index = -1
            for marker in markers:
                if isinstance(marker, pmarker.PhaseMarker) \
                        and marker.get_phasename() == self.picks_phasename:

                    nsl = marker.one_nslc()[:3]
                    if nsl not in nsl_to_index:
                        index += 1
                        nsl_to_index[nsl] = index
503

Sebastian Heimann's avatar
Sebastian Heimann committed
504
                    ind = nsl_to_index[nsl]
505

Sebastian Heimann's avatar
Sebastian Heimann committed
506
507
                    picked_index.append(ind)
                    picked_times.append((marker.tmin + marker.tmax) * 0.5)
508

Sebastian Heimann's avatar
Sebastian Heimann committed
509
510
511
512
            self._picks_data = (
                nsl_to_index,
                num.array(picked_index, dtype=num.int64),
                num.array(picked_times, dtype=num.float))
513

Sebastian Heimann's avatar
Sebastian Heimann committed
514
        return self._picks_data
515
516
517
518

    def get_fsmooth(self):
        return self.fsmooth

519
520
521
522
523
    def deltat_cf_is_available(self, deltat_cf):
        return deltat_cf < 0.025 / self.get_fsmooth()

    def preprocess(self, trs, wmin, wmax, tpad_new, deltat_cf):
        fsmooth = self.get_fsmooth()
524

Sebastian Heimann's avatar
Sebastian Heimann committed
525
526
        nsl_to_index, picked_index, picked_times = self.get_picks_data()

527
528
529
530
        if not trs:
            return []

        mask = num.logical_and(
Sebastian Heimann's avatar
Sebastian Heimann committed
531
532
            wmin <= picked_times,
            wmax >= picked_times)
533

Sebastian Heimann's avatar
Sebastian Heimann committed
534
535
        picked_times = picked_times[mask]
        picked_index = picked_index[mask]
536
537
538
539
540
541
542

        trs_filt = []
        for orig_tr in trs:
            tr = orig_tr.copy()

            nsl = tr.nslc_id[:3]
            try:
Sebastian Heimann's avatar
Sebastian Heimann committed
543
                index = nsl_to_index[nsl]
544
                ts = picked_times[picked_index == index]
Sebastian Heimann's avatar
flake8    
Sebastian Heimann committed
545
                its = (num.round((ts - tr.tmin) / tr.deltat)).astype(num.int64)
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
                its = its[num.logical_and(0 <= its, its < tr.data_len())]
                ydata = num.zeros(tr.data_len())
                ydata[its] = 1.0
                tr.set_ydata(ydata)
                trs_filt.append(tr)

            except KeyError:
                pass

        trs_by_nsl = common.grouped_by(trs_filt, lambda tr: tr.nslc_id[:3])

        dataset = []
        for nsl in sorted(trs_by_nsl.keys()):
            sumtr = None
            sumn = 0
            for tr in trs_by_nsl[nsl]:
                if sumtr is None:
                    sumtr = tr.copy()
564
                    sumtr.set_codes(channel='CF_%s' % self.name)
565
566
567
568
569
570
571
572
573
574
575
576
577
578
                    sumn = 1
                else:
                    sumtr.add(tr)
                    sumn += 1

            sumtr.ydata /= sumn

            ntap = int(num.round(1./fsmooth / tr.deltat))
            taper = num.hanning(ntap)
            sumtr.shift(-(ntap/2.*tr.deltat))

            sumtr.set_ydata(
                num.convolve(sumtr.get_ydata()/len(trs_by_nsl), taper))

Sebastian Heimann's avatar
Sebastian Heimann committed
579
580
            downsample(sumtr, deltat_cf)

581
582
583
584
585
586
587
            sumtr.chop(wmin - tpad_new, wmax + tpad_new)
            if num.any(sumtr != 0.):
                dataset.append((nsl, sumtr))

        return dataset


Sebastian Heimann's avatar
Sebastian Heimann committed
588
589
590
591
__all__ = [
    'IFC',
    'WavePacketIFC',
    'OnsetIFC',
592
    'TemplateMatchingIFC',
Marius Kriegerowski's avatar
Marius Kriegerowski committed
593
    'TraceSelector',
Sebastian Heimann's avatar
Sebastian Heimann committed
594
]