implemented basic template migration method

src/config.py

@@ -104,7 +104,7 @@ class Config(Object):
                 spacing = vmin / fsmooth_max / self.autogrid_density_factor
 
                 lat0, lon0, north, east, depth = geo.bounding_box_square(
-                    *receiver.receivers_coords(receivers),
+                    *geo.points_coords(receivers),
                     scale=self.autogrid_radius_factor)
 
                 self._grid = grid.Carthesian3DGrid(

src/core.py

@@ -37,20 +37,32 @@ 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(
+                'waveform data path does not exist: %s' % data_path)
+
+    p = pile.make_pile(config.data_paths, fileformat='detect')
+    if p.is_empty():
+        raise common.LassieError('no usable waveforms found')
+
     ifcs = config.image_function_contributions
+    for ifc in ifcs:
+        ifc.prescan(p)
 
     shift_tables = []
     tshift_maxs = []
     for ifc in ifcs:
-        shift_tables.append(ifc.shifter.get_table(grid, receivers))
+        shift_tables.append(ifc.get_table(grid, receivers))
         tshift_maxs.append(shift_tables[-1].max())
 
     tshift_max = max(tshift_maxs) * 1.0
 
-    tinc = tshift_max * 10
-
     tpad = max(ifc.get_tpad() for ifc in ifcs) + tshift_max
 
+    tinc = tshift_max * 10 + 3*tpad
+
     fsmooth_min = min(ifc.get_fsmooth() for ifc in ifcs)
 
     blacklist = set(tuple(s.split('.')) for s in config.blacklist)
@@ -59,15 +71,6 @@ def scan(
         (rec.codes, i) for (i, rec) in enumerate(receivers)
         if rec.codes not in blacklist)
 
-    for data_path in config.data_paths:
-        if not op.exists(data_path):
-            raise common.LassieError(
-                'waveform data path does not exist: %s' % data_path)
-
-    p = pile.make_pile(config.data_paths, fileformat='detect')
-    if p.is_empty():
-        raise common.LassieError('no usable waveforms found')
-
     ngridpoints = grid.size()
 
     idetection = 0

src/geo.py

 import numpy as num
 from pyrocko import orthodrome as od
+from pyrocko.guts import Object, Float
+
+
+class Point(Object):
+    lat = Float.T(default=0.0)
+    lon = Float.T(default=0.0)
+    x = Float.T(default=0.0)
+    y = Float.T(default=0.0)
+    z = Float.T(default=0.0)
+
+    @property
+    def vec(self):
+        return self.lat, self.lon, self.x, self.y, self.z
+
+
+def points_coords(points, system=None):
+    lats, lons, xs, ys, zs = num.array(
+        [r.vec for r in points], dtype=num.float).T
+
+    if system is None:
+        return (lats, lons, xs, ys, zs)
+
+    elif system == 'latlon':
+        return od.ne_to_latlon(lats, lons, xs, ys)
+
+    elif system[0] == 'ne':
+        lat0, lon0 = system[1:3]
+        if num.all(lats == lat0) and num.all(lons == lon0):
+            return xs, ys
+        else:
+            elats, elons = od.ne_to_latlon(lats, lons, xs, ys)
+            return od.latlon_to_ne_numpy(lat0, lon0, elats, elons)
 
 
 def float_array_broadcast(*args):

src/grid.py

@@ -3,7 +3,7 @@ import numpy as num
 from pyrocko.guts import Object, Float
 from pyrocko import orthodrome as od
 
-from lassie import receiver
+from lassie import geo
 
 guts_prefix = 'lassie'
 
@@ -68,7 +68,7 @@ class Carthesian3DGrid(Grid):
         nx, ny, nz = self._shape()
         x, y, z = self._get_coords()
 
-        rx, ry = receiver.receivers_coords(
+        rx, ry = geo.points_coords(
             receivers, system=('ne', self.lat, self.lon))
 
         rz = num.array([r.z for r in receivers], dtype=num.float)

src/ifc.py

 import sys
 import logging
+from collections import defaultdict
 import numpy as num
-from pyrocko.guts import Object, String, Float
-from pyrocko import trace, autopick, util
-from lassie import shifter, common
+from pyrocko.guts import Object, String, Float, Bool, StringChoice
+from pyrocko import trace, autopick, util, model, gui_util
+from lassie import shifter, common, geo
 
 logger = logging.getLogger('lassie.ifc')
 
@@ -19,9 +20,18 @@ class IFC(Object):
     fmax = Float.T()
     shifter = shifter.Shifter.T()
 
+    def get_table(self, grid, receivers):
+        return self.shifter.get_table(grid, receivers)
+
     def get_tpad(self):
         return 4. / self.fmin
 
+    def get_fsmooth(self):
+        return self.fmin
+
+    def prescan(self, p):
+        pass
+
     def preprocess(self, trs, wmin, wmax, tpad_new):
         pass
 
@@ -184,8 +194,149 @@ class OnsetIFC(IFC):
         return dataset
 
 
+class TemplateMatchingIFC(IFC):
+
+    template_event_path = String.T(
+        help='Event parameters of the template')
+
+    template_markers_path = String.T(
+        optional=True,
+        help='File with markers defining the template')
+
+    sum_square = Bool.T(
+        default=False,
+        help='Sum square of correlation')
+
+    normalization = StringChoice.T(
+        default='gliding',
+        choices=['off', 'normal', 'gliding'])
+
+    downsample_rate = Float.T(
+        optional=True,
+        help='If set, downsample to this sampling rate before processing [Hz]')
+
+    def get_tpad(self):
+        tmin_masters = min(tr.tmin for tr in self.masters.values())
+        tmax_masters = max(tr.tmax for tr in self.masters.values())
+        tmaster = tmax_masters - tmin_masters
+        return tmaster
+
+    def get_template_origin(self):
+        event = model.load_one_event(self.template_event_path)
+        origin = geo.Point(
+            lat=event.lat,
+            lon=event.lon,
+            z=event.depth)
+
+        return origin
+
+    def get_table(self, grid, receivers):
+        origin = self.get_template_origin()
+        return self.shifter.get_offset_table(grid, receivers, origin)
+
+    def extract_template(self, p):
+
+        markers = gui_util.load_markers(self.template_markers_path)
+
+        trace_selector_global = lambda tr: True
+
+        period_highpass = 1./self.fmin
+        tpad = 2 * period_highpass
+
+        master_traces = []
+        for marker in markers:
+            if not marker.nslc_ids:
+                trace_selector = trace_selector_global
+            else:
+                trace_selector = lambda tr: (
+                    marker.match_nslc(tr.nslc_id) and
+                    trace_selector_global(tr))
+
+            master_traces.extend(p.all(
+                tmin=marker.tmin,
+                tmax=marker.tmax,
+                trace_selector=trace_selector,
+                tpad=tpad))
+
+        masters = {}
+        for xtr in master_traces:
+            tr = xtr.copy()
+            if self.downsample_rate is not None:
+                tr.downsample_to(
+                    1./self.downsample_rate, demean=False, snap=True)
+
+            tr.highpass(4, self.fmin, demean=False)
+            tr.lowpass(4, self.fmax, demean=False)
+            smin = round(xtr.wmin / tr.deltat) * tr.deltat
+            smax = round(xtr.wmax / tr.deltat) * tr.deltat
+            tr.chop(smin, smax)
+            if tr.nslc_id in masters:
+                raise common.LassieError(
+                    'more than one waveform selected on trace with id "%s"'
+                    % '.'.join(tr.nslc_id))
+
+            masters[tr.nslc_id] = tr
+
+        return masters
+
+    def prescan(self, p):
+        self.masters = self.extract_template(p)
+
+    def preprocess(self, trs, wmin, wmax, tpad_new):
+
+        tmin_masters = min(tr.tmin for tr in self.masters.values())
+        tmax_masters = max(tr.tmax for tr in self.masters.values())
+        tmaster = tmax_masters - tmin_masters
+        tref = tmin_masters
+
+        nsl_to_traces = defaultdict(list)
+        for orig_b in trs:
+
+            b = orig_b.copy()
+
+            nslc = b.nslc_id
+            a = self.masters[nslc]
+            if self.downsample_rate is not None:
+                b.downsample_to(
+                    1./self.downsample_rate, demean=False, snap=True)
+
+            b.highpass(4, self.fmin, demean=False)
+            b.lowpass(4, self.fmax, demean=False)
+            smin = round((wmin - tmaster) / b.deltat) * b.deltat
+            smax = round((wmax + tmaster) / b.deltat) * b.deltat
+            b.chop(smin, smax)
+
+            normalization = self.normalization
+            if normalization == 'off':
+                normalization = None
+
+            c = trace.correlate(
+                a, b, mode='valid', normalization=normalization)
+
+            c.shift(-c.tmin + b.tmin - (a.tmin - tref))
+            c.meta = {'tabu': True}
+            #c.set_codes(location='cc')
+            if self.sum_square:
+                c.ydata = c.ydata**2
+
+            c.chop(wmin - tpad_new, wmax + tpad_new)
+
+            nsl_to_traces[nslc[:3]].append(c)
+
+        dataset = []
+        for nsl, cs in nsl_to_traces.iteritems():
+            csum = cs[0]
+            for c in cs[1:]:
+                csum.add(c)
+
+            dataset.append((nsl, csum))
+
+        return dataset
+
+
 __all__ = [
     'IFC',
     'WavePacketIFC',
     'OnsetIFC',
+    'TemplateMatchingIFC',
 ]

src/plot.py

 import time
 import numpy as num
 from pyrocko import automap
-from lassie import receiver, grid as gridmod, geo
+from lassie import grid as gridmod, geo
 
 
 def map_gmt(
@@ -50,7 +50,7 @@ def map_geometry(config, output_path):
     grid = config.get_grid()
 
     lat0, lon0, north, east, depth = geo.bounding_box_square(
-        *receiver.receivers_coords(receivers),
+        *geo.points_coords(receivers),
         scale=config.autogrid_radius_factor)
 
     radius = max((north[1] - north[0]), (east[1] - east[0]))
@@ -62,7 +62,7 @@ def map_geometry(config, output_path):
 
 
 def plot_receivers(axes, receivers, system='latlon'):
-    x, y = receiver.receivers_coords(receivers, system=system)
+    x, y = geo.points_coords(receivers, system=system)
     axes.plot(y, x, '^', color='black', ms=10.)
 
 
@@ -253,6 +253,9 @@ def plot_detection(
 
     amax = frames[imax, iframe]
 
+    axes.set_xlim(grid.ymin, grid.ymax)
+    axes.set_ylim(grid.xmin, grid.xmax)
+
     cmap = cm.YlOrBr
     system = ('ne', grid.lat, grid.lon)
     if movie:
@@ -284,6 +287,7 @@ def plot_detection(
     else:
         frame = num.max(frames[:, iframe_min:iframe_max+1], axis=1)
         grid.plot(axes, frame, amin=0.0, amax=amax, cmap=cmap, system=system)
+
         plot_receivers(axes, receivers, system=system)
         axes.plot(
             ypeak, xpeak, '*',

src/receiver.py

-import numpy as num
-from pyrocko.guts import Object, Tuple, Float, String
-import pyrocko.orthodrome as od
+from pyrocko.guts import Tuple, String
+from lassie import geo
 
 guts_prefix = 'lassie'
 
 
-class Receiver(Object):
+class Receiver(geo.Point):
     codes = Tuple.T(3, String.T())
-    lat = Float.T(default=0.0)
-    lon = Float.T(default=0.0)
-    x = Float.T(default=0.0)
-    y = Float.T(default=0.0)
-    z = Float.T(default=0.0)
 
-    @property
-    def vec(self):
-        return self.lat, self.lon, self.x, self.y, self.z
-
-
-def receivers_coords(receivers, system=None):
-    lats, lons, xs, ys, zs = num.array(
-        [r.vec for r in receivers], dtype=num.float).T
-
-    if system is None:
-        return (lats, lons, xs, ys, zs)
-
-    elif system == 'latlon':
-        return od.ne_to_latlon(lats, lons, xs, ys)
-
-    elif system[0] == 'ne':
-        lat0, lon0 = system[1:3]
-        if num.all(lats == lat0) and num.all(lons == lon0):
-            return xs, ys
-        else:
-            elats, elons = od.ne_to_latlon(lats, lons, xs, ys)
-            return od.latlon_to_ne_numpy(lat0, lon0, elats, elons)
 
 __all__ = [
     'Receiver',
-    'receivers_coords',
 ]

src/shifter.py

+import numpy as num
 from pyrocko.guts import Object, Float
+from lassie import geo
 
 guts_prefix = 'lassie'
 
@@ -14,6 +16,26 @@ class VelocityShifter(Shifter):
         distances = grid.distances(receivers)
         return distances / self.velocity
 
+    def get_offset_table(self, grid, receivers, origin):
+        distances = grid.distances(receivers)
+
+        x, y = geo.points_coords([origin], system=('ne', grid.lat, grid.lon))
+        x = x[0]
+        y = y[0]
+        z = origin.z
+
+        rx, ry = geo.points_coords(
+            receivers, system=('ne', grid.lat, grid.lon))
+
+        rz = num.array([r.z for r in receivers], dtype=num.float)
+
+        distances_origin = num.sqrt((x - rx)**2 + (y - ry)**2 + (z - rz)**2)
+
+        na = num.newaxis
+        offset_distances = distances - distances_origin[na, :]
+
+        return offset_distances / self.velocity
+
     def get_vmin(self):
         return self.velocity