[fix] restitute waveform data prior to Wood-Anderson simulation

QtPyLoT.py

@@ -126,6 +126,11 @@ class MainWindow(QMainWindow):
         # load and display waveform data
         self.dirty = False
         self.load_data()
+        finv = settings.value("inventoryFile", None)
+        if finv is not None:
+            self._metadata = read_metadata(finv)
+        else:
+            self._metadata = None
         if self.loadWaveformData():
             self.updateFilterOptions()
         else:
@@ -370,6 +375,17 @@ class MainWindow(QMainWindow):
 
         self.setCentralWidget(_widget)
 
+
+    @property
+    def metadata(self):
+        return self._metadata
+
+
+    @metadata.setter
+    def metadata(self, value):
+        self._metadata = value
+
+
     def updateFileMenu(self):
 
         self.fileMenu.clear()
@@ -996,14 +1012,14 @@ class MainWindow(QMainWindow):
                 wf = select_for_phase(self.get_data().getWFData().select(
                     station=station), a.phase)
                 delta = degrees2kilometers(a.distance)
-                wapp = calc_woodanderson_pp(wf, pick.time, self.inputs.get(
+                wapp = calc_woodanderson_pp(wf, self.metadata, pick.time,
-                    'sstop'))
+                                            self.inputs.get('sstop'))
                 # using standard Gutenberg-Richter relation
                 # TODO make the ML calculation more flexible by allowing
                 # use of custom relation functions
                 mag = np.log10(wapp) + gutenberg_richter_relation(delta)
                 mags[station] = mag
-            mag = np.median([M[1] for M in mags.values()])
+            mag = np.median([M for M in mags.values()])
             # give some information on the processing
             print('number of stations used: {0}\n'.format(len(mags.values())))
             print('stations used:\n')
@@ -1020,7 +1036,7 @@ class MainWindow(QMainWindow):
             o = e.origins[0]
             mags = dict()
             fninv = settings.value("inventoryFile", None)
-            if fninv is None:
+            if fninv is None and not self.metadata:
                 fninv, _ = QFileDialog.getOpenFileName(self, self.tr(
                     "Select inventory..."), self.tr("Select file"))
                 ans = QMessageBox.question(self, self.tr("Make default..."),
@@ -1032,7 +1048,7 @@ class MainWindow(QMainWindow):
                 if ans == QMessageBox.Yes:
                     settings.setValue("inventoryFile", fninv)
                     settings.sync()
-            metadata = read_metadata(fninv)
+                self.metadata = read_metadata(fninv)
             for a in o.arrivals:
                 if a.phase in 'sS':
                     continue

pylot/core/analysis/magnitude.py

@@ -271,7 +271,7 @@ class M0Mw(Magnitude):
                 self.picdic = picks
 
 
-def calc_woodanderson_pp(st, T0, win=10., verbosity=False):
+def calc_woodanderson_pp(st, metadata, T0, win=10., verbosity=False):
     if verbosity:
         print ("Getting Wood-Anderson peak-to-peak amplitude ...")
         print ("Simulating Wood-Anderson seismograph ...")
@@ -288,17 +288,27 @@ def calc_woodanderson_pp(st, T0, win=10., verbosity=False):
     stime, etime = common_range(st)
     st.trim(stime, etime)
 
+    invtype, inventory = metadata
+
+    st = restitute_data(st, invtype, inventory)
+
     dt = st[0].stats.delta
 
     st.simulate(paz_remove=None, paz_simulate=paz_wa)
 
     # get RMS of both horizontal components
-    sqH = np.sqrt(np.sum([np.power(tr.data, 2) for tr in st if
+    power = [np.power(tr.data, 2) for tr in st if tr.stats.channel[-1] not
-                          tr.stats.channel[-1] not in 'Z3']))
+              in 'Z3']
+    if len(power) != 2:
+        raise ValueError('Wood-Anderson amplitude defintion only valid for '
+                         'two horizontals: {0} given'.format(len(power)))
+    power_sum = power[0] + power[1]
+    sqH = np.sqrt(power_sum)
+
     # get time array
     th = np.arange(0, len(sqH) * dt, dt)
     # get maximum peak within pick window
-    iwin = getsignalwin(th, T0, win)
+    iwin = getsignalwin(th, T0 - stime, win)
     wapp = np.max(sqH[iwin])
     if verbosity:
         print("Determined Wood-Anderson peak-to-peak amplitude: {0} "

pylot/core/pick/autopick.py

@@ -604,7 +604,7 @@ def autopickstation(wfstream, pickparam, verbose=False):
             hdat += ndat
             h_copy = hdat.copy()
             [cordat, restflag] = restitute_data(h_copy, invtype, inventory)
-            if restflag == 1:
+            if restflag is True:
                 # calculate WA-peak-to-peak amplitude
                 # using subclass WApp of superclass Magnitude
                 wapp = WApp(cordat, mpickP, mpickP + sstop + (0.5 * (mpickP