snuffling: compare detections agains reference catalog

Comes handy to test lassie performance with synthetic datasets or to against another existing catalog

snuffling: compare detections agains reference catalog
Comes handy to test lassie performance with synthetic datasets or to against another existing catalog
bc956ec2 · Marius Kriegerowski · 4fce023d · bc956ec2
Commit bc956ec2 authored Nov 15, 2016 by Marius Kriegerowski
--- a/src/snuffling.py
+++ b/src/snuffling.py
@@ -4,6 +4,28 @@ from pyrocko import util, gui_util, model, orthodrome, pile, snuffler
 from pyrocko.snuffling import Snuffling, Param, Choice, Switch


+kind_default = '1 (green)'
+
+
+def detections_to_event_markers(fn_detections):
+    markers = []
+
+    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(ypeak), float(xpeak))
+            t = util.str_to_time("%s %s" % (t_d, t_t))
+            label = "%s-%s" % (i, apeak)
+            e = model.Event(lat=lat, lon=lon, depth=float(zpeak)*1000.,
+                            name=label, time=t)
+            m = gui_util.EventMarker(e, kind=int(kind_default[0]))
+            markers.append(m)
+
+    return markers
+
+
 class LassieSnuffling(Snuffling):

    @property
@@ -16,23 +38,29 @@ class LassieSnuffling(Snuffling):
        </style>
        </head>
        <body>
-        <h2 style="test-align:center"> Scrutinize Lassie Detections and
+        <h2 style="test-align:center"> Scrutinize Lassie Performance and
        Re-Detect</h2>
        <p>
-        Adjust the detector <i>threshold</i>, press <i>run</i>. From every 
-        instant, where the signal rises above <i>threshold</i>, a time length 
+        Adjust the detector <i>threshold</i>, press <i>run</i>. From every
+        instant, where the signal rises above <i>threshold</i>, a time length
        of <i>tsearch</i> seconds is searched for a maximum. Detections are
        added as event markers to the viewer.  </p>

        <p>
-        If you want to save the updated detections, it might be helpful to use the marker table 
+        If you want to save the updated detections, it might be helpful to use
+        the marker table
        (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>
+        <h3 style="test-align:center"> Compare Lassie Detections with Reference
+        Catalog</h3>
+        <p>
+        %s
+        </p>
        </body>
        </html>

-        '''
+        ''' % self.show_comparison.__doc__
        return s

    def __init__(self):
@@ -44,26 +72,49 @@ class LassieSnuffling(Snuffling):
            detector_default = self.config.detector_threshold
        else:
            detector_default = 100.
+
        self.set_name('Lassie investigate')
        self.add_parameter(Param('Tsearch', 'tsearch', 1., 0.01, 100))
        self.add_parameter(Param(
-            'Threshold', 'detector_threshold', detector_default, 100., 10000.))
+            'Detector threshold', 'detector_threshold', detector_default, 100.,
+            10000.))
+        self.add_parameter(Param(
+            'Detector deadtime [s]', 'detector_deadtime', 0, 0., 20.,
+            low_is_none=True))
        self.add_parameter(Switch('Level Trace', 'level_trace', False))
+        self.add_parameter(Switch(
+            'Hold comparison figure', 'hold_figure', False))
        self.add_parameter(Choice(
-            'new marker kind', 'marker_kind', '1 (green)',
+            'new marker kind', 'marker_kind', kind_default,
            ['1 (green)', '2 (blue)', '3 (orange)', '4 (purple)', '5 (brown)',
             '0 (red)']))
-        self.set_live_update(False)
+        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
+
+    def mycleanup(self):
+        viewer = self.get_viewer()
+        viewer.release_data(self._tickets)
+        viewer.remove_markers(self.detections)
+        self._tickets = []
+        self._markers = []

    def call(self):
-        self.cleanup()
+        self.mycleanup()
        self.detections = []
+
        for traces in self.chopper_selected_traces(
                trace_selector=lambda x: x.station == "SMAX",
                fallback=True):
            for tr_stackmax in traces:
                tpeaks, apeaks = tr_stackmax.peaks(
-                    self.detector_threshold, self.tsearch)
+                    self.detector_threshold, self.tsearch,
+                    deadtime=self.detector_deadtime)
+
                if self.level_trace:
                    ltrace = tr_stackmax.copy(data=False)
                    ltrace.set_ydata(
@@ -77,6 +128,101 @@ class LassieSnuffling(Snuffling):
                            event=e, kind=int(self.marker_kind[0])))

        self.add_markers(self.detections)
+        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
+        candidates in a time window of +- 2 seconds around the reference.
+
+        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.
+        '''
+        scan_time = 2.
+        select_by = 'first'
+
+        if not self.markers_compare:
+            self.fail('No catalog to compare to')
+
+        markers_compare = self.filter_visible(self.markers_compare)
+        not_detected = []
+        detected_compare = []
+        detections_success = []
+        detections = copy.deepcopy(self.filter_visible(self.detections))
+        for i_m, m in enumerate(markers_compare):
+            marker_times = num.array([m.tmin for m in detections])
+            i_want = num.where(num.abs(marker_times-m.tmin) < scan_time / 2.)
+            if len(i_want[0]) == 0:
+                not_detected.append(m)
+                continue
+
+            candidates = [detections[i] for i in i_want[0]]
+
+            # 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))
+
+            detected_compare.append(m)
+            detections_success.append((matched_marker, i_m))
+
+            for c in candidates:
+                detections.remove(c)
+
+        if self.hold_figure and self.fig:
+            self.fig.clf()
+        else:
+            self.fig = self.figure()
+
+        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)
+        n_detected = len(detections)
+        n_undetected = len(not_detected)
+
+        ax.text(
+            0.05, 0.95, 'Other detections: %s\nNot detected: %s (%1.1f %%)' %
+            (n_detected, n_undetected,
+             (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()


 def __snufflings__():
@@ -92,6 +238,9 @@ def snuffle(config):
        s.config = _lassie_config
        s.setup()
        win.pile_viewer.viewer.add_snuffling(s, reloaded=True)
+        s.detections = detections_to_event_markers(
+            _lassie_config.detections_path)
+        s.add_markers(s.detections)

    receivers = config.get_receivers()
    stations = set()
@@ -99,25 +248,10 @@ def snuffle(config):
        n, s, l = r.codes[:3]
        stations.add(model.Station(lat=r.lat, lon=r.lon,
                                   network=n, station=s, location=l))
-
-    markers = []
-    with open(config.detections_path, '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(ypeak), float(xpeak))
-            t = util.str_to_time("%s %s" % (t_d, t_t))
-            label = "%s-%s" % (i, apeak)
-            e = model.Event(lat=lat, lon=lon, depth=float(zpeak)*1000.,
-                            name=label, time=t)
-            m = gui_util.EventMarker(e)
-            markers.append(m)
-
    paths = config.data_paths
    paths.append('stackmax')
    p = pile.make_pile(paths=paths, fileformat='detect')
-    snuffler.snuffle(p, markers=markers, stations=stations,
+    snuffler.snuffle(p, stations=stations,
                     launch_hook=load_snuffling)