add ManualPickIFC and options to ease use of it

src/config.py

@@ -24,6 +24,14 @@ class Config(Object):
         String.T(),
         help='list of directories paths to search for data')
 
+    events_path = String.T(
+        optional=True,
+        help='limit processing to time windows around given events')
+
+    event_time_window_factor = Float.T(
+        default=2.,
+        help='controls length of time windows for event-wise processing')
+
     blacklist = List.T(
         String.T(),
         help='codes in the form NET.STA.LOC of receivers to be excluded')
@@ -70,6 +78,7 @@ class Config(Object):
         Object.__init__(self, *args, **kwargs)
         self._receivers = None
         self._grid = None
+        self._events = None
 
     def get_receivers(self):
         '''Aggregate receivers from different sources.'''
@@ -86,6 +95,15 @@ class Config(Object):
 
         return self._receivers
 
+    def get_events(self):
+        if self.events_path is None:
+            return None
+
+        if self._events is None:
+            self._events = model.load_events(self.events_path)
+
+        return self._events
+
     def get_grid(self):
         '''Get grid or make default grid.'''
 

src/core.py

-import sys
 import os
 import logging
 import os.path as op
@@ -28,7 +27,8 @@ def scan(
                 os.unlink(config.detections_path)
             else:
                 raise common.LassieError(
-                    'detections file already exists: %s' % config.detections_path)
+                    'detections file already exists: %s'
+                    % config.detections_path)
 
         util.ensuredirs(config.detections_path)
 
@@ -37,7 +37,6 @@ def scan(
 
     norm_map = gridmod.geometrical_normalization(grid, receivers)
 
-
     for data_path in config.data_paths:
         if not op.exists(data_path):
             raise common.LassieError(
@@ -79,154 +78,174 @@ def scan(
 
     tmin = override_tmin or config.tmin
     tmax = override_tmax or config.tmax
-    for trs in p.chopper(
-            tmin=tmin, tmax=tmax, tinc=tinc, tpad=tpad,
-            want_incomplete=False,
-            trace_selector=lambda tr: tr.nslc_id[:3] in station_index):
-
-        if not trs:
-            continue
-
-        logger.info('processing time window %s - %s' % (
-            util.time_to_str(trs[0].wmin),
-            util.time_to_str(trs[0].wmax)))
-
-        wmin = trs[0].wmin
-        wmax = trs[0].wmax
-
-        frames = None
-        pdata = []
-        trs_debug = []
-        shift_maxs = []
-        for iifc, ifc in enumerate(ifcs):
-            dataset = ifc.preprocess(trs, wmin, wmax, tshift_max)
-            if not dataset:
+
+    events = config.get_events()
+    twindows = []
+    if events is not None:
+        for ev in events:
+            if tmin <= ev.time <= tmax:
+                twindows.append((
+                    ev.time - tshift_max * config.event_time_window_factor,
+                    ev.time + tshift_max * config.event_time_window_factor))
+
+    else:
+        twindows.append((tmin, tmax))
+
+    for tmin_win, tmax_win in twindows:
+        for trs in p.chopper(
+                tmin=tmin_win, tmax=tmax_win, tinc=tinc, tpad=tpad,
+                want_incomplete=False,
+                trace_selector=lambda tr: tr.nslc_id[:3] in station_index):
+
+            if not trs:
                 continue
 
-            nstations_selected = len(dataset)
+            logger.info('processing time window %s - %s' % (
+                util.time_to_str(trs[0].wmin),
+                util.time_to_str(trs[0].wmax)))
 
-            nsls_selected, trs_selected = zip(*dataset)
+            wmin = trs[0].wmin
+            wmax = trs[0].wmax
 
-            trs_debug.extend(trs + list(trs_selected))
+            frames = None
+            pdata = []
+            trs_debug = []
+            shift_maxs = []
+            for iifc, ifc in enumerate(ifcs):
+                dataset = ifc.preprocess(trs, wmin, wmax, tshift_max)
+                if not dataset:
+                    continue
 
-            deltat_cf = trs_selected[0].deltat
+                nstations_selected = len(dataset)
 
-            t0 = (wmin / deltat_cf) * deltat_cf
+                nsls_selected, trs_selected = zip(*dataset)
 
-            istations_selected = num.array(
-                [station_index[nsl] for nsl in nsls_selected], dtype=num.int)
+                trs_debug.extend(trs + list(trs_selected))
 
-            arrays = [tr.ydata.astype(num.float) for tr in trs_selected]
+                deltat_cf = trs_selected[0].deltat
 
-            offsets = num.array(
-                [int(round((tr.tmin-t0) / deltat_cf)) for tr in trs_selected],
-                dtype=num.int32)
+                t0 = (wmin / deltat_cf) * deltat_cf
 
-            w = num.array(
-                [station_weights.get(nsl, 1.0) for nsl in nsls_selected],
-                dtype=num.float)
+                istations_selected = num.array(
+                    [station_index[nsl] for nsl in nsls_selected],
+                    dtype=num.int)
 
+                arrays = [tr.ydata.astype(num.float) for tr in trs_selected]
 
-            weights = num.ones((ngridpoints, nstations_selected))
-            weights *= w[num.newaxis, :]
-            weights *= ifc.weight
+                offsets = num.array(
+                    [int(round((tr.tmin-t0) / deltat_cf))
+                     for tr in trs_selected],
+                    dtype=num.int32)
 
-            shift_table = shift_tables[iifc]
+                w = num.array(
+                    [station_weights.get(nsl, 1.0) for nsl in nsls_selected],
+                    dtype=num.float)
 
-            shifts = -num.round(
-                shift_table[:, istations_selected] /
-                deltat_cf).astype(num.int32)
+                weights = num.ones((ngridpoints, nstations_selected))
+                weights *= w[num.newaxis, :]
+                weights *= ifc.weight
 
-            pdata.append((list(trs_selected), shift_table, ifc))
+                shift_table = shift_tables[iifc]
 
-            iwmin = int(round((wmin-t0) / deltat_cf))
-            iwmax = int(round((wmax-t0) / deltat_cf))
+                shifts = -num.round(
+                    shift_table[:, istations_selected] /
+                    deltat_cf).astype(num.int32)
 
-            lengthout = iwmax - iwmin + 1
+                pdata.append((list(trs_selected), shift_table, ifc))
 
-            shift_maxs.append(num.max(-shifts) * deltat_cf)
+                iwmin = int(round((wmin-t0) / deltat_cf))
+                iwmax = int(round((wmax-t0) / deltat_cf))
 
-            frames, ioff = parstack(
-                arrays, offsets, shifts, weights, 0,
-                offsetout=iwmin,
-                lengthout=lengthout,
-                result=frames,
-                impl='openmp')
+                lengthout = iwmax - iwmin + 1
 
-        shift_max = max(shift_maxs)
+                shift_maxs.append(num.max(-shifts) * deltat_cf)
 
-        if config.sharpness_normalization:
-            frame_maxs = frames.max(axis=0)
-            frame_means = num.abs(frames).mean(axis=0)
-            frames *= (frame_maxs / frame_means)[num.newaxis, :]
-            frames *= norm_map[:, num.newaxis]
+                frames, ioff = parstack(
+                    arrays, offsets, shifts, weights, 0,
+                    offsetout=iwmin,
+                    lengthout=lengthout,
+                    result=frames,
+                    impl='openmp')
 
-        frame_maxs = frames.max(axis=0)
+            shift_max = max(shift_maxs)
 
-        tmin_frames = t0 + ioff * deltat_cf
+            if config.sharpness_normalization:
+                frame_maxs = frames.max(axis=0)
+                frame_means = num.abs(frames).mean(axis=0)
+                frames *= (frame_maxs / frame_means)[num.newaxis, :]
+                frames *= norm_map[:, num.newaxis]
 
-        tr_stackmax = trace.Trace(
-            '', 'SMAX', '', '',
-            tmin=tmin_frames,
-            deltat=deltat_cf,
-            ydata=frame_maxs)
+            frame_maxs = frames.max(axis=0)
 
-        trs_debug.append(tr_stackmax)
+            tmin_frames = t0 + ioff * deltat_cf
 
-        # trace.snuffle(trs_debug)
+            tr_stackmax = trace.Trace(
+                '', 'SMAX', '', '',
+                tmin=tmin_frames,
+                deltat=deltat_cf,
+                ydata=frame_maxs)
 
-        ydata_window = tr_stackmax.chop(wmin, wmax, inplace=False).get_ydata()
+            trs_debug.append(tr_stackmax)
 
-        logger.info('CF stats: min %g, max %g, median %g' % (
-            num.min(ydata_window),
-            num.max(ydata_window),
-            num.median(ydata_window)))
+            # trace.snuffle(trs_debug)
 
-        tpeaks, apeaks = tr_stackmax.peaks(
-            config.detector_threshold, shift_max + 1.0/fsmooth_min)
+            ydata_window = tr_stackmax.chop(
+                wmin, wmax, inplace=False).get_ydata()
 
-        for (tpeak, apeak) in zip(tpeaks, apeaks):
-            if not (wmin <= tpeak and tpeak < wmax):
-                continue
+            logger.info('CF stats: min %g, max %g, median %g' % (
+                num.min(ydata_window),
+                num.max(ydata_window),
+                num.median(ydata_window)))
+
+            tpeaks, apeaks = tr_stackmax.peaks(
+                config.detector_threshold, shift_max + 1.0/fsmooth_min)
+
+            for (tpeak, apeak) in zip(tpeaks, apeaks):
+                if not (wmin <= tpeak and tpeak < wmax):
+                    continue
 
-            logger.info('detection: %s %g' % (
-                util.time_to_str(tpeak),
-                apeak))
+                logger.info('detection: %s %g' % (
+                    util.time_to_str(tpeak),
+                    apeak))
 
-            iframe = int(round(((tpeak-t0) - ioff*deltat_cf) / deltat_cf))
-            frame = frames[:, iframe]
-            imax = num.argmax(frame)
+                iframe = int(round(((tpeak-t0) - ioff*deltat_cf) / deltat_cf))
+                frame = frames[:, iframe]
+                imax = num.argmax(frame)
 
-            latpeak, lonpeak, xpeak, ypeak, zpeak = grid.index_to_location(imax)
+                latpeak, lonpeak, xpeak, ypeak, zpeak = \
+                    grid.index_to_location(imax)
 
-            idetection += 1
+                idetection += 1
 
-            if config.detections_path:
-                f = open(config.detections_path, 'a')
-                f.write('%06i %s %g %g %g %g %g %g\n' % (
-                    idetection,
-                    util.time_to_str(tpeak, format='%Y-%m-%d %H:%M:%S.6FRAC'),
-                    apeak,
-                    latpeak, lonpeak, xpeak, ypeak, zpeak))
+                if config.detections_path:
+                    f = open(config.detections_path, 'a')
+                    f.write('%06i %s %g %g %g %g %g %g\n' % (
+                        idetection,
+                        util.time_to_str(
+                            tpeak,
+                            format='%Y-%m-%d %H:%M:%S.6FRAC'),
+                        apeak,
+                        latpeak, lonpeak, xpeak, ypeak, zpeak))
 
-                f.close()
+                    f.close()
 
-            if show_detections:
-                fmin = min(ifc.fmin for ifc in ifcs)
-                fmax = min(ifc.fmax for ifc in ifcs)
-                plot.plot_detection(
-                    grid, receivers, frames, tmin_frames,
-                    deltat_cf, imax, iframe, fsmooth_min, xpeak, ypeak, zpeak,
-                    tr_stackmax, tpeaks, apeaks, config.detector_threshold,
-                    pdata, trs, fmin, fmax, idetection,
-                    grid_station_shift_max=shift_max,
-                    movie=show_movie)
+                if show_detections:
+                    fmin = min(ifc.fmin for ifc in ifcs)
+                    fmax = min(ifc.fmax for ifc in ifcs)
+                    plot.plot_detection(
+                        grid, receivers, frames, tmin_frames,
+                        deltat_cf, imax, iframe, fsmooth_min, xpeak, ypeak,
+                        zpeak,
+                        tr_stackmax, tpeaks, apeaks, config.detector_threshold,
+                        pdata, trs, fmin, fmax, idetection,
+                        grid_station_shift_max=shift_max,
+                        movie=show_movie)
 
-            if stop_after_first:
-                return
+                if stop_after_first:
+                    return
 
-        tr_stackmax.chop(wmin, wmax)
-        io.save([tr_stackmax], 'stackmax/trace_%(tmin)s.mseed')
+            tr_stackmax.chop(wmin, wmax)
+            io.save([tr_stackmax], 'stackmax/trace_%(tmin)s.mseed')
 
 
 __all__ = [

src/ifc.py

@@ -4,6 +4,7 @@ from collections import defaultdict
 import numpy as num
 from pyrocko.guts import Object, String, Float, Bool, StringChoice
 from pyrocko import trace, autopick, util, model, gui_util
+from pyrocko import marker as pmarker
 from lassie import shifter, common, geo
 
 logger = logging.getLogger('lassie.ifc')
@@ -170,7 +171,8 @@ class OnsetIFC(IFC):
 
             sumtr.set_ydata(sumtr.get_ydata().astype(num.float32))
             autopick.recursive_stalta(swin, lwin, 1.0, 4.0, 3.0, sumtr)
-            sumtr.shift(-(swin + n/2.*sumtr.deltat))
+            #sumtr.shift(-(swin + n/2.*sumtr.deltat))
+            sumtr.shift(-swin)
 
             ntap = int(num.round(1./fsmooth / tr.deltat))
             taper = num.hanning(ntap)
@@ -315,7 +317,7 @@ class TemplateMatchingIFC(IFC):
 
             c.shift(-c.tmin + b.tmin - (a.tmin - tref))
             c.meta = {'tabu': True}
-            #c.set_codes(location='cc')
+            # c.set_codes(location='cc')
             if self.sum_square:
                 c.ydata = c.ydata**2
 
@@ -334,6 +336,112 @@ class TemplateMatchingIFC(IFC):
         return dataset
 
 
+class ManualPickIFC(IFC):
+
+    '''
+    An image function based on manual picks.
+    '''
+
+    fsmooth = Float.T(default=0.1)
+    picks_path = String.T()
+    picks_phasename = String.T()
+
+    def __init__(self, *args, **kwargs):
+        IFC.__init__(self, *args, **kwargs)
+
+        markers = pmarker.load_markers(self.picks_path)
+        nsl_to_index = {}
+        picked_index = []
+        picked_times = []
+        index = -1
+        for marker in markers:
+            if isinstance(marker, pmarker.PhaseMarker) \
+                    and marker.get_phasename() == self.picks_phasename:
+
+                nsl = marker.one_nslc()[:3]
+                if nsl not in nsl_to_index:
+                    index += 1
+                    nsl_to_index[nsl] = index
+
+                ind = nsl_to_index[nsl]
+
+                picked_index.append(ind)
+                picked_times.append((marker.tmin + marker.tmax) * 0.5)
+
+        self._nsl_to_index = nsl_to_index
+        self._picked_index = num.array(picked_index, dtype=num.int64)
+        self._picked_times = num.array(picked_times, dtype=num.float)
+
+    def get_fsmooth(self):
+        return self.fsmooth
+
+    def preprocess(self, trs, wmin, wmax, tpad_new):
+        fsmooth = self.fsmooth
+
+        if not trs:
+            return []
+
+        deltat_cf = trs[0].deltat
+        while deltat_cf * 2 < 0.25 / self.fmin:
+            deltat_cf *= 2
+
+        mask = num.logical_and(
+            wmin <= self._picked_times,
+            wmax >= self._picked_times)
+
+        picked_times = self._picked_times[mask]
+        picked_index = self._picked_index[mask]
+
+        trs_filt = []
+        for orig_tr in trs:
+            tr = orig_tr.copy()
+
+            nsl = tr.nslc_id[:3]
+            try:
+                index = self._nsl_to_index[nsl]
+                ts = picked_times[picked_index == index]
+                its = (num.round((ts - tr.tmin)/ tr.deltat)).astype(num.int64)
+                its = its[num.logical_and(0 <= its, its < tr.data_len())]
+                ydata = num.zeros(tr.data_len())
+                ydata[its] = 1.0
+                tr.set_ydata(ydata)
+                trs_filt.append(tr)
+
+            except KeyError:
+                pass
+
+        trs_by_nsl = common.grouped_by(trs_filt, lambda tr: tr.nslc_id[:3])
+
+        dataset = []
+        for nsl in sorted(trs_by_nsl.keys()):
+            sumtr = None
+            sumn = 0
+            for tr in trs_by_nsl[nsl]:
+                if sumtr is None:
+                    sumtr = tr.copy()
+                    sumtr.set_codes(channel='CF')
+                    sumn = 1
+                else:
+                    sumtr.add(tr)
+                    sumn += 1
+
+            sumtr.ydata /= sumn
+
+            ntap = int(num.round(1./fsmooth / tr.deltat))
+            taper = num.hanning(ntap)
+            sumtr.shift(-(ntap/2.*tr.deltat))
+
+            sumtr.set_ydata(
+                num.convolve(sumtr.get_ydata()/len(trs_by_nsl), taper))
+
+            sumtr.downsample_to(deltat_cf, snap=True, demean=False)
+            sumtr.chop(wmin - tpad_new, wmax + tpad_new)
+            if num.any(sumtr != 0.):
+                dataset.append((nsl, sumtr))
+
+        return dataset
+
+
 __all__ = [
     'IFC',
     'WavePacketIFC',

src/plot.py

@@ -226,25 +226,25 @@ def plot_detection(
                 dist + scalefactor*amp + ishifter*scalefactor*0.1,
                 color=color)
 
-    for tr in trs_raw:
-        istation = station_index[tr.nslc_id[:3]]
-        dist = distances[imax, istation]
-        shift = shift_table[imax, istation]
+        for tr in trs_raw:
+            istation = station_index[tr.nslc_id[:3]]
+            dist = distances[imax, istation]
+            shift = shift_table[imax, istation]
 
-        tr = tr.copy()
+            tr = tr.copy()
 
-        tr.highpass(4, fmin, demean=True)
-        tr.lowpass(4, fmax, demean=False)
+            tr.highpass(4, fmin, demean=True)
+            tr.lowpass(4, fmax, demean=False)
 
-        tr.chop(
-            tpeak_current - 0.5*tduration + shift_min,
-            tpeak_current + 0.5*tduration + shift_max)
+            tr.chop(
+                tpeak_current - 0.5*tduration + shift_min,
+                tpeak_current + 0.5*tduration + shift_max)
 
-        t = tr.get_xdata()
-        amp = tr.get_ydata().astype(num.float)
-        amp = amp / num.max(num.abs(amp))
-        axes4.plot(t-t0, dist + scalefactor*amp, color='black')
-        axes4.plot(shift, dist, '|', mew=2, mec=color, ms=10)
+            t = tr.get_xdata()
+            amp = tr.get_ydata().astype(num.float)
+            amp = amp / num.max(num.abs(amp))
+            axes4.plot(t-t0, dist + scalefactor*amp, color='black', alpha=0.5)
+            axes4.plot(shift, dist, '|', mew=2, mec=color, ms=10)
 
     nframes = frames.shape[1]