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


8
9
10
11
12
13
14
15
16
17
18
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(
19
                float(latpeak), float(lonpeak), float(xpeak), float(ypeak))
20
21
            t = util.str_to_time("%s %s" % (t_d, t_t))
            label = "%s-%s" % (i, apeak)
22
            e = model.Event(lat=lat, lon=lon, depth=float(zpeak),
23
24
25
26
27
28
29
                            name=label, time=t)
            m = gui_util.EventMarker(e, kind=int(kind_default[0]))
            markers.append(m)

    return markers


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

    @property
    def __doc__(self):
        s = '''
        <html>
        <head>
        <style type="text/css">
            body { margin-left:10px };
        </style>
        </head>
        <body>
42
        <h2 style="test-align:center"> Scrutinize Lassie Performance and
Marius Kriegerowski's avatar
Marius Kriegerowski committed
43
44
        Re-Detect</h2>
        <p>
45
46
        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
47
48
49
50
        of <i>tsearch</i> seconds is searched for a maximum. Detections are
        added as event markers to the viewer.  </p>

        <p>
51
52
        If you want to save the updated detections, it might be helpful to use
        the marker table
Marius Kriegerowski's avatar
Marius Kriegerowski committed
53
54
55
        (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>
56
57
58
59
60
        <h3 style="test-align:center"> Compare Lassie Detections with Reference
        Catalog</h3>
        <p>
        %s
        </p>
Marius Kriegerowski's avatar
Marius Kriegerowski committed
61
62
63
        </body>
        </html>

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

67
68
69
70
71
72
73
74
75
    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.
76

77
        self.set_name('Lassie investigate')
78
        self.add_parameter(Param('Tsearch', 'tsearch', 20., 0.01, 100))
79
        self.add_parameter(Param(
80
            'Detector threshold', 'detector_threshold', detector_default, 1.,
81
            10000.))
82
        self.add_parameter(Switch('Level Trace', 'level_trace', False))
83
84
        self.add_parameter(Switch(
            'Hold comparison figure', 'hold_figure', False))
85
        self.add_parameter(Choice(
86
            'new marker kind', 'marker_kind', kind_default,
87
88
            ['1 (green)', '2 (blue)', '3 (orange)', '4 (purple)', '5 (brown)',
             '0 (red)']))
89
90
91
92
93
94
95
        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
96
        self.fframe = None
97
        self.grid = self.config.get_grid()
98
99
100
101
102
103
104

    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
        i_detection = 0
110
111
112
        zpeak = 0.
        lat = 0.
        lon = 0.
113
        for traces in self.chopper_selected_traces(
114
                mode='all',
115
116
                trace_selector=lambda x: x.station == "SMAX",
                fallback=True):
117
            tr_smax = filter(lambda x: x.location == '', traces)
118
119
120
121
            tr_i = filter(lambda x: x.location == 'i', traces)
            if not tr_i:
                tr_i = [None] * len(tr_smax)

122
            for tr_i, tr_stackmax in zip(tr_i, tr_smax):
123
                tpeaks, apeaks = tr_stackmax.peaks(
124
                    self.detector_threshold, self.tsearch)
125
126
127
128
129
130
131
                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):
132
133
134
135
                    if tr_i:
                        lat, lon, xpeak, ypeak, zpeak = self.grid.index_to_location(
                            tr_i(t)[1])
                        lat, lon = orthodrome.ne_to_latlon(lat, lon, xpeak, ypeak)
136
137
138
                    e = model.Event(
                        time=t, name="%s-%s" % (i_detection, a), lat=lat,
                        lon=lon, depth=zpeak)
139
140
141
                    self.detections.append(
                        gui_util.EventMarker(
                            event=e, kind=int(self.marker_kind[0])))
142
                    i_detection += 1
143
        self.add_markers(self.detections)
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
        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
172
        candidates in a time window of +- 1.5 seconds around the reference.
173
174
175
176
177
178
179

        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.
        '''
180
        scan_time = 3.
181
        # select_by = 'first'
182
183
184
185
186
187
188
189

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

        markers_compare = self.filter_visible(self.markers_compare)
        not_detected = []
        detections_success = []
        detections = copy.deepcopy(self.filter_visible(self.detections))
190
191
        for i_m, mcompare in enumerate(markers_compare):
            detection_times = num.array([d.tmin for d in detections])
192
            i_want = num.where(num.abs(detection_times - mcompare.tmin)
193
194
195
                               < (scan_time / 2.))[0]
            if len(i_want) == 0:
                not_detected.append(mcompare)
196
197
                continue

198
            candidates = [detections[i] for i in i_want]
199
200
201
202
203
204
205
206
207
208
209
210
211
212

            # 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))

            detections_success.append((matched_marker, i_m))

            for c in candidates:
                detections.remove(c)

213
        if self.hold_figure and self.fframe and not self.fframe.closed:
214
215
            self.fig.clf()
        else:
216
217
            self.fframe = self.pylab('Lassie', get='figure_frame')
            self.fig = self.fframe.gcf()
218
219
220
221
222
223
224
225
226

        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)
227
        n_leftover_detections = len(detections)
228
229
230
231
        n_undetected = len(not_detected)

        ax.text(
            0.05, 0.95, 'Other detections: %s\nNot detected: %s (%1.1f %%)' %
232
            (n_leftover_detections, n_undetected,
233
234
235
236
237
238
239
             (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()
240
241
242
243
244
245
246
247
248
249
250
251
252
253


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()
254
        win.pile_viewer.viewer.add_snuffling(s, reloaded=True)
255
256
        win.pile_viewer.viewer.add_blacklist_pattern('*.SMAX.i.*')

257
258
259
        s.detections = detections_to_event_markers(
            _lassie_config.detections_path)
        s.add_markers(s.detections)
260
261
262

    receivers = config.get_receivers()
    stations = set()
263
264
265
266
267
268
    lats, lons = receiver.receivers_coords(receivers, system='latlon')
    for ir, (lat, lon) in enumerate(zip(lats, lons)):
        n, s, l = receivers[ir].codes[:3]
        stations.add(model.Station(
            lat=lat, lon=lon, network=n, station=s, location=l))

269
    paths = config.data_paths
270
271
272
    if config.stackmax_path:
        paths.append(config.stackmax_path)

273
    p = pile.make_pile(paths=paths, fileformat='detect')
274
    snuffler.snuffle(p, stations=stations,
275
276
277
278
279
                     launch_hook=load_snuffling)


__all__ = [
    'snuffle']