Merged in rms (pull request #18)

Rms

backtrackbb/bp_types.py

@@ -31,7 +31,8 @@ def __init__(self,
                  max_grid=None,
                  ntraces=None,
                  beg_win=None, end_win=None,
-                 center_win=None):
+                 center_win=None,
+                 RMS_P=None, RMS_S=None):
         self.eventid = None
         self.picks = []
         self.x = x
@@ -50,6 +51,8 @@ def __init__(self,
         self.lon = None
         self.origin_time = None
         self.valid = True
+        self.rms_p = RMS_P
+        self.rms_s = RMS_S
 
     def __str__(self):
         s = '%s ' % self.eventid
@@ -67,6 +70,10 @@ def __str__(self):
             s += ' LAT %.5f LON %.5f' % (self.lat, self.lon)
         if (self.origin_time is not None):
             s += ' T_ORIG %s' % (self.origin_time)
+        if (self.rms_p is not None):
+            s += ' RMS-P %.2f' % self.rms_p
+        if (self.rms_s is not None):
+            s += ' RMS-S %.2f' % self.rms_s
         return s
 
     def from_str(self, string):
@@ -88,6 +95,14 @@ def from_str(self, string):
         self.lat = float(word[16])
         self.lon = float(word[18])
         self.origin_time = UTCDateTime(word[20])
+        try:
+            self.rms_p = float(word[22])
+        except IndexError:
+            self.rms_p = None
+        try:
+            self.rms_s = float(word[24])
+        except IndexError:
+            self.rms_s = None
 
     def add_pick(self, pick):
         self.picks.append(pick)
@@ -145,6 +160,18 @@ def compute_origin_time(self, dt_min, update_evid=False):
             pick.theor_time = pick.theor_time_absolute - otime
         self.origin_time = otime
 
+    def compute_rms(self, wave_type):
+        for wave in wave_type:
+            sq_mean = []
+            for pick in self.picks:
+                if pick.arrival_type == wave and pick.pick_time != -10.0:
+                    sq_mean.append((pick.pick_time - pick.theor_time)**2)
+            rms = np.sqrt(np.mean(sq_mean))
+            if wave == 'P':
+                self.rms_p = rms
+            elif wave == 'S':
+                self.rms_s = rms
+
     def set_eventid(self, eventid=None):
         if eventid is None:
             try:

backtrackbb/configspec.conf

@@ -54,6 +54,7 @@ start_time = string(default=None)
 end_time = string(default=None)
 #------------------------------------------------------------------------
 ignore_noisy_CF = boolean(default=False)
+n_ignore_noisy_CF = integer(default=2)
 grid_power = string(default='1')
 grid_power_ellipsoid = string(default=None)
 trigger = float(default=None)

backtrackbb/mod_btbb.py

@@ -76,10 +76,10 @@ def _run_btbb(config, st, st_CF, t_begin,
 
         min_sum = min(s[2] for s in sums)
 
-        # A noisy CF has an integral at least two times larger than
-        # the smaller one
-        # TODO: parametrize?
-        noisy_sta_wave = [(s[0], s[1]) for s in sums if s[2] >= 2*min_sum]
+        # A noisy CF has an integral at least N (=config.n_ignore_noisy_CF)
+        # times larger than the smaller one
+        noisy_sta_wave = [(s[0], s[1]) for s in sums if s[2] >=
+                          config.n_ignore_noisy_CF*min_sum]
     else:
         noisy_sta_wave = []
 
@@ -272,6 +272,7 @@ def _run_btbb(config, st, st_CF, t_begin,
             phase2 = list(set(config.grid_type) - set(config.wave_type))
             trigger.make_picks(config.stations, phase2, grids=GRD_sta)
         trigger.compute_origin_time(config.dt_min)
+        trigger.compute_rms(config.wave_type)
         trigger.set_eventid()
         trigger.check_validity()
 

backtrackbb/mod_filter_picker.py

@@ -188,8 +188,8 @@ def MBfilter_CF(st, frequencies,
             # Define the decay constant
             rosenberger_decay_constant = 1 / rosenberger_decay_nsmps
 
-            print('Rosenberger in process {}/{}\r'.format(n+1, Nb),
-                  sys.stdout.flush())
+            # print('Rosenberger in process {}/{}\r'.format(n+1, Nb),
+            #       sys.stdout.flush())
 
             # third value returned by rosenberger() is the polarizaion filter,
             # which we do not use here

backtrackbb/mod_group_trigs.py

@@ -12,18 +12,29 @@ def trig_dist(trg1, trg2):
 
 
 def group_triggers(config, triggers):
+    # Sorting triggers in the 
+    sorted_trg = sorted(triggers, key=lambda x: x.max_grid,
+                        reverse=True)
 
-    for n, trg1 in enumerate(triggers):
-        # we assume (for speed) than no more than
-        # 10 triggers need to be grouped
-        for trg2 in triggers[n+1:n+10]:
+    for n, trg1 in enumerate(sorted_trg):
+        # Currently, we compare all the triggers, 
+        # because trigger list was sorted.
+        #
+        # For speed, number of triggers to be grouped can be reduced 
+        # to a number N (ex N=50) 
+        # for trg2 in sorted_trg[n+1:n+50]:
+
+        for trg2 in sorted_trg[n+1:]:
+
             dist, time_diff = trig_dist(trg1, trg2)
+
             if (dist <= config.group_min_dist and
                     time_diff <= config.group_min_time_diff):
+
                 if trg2.max_grid <= trg1.max_grid:
-                    triggers.remove(trg2)
+                    sorted_trg.remove(trg2)
                 else:
-                    triggers.remove(trg1)
+                    sorted_trg.remove(trg1)
                     break
 
-    return triggers
+    return sorted(sorted_trg, key=lambda x: x.eventid)

backtrackbb/scripts/bt2eventdata.py

@@ -74,7 +74,7 @@ def main():
             for pick in trigger.picks:
                 f.write('%-6s ?    ?    ? %-6s ? ' %
                         (pick.station, pick.arrival_type))
-                if pick.arrival_type is 'P':
+                if pick.arrival_type == 'P':
                     time = trigger.origin_time + pick.pick_time
                 else:
                     time = trigger.origin_time + pick.theor_time
@@ -92,7 +92,7 @@ def main():
 
         # -----cutting events and saving data in specified format-----
         for pick in trigger.picks:
-            if pick.arrival_type is not 'P':
+            if pick.arrival_type != 'P':
                 continue
 
             trim_starttime = (trigger.origin_time + pick.theor_time