Skip to content
GitLab
Switch branch/tag
Find file Normal viewHistoryPermalink
snuffling.py 9.98 KB
Newer Older
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
1 
import copy
Marius Kriegerowski's avatar
import os    
Marius Kriegerowski committed
2 
import os
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
3
4
5 
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
6 
from lassie import geo
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
7
8
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
9
10
11
12
13 
kind_default = '1 (green)'


def detections_to_event_markers(fn_detections):
    markers = []
Marius's avatar
ignore missing detections.txt    
Marius committed
14
15
16
17
18
19
20
21
22
23
24
25
26 
    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
27
28
29
30 

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

    @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
43 
        <h2 style="test-align:center"> Scrutinize Lassie Performance and
Marius Kriegerowski's avatar
snuffling: doc string    
Marius Kriegerowski committed
44
45 
        Re-Detect</h2>
        <p>
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
46
47 
        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
48
49
50
51 
        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
52
53 
        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
54
55
56 
        (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
57
58
59
60
61 
        <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
62
63
64 
        </body>
        </html>
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
65 
        ''' % self.show_comparison.__doc__
Marius Kriegerowski's avatar
snuffling: doc string    
Marius Kriegerowski committed
66
67 
        return s
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
68
69
70
71
72
73
74
75
76 
    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
77
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
78 
        self.set_name('Lassie investigate')
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
79 
        self.add_parameter(Param('Tsearch', 'tsearch', 20., 0.01, 100))
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
80 
        self.add_parameter(Param(
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
81 
            'Detector threshold', 'detector_threshold', detector_default, 1.,
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
82 
            10000.))
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
83 
        self.add_parameter(Switch('Level Trace', 'level_trace', False))
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
84
85 
        self.add_parameter(Switch(
            'Hold comparison figure', 'hold_figure', False))
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
86 
        self.add_parameter(Choice(
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
87 
            'new marker kind', 'marker_kind', kind_default,
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
88
89 
            ['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
90
91
92
93
94
95
96 
        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
97 
        self.fframe = None
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
98 
        self.grid = self.config.get_grid()
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
99
100
101
102
103
104
105 

    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
106
107 

    def call(self):
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
108 
        self.mycleanup()
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
109 
        self.detections = []
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
110 
        i_detection = 0
Marius Kriegerowski's avatar
snufflings: backward compatibility    
Marius Kriegerowski committed
111
112
113 
        zpeak = 0.
        lat = 0.
        lon = 0.
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
114 
        for traces in self.chopper_selected_traces(
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
115 
                mode='all',
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
116
117 
                trace_selector=lambda x: x.station == "SMAX",
                fallback=True):
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
118 
            tr_smax = filter(lambda x: x.location == '', traces)
Marius Kriegerowski's avatar
snufflings: backward compatibility    
Marius Kriegerowski committed
119
120
121
122 
            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
123 
            for tr_i, tr_stackmax in zip(tr_i, tr_smax):
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
124 
                tpeaks, apeaks = tr_stackmax.peaks(
Marius Kriegerowski's avatar
snufflings: remove detector deadtime, flake8    
Marius Kriegerowski committed
125 
                    self.detector_threshold, self.tsearch)
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
126
127
128
129
130
131
132 
                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
133
134
135
136 
                    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
137
138
139 
                    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
140
141
142 
                    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
143 
                    i_detection += 1
Marius Kriegerowski's avatar
lassie detection snuffling  
Marius Kriegerowski committed
144 
        self.add_markers(self.detections)
Marius Kriegerowski's avatar
snufflings: hold figure after closing figure frame    
Marius Kriegerowski committed
145
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
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
172 
        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
173 
        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
174
175
176
177
178
179
180 

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

        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
191
192 
        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
193 
            i_want = num.where(num.abs(detection_times - mcompare.tmin)
Marius Kriegerowski's avatar
snuffling: fix comparison    
Marius Kriegerowski committed
194
195
196 
                               < (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
197
198 
                continue
Marius Kriegerowski's avatar
snuffling: fix comparison    
Marius Kriegerowski committed
199 
            candidates = [detections[i] for i in i_want]
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
200
201
202
203
204
205
206
207
208
209
210
211
212
213 

            # 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
214 
        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
215
216 
            self.fig.clf()
        else:
Marius Kriegerowski's avatar
snufflings: hold figure after closing figure frame    
Marius Kriegerowski committed
217
218 
            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
219
220
221
222
223
224
225
226
227 

        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
228 
        n_leftover_detections = len(detections)
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
229
230
231
232 
        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
233 
            (n_leftover_detections, n_undetected,
Marius Kriegerowski's avatar
snuffling: compare detections agains reference catalog    
Marius Kriegerowski committed
234
235
236
237
238
239
240 
             (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
241
242
243
244
245
246
247
248
249
250
251
252
253
254 


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
255 
        win.pile_viewer.viewer.add_snuffling(s, reloaded=True)
Marius Kriegerowski's avatar
snuffling: read and evaluate argmax indices trace    
Marius Kriegerowski committed
256 
        win.pile_viewer.viewer.add_blacklist_pattern('*.SMAX.i.*')
Marius Kriegerowski's avatar
respect blacklist    
Marius Kriegerowski committed
257
258 
        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
259
Marius Kriegerowski's avatar
respect blacklist    
Marius Kriegerowski committed
260
261
262
263 
        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
264
265
266 

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


__all__ = [
    'snuffle']