snuffling.py 9.98 KB
Newer Older
1
import copy
Marius Kriegerowski's avatar
Marius Kriegerowski committed
2
import os
3
4
5
import numpy as num
from pyrocko import util, gui_util, model, orthodrome, pile, snuffler
from pyrocko.snuffling import Snuffling, Param, Choice, Switch
6
from lassie import geo
7
8


9
10
11
12
13
kind_default = '1 (green)'


def detections_to_event_markers(fn_detections):
    markers = []
Marius's avatar
Marius committed
14
15
16
17
18
19
20
21
22
23
24
25
26
    if fn_detections:
        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(xpeak), float(ypeak))
                t = util.str_to_time("%s %s" % (t_d, t_t))
                label = "%s-%s" % (apeak, i)
                e = model.Event(lat=lat, lon=lon, depth=float(zpeak),
                                name=label, time=t)
                m = gui_util.EventMarker(e, kind=int(kind_default[0]))
                markers.append(m)
27
28
29
30

    return markers


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

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

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

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

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

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

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

    def call(self):
108
        self.mycleanup()
109
        self.detections = []
110
        i_detection = 0
111
112
113
        zpeak = 0.
        lat = 0.
        lon = 0.
114
        for traces in self.chopper_selected_traces(
115
                mode='all',
116
117
                trace_selector=lambda x: x.station == "SMAX",
                fallback=True):
118
            tr_smax = filter(lambda x: x.location == '', traces)
119
120
121
122
            tr_i = filter(lambda x: x.location == 'i', traces)
            if not tr_i:
                tr_i = [None] * len(tr_smax)

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

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

        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))
191
192
        for i_m, mcompare in enumerate(markers_compare):
            detection_times = num.array([d.tmin for d in detections])
193
            i_want = num.where(num.abs(detection_times - mcompare.tmin)
194
195
196
                               < (scan_time / 2.))[0]
            if len(i_want) == 0:
                not_detected.append(mcompare)
197
198
                continue

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

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

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

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

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


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()
255
        win.pile_viewer.viewer.add_snuffling(s, reloaded=True)
256
        win.pile_viewer.viewer.add_blacklist_pattern('*.SMAX.i.*')
Marius Kriegerowski's avatar
Marius Kriegerowski committed
257
258
        for bl in _lassie_config.blacklist:
            win.pile_viewer.viewer.add_blacklist_pattern('%s.*' % bl)
259

Marius Kriegerowski's avatar
Marius Kriegerowski committed
260
261
262
263
        if os.path.exists(_lassie_config.detections_path):
            s.detections = detections_to_event_markers(
                _lassie_config.detections_path)
            s.add_markers(s.detections)
264
265
266

    receivers = config.get_receivers()
    stations = set()
267
    lats, lons = geo.points_coords(receivers, system='latlon')
268
269
270
271
272
    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))

273
    paths = config.data_paths
274
275
276
    if config.stackmax_path:
        paths.append(config.stackmax_path)

277
    p = pile.make_pile(paths=paths, fileformat='detect')
278
    snuffler.snuffle(p, stations=stations,
279
280
281
282
283
                     launch_hook=load_snuffling)


__all__ = [
    'snuffle']