Skip to content
GitLab
Switch branch/tag
Find file Normal viewHistoryPermalink
snuffling.py 9.97 KB
Newer Older
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
1
2
3
4 
import copy
import numpy as num
from pyrocko import util, gui_util, model, orthodrome, pile, snuffler
from pyrocko.snuffling import Snuffling, Param, Choice, Switch
Sebastian Heimann's avatar
fix snuffling to use new api call    
Sebastian Heimann committed
5 
from lassie import geo
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
6
7
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
8
9
10
11
12 
kind_default = '1 (green)'


def detections_to_event_markers(fn_detections):
    markers = []
Marius's avatar
ignore missing detections.txt    
Marius committed
13
14
15
16
17
18
19
20
21
22
23
24
25 
    if fn_detections:
        with open(fn_detections, 'r') as f:
            for line in f.readlines():
                data = line.split()
                i, t_d, t_t, apeak, latpeak, lonpeak, xpeak, ypeak, zpeak = data
                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)
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
26
27
28
29 

    return markers
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
30 
class LassieSnuffling(Snuffling):
Marius Kriegerowski's avatar
snuffling: doc string    
Marius Kriegerowski committed
31
32
33
34
35
36
37
38
39
40
41 

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

        <p>
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
51
52 
        If you want to save the updated detections, it might be helpful to use
        the marker table
Marius Kriegerowski's avatar
snuffling: doc string    
Marius Kriegerowski committed
53
54
55 
        (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>
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
56
57
58
59
60 
        <h3 style="test-align:center"> Compare Lassie Detections with Reference
        Catalog</h3>
        <p>
        %s
        </p>
Marius Kriegerowski's avatar
snuffling: doc string    
Marius Kriegerowski committed
61
62
63 
        </body>
        </html>
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
64 
        ''' % self.show_comparison.__doc__
Marius Kriegerowski's avatar
snuffling: doc string    
Marius Kriegerowski committed
65
66 
        return s
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
67
68
69
70
71
72
73
74
75 
    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.
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
76
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
77 
        self.set_name('Lassie investigate')
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
78 
        self.add_parameter(Param('Tsearch', 'tsearch', 20., 0.01, 100))
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
79 
        self.add_parameter(Param(
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
80 
            'Detector threshold', 'detector_threshold', detector_default, 1.,
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
81 
            10000.))
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
82 
        self.add_parameter(Switch('Level Trace', 'level_trace', False))
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
83
84 
        self.add_parameter(Switch(
            'Hold comparison figure', 'hold_figure', False))
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
85 
        self.add_parameter(Choice(
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
86 
            'new marker kind', 'marker_kind', kind_default,
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
87
88 
            ['1 (green)', '2 (blue)', '3 (orange)', '4 (purple)', '5 (brown)',
             '0 (red)']))
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
89
90
91
92
93
94
95 
        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
Marius Kriegerowski's avatar
snufflings: hold figure after closing figure frame    
Marius Kriegerowski committed
96 
        self.fframe = None
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
97 
        self.grid = self.config.get_grid()
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
98
99
100
101
102
103
104 

    def mycleanup(self):
        viewer = self.get_viewer()
        viewer.release_data(self._tickets)
        viewer.remove_markers(self.detections)
        self._tickets = []
        self._markers = []
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
105
106 

    def call(self):
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
107 
        self.mycleanup()
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
108 
        self.detections = []
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
109 
        i_detection = 0
Marius Kriegerowski's avatar
snufflings: backward compatibility    
Marius Kriegerowski committed
110
111
112 
        zpeak = 0.
        lat = 0.
        lon = 0.
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
113 
        for traces in self.chopper_selected_traces(
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
114 
                mode='all',
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
115
116 
                trace_selector=lambda x: x.station == "SMAX",
                fallback=True):
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
117 
            tr_smax = filter(lambda x: x.location == '', traces)
Marius Kriegerowski's avatar
snufflings: backward compatibility    
Marius Kriegerowski committed
118
119
120
121 
            tr_i = filter(lambda x: x.location == 'i', traces)
            if not tr_i:
                tr_i = [None] * len(tr_smax)
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
122 
            for tr_i, tr_stackmax in zip(tr_i, tr_smax):
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
123 
                tpeaks, apeaks = tr_stackmax.peaks(
Marius Kriegerowski's avatar
snufflings: remove detector deadtime, flake8    
Marius Kriegerowski committed
124 
                    self.detector_threshold, self.tsearch)
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
125
126
127
128
129
130
131 
                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):
Marius Kriegerowski's avatar
snufflings: backward compatibility    
Marius Kriegerowski committed
132
133
134
135 
                    if tr_i:
                        lat, lon, xpeak, ypeak, zpeak = self.grid.index_to_location(
                            tr_i(t)[1])
                        lat, lon = orthodrome.ne_to_latlon(lat, lon, xpeak, ypeak)
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
136
137
138 
                    e = model.Event(
                        time=t, name="%s-%s" % (i_detection, a), lat=lat,
                        lon=lon, depth=zpeak)
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
139
140
141 
                    self.detections.append(
                        gui_util.EventMarker(
                            event=e, kind=int(self.marker_kind[0])))
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
142 
                    i_detection += 1
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
143 
        self.add_markers(self.detections)
Marius Kriegerowski's avatar
snufflings: hold figure after closing figure frame    
Marius Kriegerowski committed
144
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171 
        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
Marius Kriegerowski's avatar
snuffling: fix comparison    
Marius Kriegerowski committed
172 
        candidates in a time window of +- 1.5 seconds around the reference.
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
173
174
175
176
177
178
179 

        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.
        '''
Marius Kriegerowski's avatar
snuffling: fix comparison    
Marius Kriegerowski committed
180 
        scan_time = 3.
Marius Kriegerowski's avatar
snufflings: remove detector deadtime, flake8    
Marius Kriegerowski committed
181 
        # select_by = 'first'
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
182
183
184
185
186
187
188
189 

        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))
Marius Kriegerowski's avatar
snuffling: fix comparison    
Marius Kriegerowski committed
190
191 
        for i_m, mcompare in enumerate(markers_compare):
            detection_times = num.array([d.tmin for d in detections])
Marius Kriegerowski's avatar
snufflings: remove detector deadtime, flake8    
Marius Kriegerowski committed
192 
            i_want = num.where(num.abs(detection_times - mcompare.tmin)
Marius Kriegerowski's avatar
snuffling: fix comparison    
Marius Kriegerowski committed
193
194
195 
                               < (scan_time / 2.))[0]
            if len(i_want) == 0:
                not_detected.append(mcompare)
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
196
197 
                continue
Marius Kriegerowski's avatar
snuffling: fix comparison    
Marius Kriegerowski committed
198 
            candidates = [detections[i] for i in i_want]
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
199
200
201
202
203
204
205
206
207
208
209
210
211
212 

            # 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)
Marius Kriegerowski's avatar
snufflings: hold figure after closing figure frame    
Marius Kriegerowski committed
213 
        if self.hold_figure and self.fframe and not self.fframe.closed:
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
214
215 
            self.fig.clf()
        else:
Marius Kriegerowski's avatar
snufflings: hold figure after closing figure frame    
Marius Kriegerowski committed
216
217 
            self.fframe = self.pylab('Lassie', get='figure_frame')
            self.fig = self.fframe.gcf()
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
218
219
220
221
222
223
224
225
226 

        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)
Marius Kriegerowski's avatar
snuffling: fix comparison    
Marius Kriegerowski committed
227 
        n_leftover_detections = len(detections)
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
228
229
230
231 
        n_undetected = len(not_detected)

        ax.text(
            0.05, 0.95, 'Other detections: %s\nNot detected: %s (%1.1f %%)' %
Marius Kriegerowski's avatar
snuffling: fix comparison    
Marius Kriegerowski committed
232 
            (n_leftover_detections, n_undetected,
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
233
234
235
236
237
238
239 
             (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()
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
240
241
242
243
244
245
246
247
248
249
250
251
252
253 


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()
Marius Kriegerowski's avatar
snuffling: set reloaded to force visibility at start    
Marius Kriegerowski committed
254 
        win.pile_viewer.viewer.add_snuffling(s, reloaded=True)
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
255 
        win.pile_viewer.viewer.add_blacklist_pattern('*.SMAX.i.*')
Marius Kriegerowski's avatar
respect blacklist    
Marius Kriegerowski committed
256
257 
        for bl in _lassie_config.blacklist:
            win.pile_viewer.viewer.add_blacklist_pattern('%s.*' % bl)
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
258
Marius Kriegerowski's avatar
respect blacklist    
Marius Kriegerowski committed
259
260
261
262 
        if os.path.exists(_lassie_config.detections_path):
            s.detections = detections_to_event_markers(
                _lassie_config.detections_path)
            s.add_markers(s.detections)
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
263
264
265 

    receivers = config.get_receivers()
    stations = set()
Sebastian Heimann's avatar
fix snuffling to use new api call    
Sebastian Heimann committed
266 
    lats, lons = geo.points_coords(receivers, system='latlon')
Marius Kriegerowski's avatar
snuffling: correct receiver locations    
Marius Kriegerowski committed
267
268
269
270
271 
    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))
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
272 
    paths = config.data_paths
Daniela Kühn's avatar
Merged stackmax_path and master    
Daniela Kühn committed
273
274
275 
    if config.stackmax_path:
        paths.append(config.stackmax_path)
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
276 
    p = pile.make_pile(paths=paths, fileformat='detect')
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
277 
    snuffler.snuffle(p, stations=stations,
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
278
279
280
281
282 
                     launch_hook=load_snuffling)


__all__ = [
    'snuffle']