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
        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
98
        self.detections = []
99
        self.fig = None
100
        self.fframe = None
101
        self.grid = self.config.get_grid()
102
103
104
105
106
107
108

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

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

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

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

143
144
145
                    e = model.Event(
                        time=t, name="%s-%s" % (i_detection, a), lat=lat,
                        lon=lon, depth=zpeak)
146
147
148
                    self.detections.append(
                        gui_util.EventMarker(
                            event=e, kind=int(self.marker_kind[0])))
149
                    i_detection += 1
150
        self.add_markers(self.detections)
151

152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
        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
174
        return [x for x in markers if vtmin < x.tmin < vtmax]
175
176
177
178

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

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

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

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

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

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

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

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


def __snufflings__():
    return [LassieSnuffling()]


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

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

Sebastian Heimann's avatar
Sebastian Heimann committed
268
269
270
271
        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
272
            s.add_markers(s.detections)
273

274
275
276
277
278
        for _ifc in s.config.image_function_contributions:
            if isinstance(_ifc, ifc.ManualPickIFC):
                s.add_markers(
                        pmarker.load_markers(_ifc.picks_path))

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

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

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

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


__all__ = [
    'snuffle']