[new] added richter magnitude calculation (to be tested)

QtPyLoT.py

@@ -42,12 +42,13 @@ from obspy import UTCDateTime
 from obspy.geodetics import degrees2kilometers
 from obspy.core.event import Magnitude
 
-from pylot.core.analysis.magnitude import calcsourcespec, calcMoMw
+from pylot.core.analysis.magnitude import calcsourcespec, calcMoMw, \
+    calc_woodanderson_pp, gutenberg_richter_relation
 from pylot.core.io.data import Data
 from pylot.core.io.inputs import FilterOptions, AutoPickParameter
 from pylot.core.pick.autopick import autopickevent
 from pylot.core.pick.compare import Comparison
-from pylot.core.pick.utils import symmetrize_error
+from pylot.core.pick.utils import symmetrize_error, select_for_phase
 from pylot.core.io.phases import picksdict_from_picks
 import pylot.core.loc.nll as nll
 from pylot.core.util.defaults import FILTERDEFAULTS, COMPNAME_MAP, \
@@ -972,7 +973,47 @@ class MainWindow(QMainWindow):
         self.get_data().applyEVTData(lt.read_location(locpath), type='event')
         self.get_data().get_evt_data().magnitudes.append(self.calc_magnitude())
 
-    def calc_magnitude(self):
+
+    def calc_magnitude(self, type='ML'):
+        if type == 'ML':
+            return self.calc_richter_magnitude()
+        elif type == 'Mw':
+            return self.calc_moment_magnitude()
+        else:
+            return None
+
+
+    def calc_richter_magnitude(self):
+        e = self.get_data().get_evt_data()
+        if e.origins:
+            o = e.origins[0]
+            mags = dict()
+            for a in o.arrivals:
+                if a.phase not in 'sS':
+                    continue
+                pick = a.pick_id.get_referred_object()
+                station = pick.waveform_id.station_code
+                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(
+                    '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()])
+            # give some information on the processing
+            print('number of stations used: {0}\n'.format(len(mags.values())))
+            print('stations used:\n')
+            for s in mags.keys(): print('\t{0}'.format(s))
+
+            return Magnitude(mag=mag, magnitude_type='ML')
+        else:
+            return None
+
+    def calc_moment_magnitude(self):
         e = self.get_data().get_evt_data()
         settings = QSettings()
         if e.origins:

pylot/core/analysis/magnitude.py

@@ -271,6 +271,41 @@ class M0Mw(Magnitude):
                 self.picdic = picks
 
 
+def calc_woodanderson_pp(st, T0, win=10., verbosity=False):
+    if verbosity:
+        print ("Getting Wood-Anderson peak-to-peak amplitude ...")
+        print ("Simulating Wood-Anderson seismograph ...")
+
+    # poles, zeros and sensitivity of WA seismograph
+    # (see Uhrhammer & Collins, 1990, BSSA, pp. 702-716)
+    paz_wa = {
+        'poles': [5.6089 - 5.4978j, -5.6089 - 5.4978j],
+        'zeros': [0j, 0j],
+        'gain': 2080,
+        'sensitivity': 1
+    }
+
+    stime, etime = common_range(st)
+    st.trim(stime, etime)
+
+    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
+                          tr.stats.channel[-1] not in 'Z3']))
+    # get time array
+    th = np.arange(0, len(sqH) * dt, dt)
+    # get maximum peak within pick window
+    iwin = getsignalwin(th, T0, win)
+    wapp = np.max(sqH[iwin])
+    if verbosity:
+        print("Determined Wood-Anderson peak-to-peak amplitude: {0} "
+              "mm".format(wapp))
+    return wapp
+
+
 def calcMoMw(wfstream, w0, rho, vp, delta):
     '''
     Subfunction of run_calcMoMw to calculate individual
@@ -358,7 +393,6 @@ def calcsourcespec(wfstream, onset, metadata, vp, delta, azimuth, incidence,
 
     fc = None
     w0 = None
-    data = Data()
     wf_copy = wfstream.copy()
 
     invtype, inventory = metadata