snuffling.py 10 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
    if fn_detections:
        with open(fn_detections, 'r') as f:
            for line in f.readlines():
                data = line.split()
Sebastian Heimann's avatar
flake8    
Sebastian Heimann committed
18
19
                i, t_d, t_t, apeak, latpeak, lonpeak, xpeak, ypeak, zpeak = \
                    data
Marius's avatar
Marius committed
20
21
22
23
24
25
26
27
                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)
28
29
30
31

    return markers


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

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

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

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

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

79
        self.set_name('Lassie investigate')
80
        self.add_parameter(Param('Tsearch', 'tsearch', 20., 0.01, 100))
81
        self.add_parameter(Param(
82
            'Detector threshold', 'detector_threshold', detector_default, 1.,
83
            10000.))
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
        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
98
        self.fframe = None
99
        self.grid = self.config.get_grid()
100
101
102
103
104
105
106

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

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

124
            for tr_i, tr_stackmax in zip(tr_i, tr_smax):
125
                tpeaks, apeaks = tr_stackmax.peaks(
126
                    self.detector_threshold, self.tsearch)
127
128
129
130
131
132
133
                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):
134
                    if tr_i:
Sebastian Heimann's avatar
flake8    
Sebastian Heimann committed
135
136
137
138
139
140
                        lat, lon, xpeak, ypeak, zpeak = \
                            self.grid.index_to_location(tr_i(t)[1])

                        lat, lon = orthodrome.ne_to_latlon(
                            lat, lon, xpeak, ypeak)

141
142
143
                    e = model.Event(
                        time=t, name="%s-%s" % (i_detection, a), lat=lat,
                        lon=lon, depth=zpeak)
144
145
146
                    self.detections.append(
                        gui_util.EventMarker(
                            event=e, kind=int(self.marker_kind[0])))
147
                    i_detection += 1
148
        self.add_markers(self.detections)
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
        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
177
        candidates in a time window of +- 1.5 seconds around the reference.
178
179
180
181
182
183
184

        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.
        '''
185
        scan_time = 3.
186
        # select_by = 'first'
187
188
189
190
191
192
193
194

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

203
            candidates = [detections[i] for i in i_want]
204
205
206
207
208
209
210
211
212
213
214
215
216
217

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

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

        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)
232
        n_leftover_detections = len(detections)
233
234
235
236
        n_undetected = len(not_detected)

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


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

Sebastian Heimann's avatar
Sebastian Heimann committed
264
265
266
267
        detections_path = _lassie_config.get_detections_path()

        if os.path.exists(detections_path):
            s.detections = detections_to_event_markers(detections_path)
Marius Kriegerowski's avatar
Marius Kriegerowski committed
268
            s.add_markers(s.detections)
269
270
271

    receivers = config.get_receivers()
    stations = set()
272
    lats, lons = geo.points_coords(receivers, system='latlon')
273
274
275
276
277
    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))

Sebastian Heimann's avatar
Sebastian Heimann committed
278
279
    paths = config.expand_path(config.data_paths)
    paths.append(config.get_ifm_dir_path())
280

281
    p = pile.make_pile(paths=paths, fileformat='detect')
282
    snuffler.snuffle(p, stations=stations,
283
284
285
286
287
                     launch_hook=load_snuffling)


__all__ = [
    'snuffle']