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


12
13
14
15
16
kind_default = '1 (green)'


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

    return markers


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

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

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

69
        ''' % self.show_comparison.__doc__
Marius Kriegerowski's avatar
Marius Kriegerowski committed
70
71
        return s

72
73
74
    def __init__(self):
        Snuffling.__init__(self)
        self.config = None
75
        self.detections = []
76
77
78
79
80
81

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

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

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

    def call(self):
114
        self.mycleanup()
115
        self.detections = []
116
        i_detection = 0
117
118
119
        zpeak = 0.
        lat = 0.
        lon = 0.
120
        for traces in self.chopper_selected_traces(
121
                mode='all',
122
123
                trace_selector=lambda x: x.station == "SMAX",
                fallback=True):
Marius Kriegerowski's avatar
wip py3    
Marius Kriegerowski committed
124
125
            tr_smax = [tr for tr in traces if tr.location == '']
            tr_i = [tr for tr in traces if tr.location == 'i']
126
127
128
            if not tr_i:
                tr_i = [None] * len(tr_smax)

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

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

146
147
148
                    e = model.Event(
                        time=t, name="%s-%s" % (i_detection, a), lat=lat,
                        lon=lon, depth=zpeak)
149
150
151
                    self.detections.append(
                        gui_util.EventMarker(
                            event=e, kind=int(self.marker_kind[0])))
152
                    i_detection += 1
153
        self.add_markers(self.detections)
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()
Marius Kriegerowski's avatar
wip py3    
Marius Kriegerowski committed
177
        return [x for x in markers if vtmin < x.tmin < vtmax]
178
179
180
181

    def show_comparison(self):
        '''
        Iterates through reference catalog and searches for lassie detection
182
        candidates in a time window of +- 1.5 seconds around the reference.
183
184
185
186
187
188
189

        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.
        '''
190
        scan_time = 3.
191
        # select_by = 'first'
192
193
194
195
196
197
198
199

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

208
            candidates = [detections[i] for i in i_want]
209
210
211
212
213
214
215
216
217
218
219
220
221
222

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

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

        ax = self.fig.add_subplot(111)
Marius Kriegerowski's avatar
wip py3    
Marius Kriegerowski committed
230
        compare_events = [x.get_event() for x in markers_compare]
231
232
233
234
235
236
        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)
237
        n_leftover_detections = len(detections)
238
239
240
241
        n_undetected = len(not_detected)

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


def __snufflings__():
    return [LassieSnuffling()]


def snuffle(config):
    global _lassie_config
    _lassie_config = copy.deepcopy(config)
Marius Kriegerowski's avatar
py3    
Marius Kriegerowski committed
259
260
    for _ifc in _lassie_config.image_function_contributions:
        _ifc.setup(config)
261
262
263
264
265

    def load_snuffling(win):
        s = LassieSnuffling()
        s.config = _lassie_config
        s.setup()
266
        win.pile_viewer.viewer.add_snuffling(s, reloaded=True)
267
        win.pile_viewer.viewer.add_blacklist_pattern('*.SMAX.i.*')
Marius Kriegerowski's avatar
Marius Kriegerowski committed
268
269
        for bl in _lassie_config.blacklist:
            win.pile_viewer.viewer.add_blacklist_pattern('%s.*' % bl)
270

Sebastian Heimann's avatar
Sebastian Heimann committed
271
272
273
274
        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
275
            s.add_markers(s.detections)
276

277
278
        for _ifc in s.config.image_function_contributions:
            if isinstance(_ifc, ifc.ManualPickIFC):
279
280
281
282
283
                markers_path_extra = _ifc.picks_path
            elif isinstance(_ifc, ifc.TemplateMatchingIFC):
                markers_path_extra = _ifc.template_markers_path

            s.add_markers(pmarker.load_markers(markers_path_extra))
284

285
286
    receivers = config.get_receivers()
    stations = set()
287
    lats, lons = geo.points_coords(receivers, system='latlon')
288
289
290
291
292
    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
293
294
    paths = config.expand_path(config.data_paths)
    paths.append(config.get_ifm_dir_path())
295

296
    p = pile.make_pile(paths=paths, fileformat='detect')
Marius Kriegerowski's avatar
Marius Kriegerowski committed
297
298
299
300
301
302
303
304

    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

305
    snuffler.snuffle(p, stations=stations,
306
307
308
309
310
                     launch_hook=load_snuffling)


__all__ = [
    'snuffle']