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


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


def detections_to_event_markers(fn_detections):
    markers = []
Marius's avatar
Marius committed
15
16
17
18
    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
19
20
                i, t_d, t_t, apeak, latpeak, lonpeak, xpeak, ypeak, zpeak = \
                    data
Marius's avatar
Marius committed
21
22
23
24
25
26
27
28
                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)
29
30
31
32

    return markers


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

Sebastian Heimann's avatar
Sebastian Heimann committed
265
266
267
268
        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
269
            s.add_markers(s.detections)
270

271
272
273
274
275
        for _ifc in s.config.image_function_contributions:
            if isinstance(_ifc, ifc.ManualPickIFC):
                s.add_markers(
                        pmarker.load_markers(_ifc.picks_path))

276
277
    receivers = config.get_receivers()
    stations = set()
278
    lats, lons = geo.points_coords(receivers, system='latlon')
279
280
281
282
283
    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
284
285
    paths = config.expand_path(config.data_paths)
    paths.append(config.get_ifm_dir_path())
286

287
    p = pile.make_pile(paths=paths, fileformat='detect')
Marius Kriegerowski's avatar
Marius Kriegerowski committed
288
289
290
291
292
293
294
295

    meta = {'tabu': True}
    for tr in p.iter_traces(trace_selector=lambda x: x.station == 'SMAX'):
        if tr.meta:
            tr.meta.update(meta)
        else:
            tr.meta = meta

296
    snuffler.snuffle(p, stations=stations,
297
298
299
300
301
                     launch_hook=load_snuffling)


__all__ = [
    'snuffle']