snuffling.py 8.84 KB
Newer Older
1
2
3
4
5
6
import copy
import numpy as num
from pyrocko import util, gui_util, model, orthodrome, pile, snuffler
from pyrocko.snuffling import Snuffling, Param, Choice, Switch


7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
kind_default = '1 (green)'


def detections_to_event_markers(fn_detections):
    markers = []

    with open(fn_detections, 'r') as f:
        for line in f.readlines():
            data = line.split()
            i, t_d, t_t, apeak, latpeak, lonpeak, xpeak, ypeak, zpeak = data
            lat, lon = orthodrome.ne_to_latlon(
                float(latpeak), float(lonpeak), float(ypeak), float(xpeak))
            t = util.str_to_time("%s %s" % (t_d, t_t))
            label = "%s-%s" % (i, apeak)
            e = model.Event(lat=lat, lon=lon, depth=float(zpeak)*1000.,
                            name=label, time=t)
            m = gui_util.EventMarker(e, kind=int(kind_default[0]))
            markers.append(m)

    return markers


29
class LassieSnuffling(Snuffling):
Marius Kriegerowski's avatar
Marius Kriegerowski committed
30
31
32
33
34
35
36
37
38
39
40

    @property
    def __doc__(self):
        s = '''
        <html>
        <head>
        <style type="text/css">
            body { margin-left:10px };
        </style>
        </head>
        <body>
41
        <h2 style="test-align:center"> Scrutinize Lassie Performance and
Marius Kriegerowski's avatar
Marius Kriegerowski committed
42
43
        Re-Detect</h2>
        <p>
44
45
        Adjust the detector <i>threshold</i>, press <i>run</i>. From every
        instant, where the signal rises above <i>threshold</i>, a time length
Marius Kriegerowski's avatar
Marius Kriegerowski committed
46
47
48
49
        of <i>tsearch</i> seconds is searched for a maximum. Detections are
        added as event markers to the viewer.  </p>

        <p>
50
51
        If you want to save the updated detections, it might be helpful to use
        the marker table
Marius Kriegerowski's avatar
Marius Kriegerowski committed
52
53
54
        (see <a href="http://pyrocko.org/v0.3/snuffler_tutorial.html"> Snuffler
        Tutorial</a> at the bottom) to sort all markers by their kind.
        </p>
55
56
57
58
59
        <h3 style="test-align:center"> Compare Lassie Detections with Reference
        Catalog</h3>
        <p>
        %s
        </p>
Marius Kriegerowski's avatar
Marius Kriegerowski committed
60
61
62
        </body>
        </html>

63
        ''' % self.show_comparison.__doc__
Marius Kriegerowski's avatar
Marius Kriegerowski committed
64
65
        return s

66
67
68
69
70
71
72
73
74
    def __init__(self):
        Snuffling.__init__(self)
        self.config = None

    def setup(self):
        if self.config:
            detector_default = self.config.detector_threshold
        else:
            detector_default = 100.
75

76
        self.set_name('Lassie investigate')
Marius Kriegerowski's avatar
Marius Kriegerowski committed
77
        self.add_parameter(Param('Tsearch', 'tsearch', 1., 0.01, 100))
78
        self.add_parameter(Param(
79
80
81
82
83
            'Detector threshold', 'detector_threshold', detector_default, 100.,
            10000.))
        self.add_parameter(Param(
            'Detector deadtime [s]', 'detector_deadtime', 0, 0., 20.,
            low_is_none=True))
84
        self.add_parameter(Switch('Level Trace', 'level_trace', False))
85
86
        self.add_parameter(Switch(
            'Hold comparison figure', 'hold_figure', False))
87
        self.add_parameter(Choice(
88
            'new marker kind', 'marker_kind', kind_default,
89
90
            ['1 (green)', '2 (blue)', '3 (orange)', '4 (purple)', '5 (brown)',
             '0 (red)']))
91
92
93
94
95
96
97
98
99
100
101
102
103
104
        self.add_trigger('load reference', self.load_comparison)
        self.add_trigger('show comparison', self.show_comparison)
        self.add_trigger('remove comparison', self.remove_comparison)
        # self.add_trigger('read Lassie config', self.fail)
        self.set_live_update(True)
        self.markers_compare = []
        self.fig = None

    def mycleanup(self):
        viewer = self.get_viewer()
        viewer.release_data(self._tickets)
        viewer.remove_markers(self.detections)
        self._tickets = []
        self._markers = []
105
106

    def call(self):
107
        self.mycleanup()
108
        self.detections = []
109

110
111
112
113
114
        for traces in self.chopper_selected_traces(
                trace_selector=lambda x: x.station == "SMAX",
                fallback=True):
            for tr_stackmax in traces:
                tpeaks, apeaks = tr_stackmax.peaks(
115
116
117
                    self.detector_threshold, self.tsearch,
                    deadtime=self.detector_deadtime)

118
119
120
121
122
123
124
125
126
127
128
129
130
                if self.level_trace:
                    ltrace = tr_stackmax.copy(data=False)
                    ltrace.set_ydata(
                        num.ones(
                            tr_stackmax.data_len()) * self.detector_threshold)
                    self.add_trace(ltrace)
                for t, a in zip(tpeaks, apeaks):
                    e = model.Event(time=t, name=str(a))
                    self.detections.append(
                        gui_util.EventMarker(
                            event=e, kind=int(self.marker_kind[0])))

        self.add_markers(self.detections)
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
        if self.hold_figure:
            self.show_comparison()

    def load_comparison(self):
        '''
        For comparison in synthetic tests.
        '''
        fn = self.input_filename(caption='Select an event catalog')
        kind_compare = 4
        compare_events = model.load_events(fn)
        markers = [gui_util.EventMarker(event=e, kind=kind_compare) for e in
                   compare_events]

        self.markers_compare = markers
        self.add_markers(self.markers_compare)

    def remove_comparison(self):
        '''Remove comparison markers from viewer.'''
        self.get_viewer().remove_markers(self.markers_compare)

    def filter_visible(self, markers):
        vtmin, vtmax = self.get_viewer().get_time_range()
        return filter(lambda x: vtmin < x.tmin < vtmax, markers)

    def show_comparison(self):
        '''
        Iterates through reference catalog and searches for lassie detection
        candidates in a time window of +- 2 seconds around the reference.

        If multiple candidates are available selects the first as the matching
        lassie detection for this reference.

        This option requires the catalog to contain only pyrocko.model.Event
        instances.
        '''
        scan_time = 2.
        select_by = 'first'

        if not self.markers_compare:
            self.fail('No catalog to compare to')

        markers_compare = self.filter_visible(self.markers_compare)
        not_detected = []
        detected_compare = []
        detections_success = []
        detections = copy.deepcopy(self.filter_visible(self.detections))
        for i_m, m in enumerate(markers_compare):
            marker_times = num.array([m.tmin for m in detections])
            i_want = num.where(num.abs(marker_times-m.tmin) < scan_time / 2.)
            if len(i_want[0]) == 0:
                not_detected.append(m)
                continue

            candidates = [detections[i] for i in i_want[0]]

            # if select_by == 'first':
            matched_marker = min(
                candidates, key=lambda x: x.get_event().time)

            # elif select_by == 'strongest':
            #     matched_marker = max(
            #         candidates, key=lambda x: float(x.get_event().name))

            detected_compare.append(m)
            detections_success.append((matched_marker, i_m))

            for c in candidates:
                detections.remove(c)

        if self.hold_figure and self.fig:
            self.fig.clf()
        else:
            self.fig = self.figure()

        ax = self.fig.add_subplot(111)
        compare_events = map(lambda x: x.get_event(), markers_compare)
        associated_events = [compare_events[a[1]] for a in detections_success]
        magnitudes = [e.get_event().magnitude for e in markers_compare]
        detected_magnitudes = [e.magnitude for e in associated_events]
        bins = num.linspace(-1, max(magnitudes), 30)
        ax.hist([detected_magnitudes, magnitudes], bins,
                label=['Lassie', 'Reference'], alpha=0.7)
        n_detected = len(detections)
        n_undetected = len(not_detected)

        ax.text(
            0.05, 0.95, 'Other detections: %s\nNot detected: %s (%1.1f %%)' %
            (n_detected, n_undetected,
             (float(n_undetected)/len(markers_compare)*100.)),
            transform=ax.transAxes)

        ax.set_xlabel('Magnitude')
        ax.set_ylabel('N detections')
        ax.legend()
        self.fig.canvas.draw()
226
227
228
229
230
231
232
233
234
235
236
237
238
239


def __snufflings__():
    return [LassieSnuffling()]


def snuffle(config):
    global _lassie_config
    _lassie_config = copy.deepcopy(config)

    def load_snuffling(win):
        s = LassieSnuffling()
        s.config = _lassie_config
        s.setup()
240
        win.pile_viewer.viewer.add_snuffling(s, reloaded=True)
241
242
243
        s.detections = detections_to_event_markers(
            _lassie_config.detections_path)
        s.add_markers(s.detections)
244
245
246
247
248
249
250
251
252
253

    receivers = config.get_receivers()
    stations = set()
    for r in receivers:
        n, s, l = r.codes[:3]
        stations.add(model.Station(lat=r.lat, lon=r.lon,
                                   network=n, station=s, location=l))
    paths = config.data_paths
    paths.append('stackmax')
    p = pile.make_pile(paths=paths, fileformat='detect')
254
    snuffler.snuffle(p, stations=stations,
255
256
257
258
259
                     launch_hook=load_snuffling)


__all__ = [
    'snuffle']