snuffling.py 11.1 KB
Newer Older
1
import copy
Marius Kriegerowski's avatar
Marius Kriegerowski committed
2
import os
3
import numpy as num
4
5
import logging

Marius Kriegerowski's avatar
Marius Kriegerowski committed
6
7
8
9
10
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
11
from lassie import geo, ifc
12
13


14
15
16
logger = logging.getLogger('lassie.snuffling')


17
18
19
20
21
kind_default = '1 (green)'


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

    return markers


40
class LassieSnuffling(Snuffling):
Marius Kriegerowski's avatar
Marius Kriegerowski committed
41
42
43
44
45
46
47
48
49
50
51

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

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

74
        ''' % self.show_comparison.__doc__
Marius Kriegerowski's avatar
Marius Kriegerowski committed
75
76
        return s

77
78
79
    def __init__(self):
        Snuffling.__init__(self)
        self.config = None
80
        self.detections = []
81
82
83
84
85
86

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

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

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

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

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

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

151
152
153
                    e = model.Event(
                        time=t, name="%s-%s" % (i_detection, a), lat=lat,
                        lon=lon, depth=zpeak)
154
155
156
                    self.detections.append(
                        gui_util.EventMarker(
                            event=e, kind=int(self.marker_kind[0])))
157
                    i_detection += 1
158
        self.add_markers(self.detections)
159

160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
        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
182
        return [x for x in markers if vtmin < x.tmin < vtmax]
183
184
185
186

    def show_comparison(self):
        '''
        Iterates through reference catalog and searches for lassie detection
187
        candidates in a time window of +- 1.5 seconds around the reference.
188
189
190
191
192
193
194

        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.
        '''
195
        scan_time = 3.
196
        # select_by = 'first'
197
198
199
200
201
202
203
204

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

213
            candidates = [detections[i] for i in i_want]
214
215
216
217
218
219
220
221
222
223
224
225
226
227

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

228
        if self.hold_figure and self.fframe and not self.fframe.closed:
229
230
            self.fig.clf()
        else:
231
232
            self.fframe = self.pylab('Lassie', get='figure_frame')
            self.fig = self.fframe.gcf()
233
234

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

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


def __snufflings__():
    return [LassieSnuffling()]


def snuffle(config):
    global _lassie_config
    _lassie_config = copy.deepcopy(config)
Marius Kriegerowski's avatar
py3    
Marius Kriegerowski committed
264
265
    for _ifc in _lassie_config.image_function_contributions:
        _ifc.setup(config)
266
267
268
269
270

    def load_snuffling(win):
        s = LassieSnuffling()
        s.config = _lassie_config
        s.setup()
271
        win.pile_viewer.viewer.add_snuffling(s, reloaded=True)
272
        win.pile_viewer.viewer.add_blacklist_pattern('*.SMAX.i.*')
Marius Kriegerowski's avatar
Marius Kriegerowski committed
273
274
        for bl in _lassie_config.blacklist:
            win.pile_viewer.viewer.add_blacklist_pattern('%s.*' % bl)
275

Sebastian Heimann's avatar
Sebastian Heimann committed
276
277
278
279
        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
280
            s.add_markers(s.detections)
281

282
283
        for _ifc in s.config.image_function_contributions:
            if isinstance(_ifc, ifc.ManualPickIFC):
284
285
286
                markers_path_extra = _ifc.picks_path
            elif isinstance(_ifc, ifc.TemplateMatchingIFC):
                markers_path_extra = _ifc.template_markers_path
287
288
            else:
                continue
289

290
291
292
293
294
            if os.path.exists(markers_path_extra):
                s.add_markers(pmarker.load_markers(markers_path_extra))
            else:
                logger.warn('No such file: %s (referenced in %s, named %s)' % (
                    markers_path_extra, _ifc.__class__.__name__, _ifc.name))
295

296
297
    receivers = config.get_receivers()
    stations = set()
298
    lats, lons = geo.points_coords(receivers, system='latlon')
299
300
301
302
303
    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
304
305
    paths = config.expand_path(config.data_paths)
    paths.append(config.get_ifm_dir_path())
306

307
    p = pile.make_pile(paths=paths, fileformat='detect')
Marius Kriegerowski's avatar
Marius Kriegerowski committed
308
309
310
311
312
313
314
315

    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

316
    snuffler.snuffle(p, stations=stations,
317
318
319
320
321
                     launch_hook=load_snuffling)


__all__ = [
    'snuffle']