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

10
11
logger = logging.getLogger('lassie.ifc')

Sebastian Heimann's avatar
Sebastian Heimann committed
12
13
14
guts_prefix = 'lassie'


Sebastian Heimann's avatar
Sebastian Heimann committed
15
16
17
18
19
20
21
22
23
24
25
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)


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

32
33
34
35
36
37
38
39
40
41
42
43
44
    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'
    )
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

    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
64
class IFC(common.HasPaths):
Sebastian Heimann's avatar
Sebastian Heimann committed
65
66
67
    '''Image function contribution.'''

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

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

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

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

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

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

93
94
95
    def get_table(self, grid, receivers):
        return self.shifter.get_table(grid, receivers)

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

99
100
101
102
103
104
    def get_fsmooth(self):
        return self.fmin

    def prescan(self, p):
        pass

105
106
107
108
    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
109
110
        pass

111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
    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
128
129
130
131
132
133
134
135
136
137
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.
    '''

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

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

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

150
151
152
153
154
    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
155

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

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

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

Sebastian Heimann's avatar
Sebastian Heimann committed
167
            tr.ydata = tr.ydata**2
168

Sebastian Heimann's avatar
Sebastian Heimann committed
169
170
            n = int(num.round(1./fsmooth / tr.deltat))
            taper = num.hanning(n)
171
            tr.set_ydata(fftconvolve(tr.get_ydata(), taper))
Sebastian Heimann's avatar
Sebastian Heimann committed
172
173
            tr.set_ydata(num.maximum(tr.ydata, 0.0))
            tr.shift(-(n/2.*tr.deltat))
Sebastian Heimann's avatar
Sebastian Heimann committed
174
            downsample(tr, deltat_cf)
Sebastian Heimann's avatar
Sebastian Heimann committed
175
176
177
178
179
180
181
182
183
184
            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()
185
                    sumtr.set_codes(channel='CF_%s' % self.name)
Sebastian Heimann's avatar
Sebastian Heimann committed
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
                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
213
    short_window = Float.T()
Sebastian Heimann's avatar
Sebastian Heimann committed
214
    window_ratio = Float.T(default=5.0)
Sebastian Heimann's avatar
Sebastian Heimann committed
215
216
217
218
219
220
221
    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
222
223

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

        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
234
235
            raise common.LassieError(
                'must use fnormalize or fnormalize_factor')
Sebastian Heimann's avatar
Sebastian Heimann committed
236
237
238
239
240

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

242
243
244
245
246
    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
247
        fnormalize = self.get_fnormalize()
Sebastian Heimann's avatar
Sebastian Heimann committed
248

249
250
251
        if self.trace_selector:
            trs = self.trace_selector(trs)

Sebastian Heimann's avatar
Sebastian Heimann committed
252
253
254
255
256
257
258
259
260
261
262
263
264
265
        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
266
        swin = self.short_window
Sebastian Heimann's avatar
Sebastian Heimann committed
267
268
269
270
271
272
273
274
        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()
275
                    sumtr.set_codes(channel='CF_%s' % self.name)
Sebastian Heimann's avatar
Sebastian Heimann committed
276
277
278
279
280
                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)
281

Sebastian Heimann's avatar
Sebastian Heimann committed
282
            sumtr.shift(-swin/2.)
283

Sebastian Heimann's avatar
Sebastian Heimann committed
284
            normtr = sumtr.copy()
Sebastian Heimann's avatar
Sebastian Heimann committed
285

Sebastian Heimann's avatar
Sebastian Heimann committed
286
287
288
289
            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)
290

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

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

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

Sebastian Heimann's avatar
Sebastian Heimann committed
302
            normtr.set_codes(channel='normtr')
303

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

Sebastian Heimann's avatar
Sebastian Heimann committed
309
            sumtr.set_ydata(sumtr.get_ydata() / normtr.get_ydata())
Sebastian Heimann's avatar
Sebastian Heimann committed
310

Sebastian Heimann's avatar
Sebastian Heimann committed
311
            downsample(sumtr, deltat_cf)
Sebastian Heimann's avatar
Sebastian Heimann committed
312

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

            dataset.append((nsl, sumtr))

        return dataset


320
321
class TemplateMatchingIFC(IFC):

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

325
    template_markers_path = common.Path.T(
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
        optional=True,
        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
348
349
350
351

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

352
353
354
355
356
357
358
359
360
361
362
363
364
        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
365
366
        markers = gui_util.load_markers(self.expand_path(
            self.template_markers_path))
367

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

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

        master_traces = []
        for marker in markers:
            if not marker.nslc_ids:
                trace_selector = trace_selector_global
            else:
Sebastian Heimann's avatar
flake8    
Sebastian Heimann committed
379
380
381
382
                def trace_selector(tr):
                    return (
                        marker.match_nslc(tr.nslc_id) and
                        trace_selector_global(tr))
383
384
385
386
387
388
389
390
391
392

            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()
Sebastian Heimann's avatar
Sebastian Heimann committed
393
            downsample(tr, 1./self.downsample_rate)
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411

            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)

412
413
414
415
    def deltat_cf_is_available(self, deltat_cf):
        return False

    def preprocess(self, trs, wmin, wmax, tpad_new, deltat_cf):
416
417
418
419
420
421
422
423
424
425
426
427
428

        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
            a = self.masters[nslc]
Sebastian Heimann's avatar
Sebastian Heimann committed
429
            downsample(b, 1./self.downsample_rate)
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463

            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


464
465
466
467
468
469
470
class ManualPickIFC(IFC):

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

    fsmooth = Float.T(default=0.1)
471
    picks_path = common.Path.T()
472
473
474
475
476
    picks_phasename = String.T()

    def __init__(self, *args, **kwargs):
        IFC.__init__(self, *args, **kwargs)

Sebastian Heimann's avatar
Sebastian Heimann committed
477
        markers = pmarker.load_markers(self.expand_path(self.picks_path))
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
        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

                ind = nsl_to_index[nsl]

                picked_index.append(ind)
                picked_times.append((marker.tmin + marker.tmax) * 0.5)

        self._nsl_to_index = nsl_to_index
        self._picked_index = num.array(picked_index, dtype=num.int64)
        self._picked_times = num.array(picked_times, dtype=num.float)

    def get_fsmooth(self):
        return self.fsmooth

503
504
505
506
507
    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()
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526

        if not trs:
            return []

        mask = num.logical_and(
            wmin <= self._picked_times,
            wmax >= self._picked_times)

        picked_times = self._picked_times[mask]
        picked_index = self._picked_index[mask]

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

            nsl = tr.nslc_id[:3]
            try:
                index = self._nsl_to_index[nsl]
                ts = picked_times[picked_index == index]
Sebastian Heimann's avatar
flake8    
Sebastian Heimann committed
527
                its = (num.round((ts - tr.tmin) / tr.deltat)).astype(num.int64)
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
                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()
546
                    sumtr.set_codes(channel='CF_%s' % self.name)
547
548
549
550
551
552
553
554
555
556
557
558
559
560
                    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
561
562
            downsample(sumtr, deltat_cf)

563
564
565
566
567
568
569
            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
570
571
572
573
__all__ = [
    'IFC',
    'WavePacketIFC',
    'OnsetIFC',
574
    'TemplateMatchingIFC',
Marius Kriegerowski's avatar
Marius Kriegerowski committed
575
    'TraceSelector',
Sebastian Heimann's avatar
Sebastian Heimann committed
576
]