Extended class MoMw for calculating source spectrum. New functions calcsourcespec, calcMoMw and run_calcMoMw implemented.

pylot/core/analysis/magnitude.py

@@ -8,10 +8,12 @@ Created August/September 2015.
 
 import matplotlib.pyplot as plt
 import numpy as np
-from obspy.core import Stream
+from obspy.core import Stream, UTCDateTime
-from pylot.core.pick.utils import getsignalwin
+from pylot.core.pick.utils import getsignalwin, crossings_nonzero_all
 from pylot.core.util.utils import getPatternLine
 from scipy.optimize import curve_fit
+from scipy import integrate
+from pylot.core.read.data import Data
 
 class Magnitude(object):
     '''
@@ -20,7 +22,8 @@ class Magnitude(object):
     and moment magnitudes.
     '''
 
-    def __init__(self, wfstream, To, pwin, iplot, NLLocfile=None, picks=None, rho=None, vp=None):
+    def __init__(self, wfstream, To, pwin, iplot, NLLocfile=None, \
+                 picks=None, rho=None, vp=None, invdir=None):
         '''
         :param: wfstream
         :type: `~obspy.core.stream.Stream
@@ -30,7 +33,7 @@ class Magnitude(object):
 
         :param: pwin, pick window [To To+pwin] to get maximum
                 peak-to-peak amplitude (WApp) or to calculate
-                source spectrum (DCfc)
+                source spectrum (DCfc) around P onset
         :type:  float
 
         :param: iplot, no. of figure window for plotting interims results
@@ -48,6 +51,9 @@ class Magnitude(object):
 
         :param: vp [m/s], P-velocity
         :param: integer
+
+        :param: invdir, path to inventory or dataless-SEED file
+        :type:  string
         '''
 
         assert isinstance(wfstream, Stream), "%s is not a stream object" % str(wfstream)
@@ -60,6 +66,7 @@ class Magnitude(object):
         self.setrho(rho)
         self.setpicks(picks)
         self.setvp(vp)
+        self.setinvdir(invdir)
         self.calcwapp()
         self.calcsourcespec()
         self.run_calcMoMw()
@@ -122,6 +129,12 @@ class Magnitude(object):
     def getfc(self):
         return self.fc
 
+    def setinvdir(self, invdir):
+        self.invdir = invdir
+
+    def getinvdir(self):
+        return self.invdir
+
     def getpicdic(self):
         return self.picdic
 
@@ -190,7 +203,8 @@ class M0Mw(Magnitude):
     Requires results of class w0fc for calculating plateau w0 
     and corner frequency fc of source spectrum, respectively. Uses 
     subfunction calcMoMw.py. Returns modified dictionary of picks including 
-    seismic moment Mo and corresponding moment magntiude Mw.
+    Dc-value, corner frequency fc, seismic moment Mo and 
+    corresponding moment magntiude Mw.
     '''
 
     def run_calcMoMw(self):
@@ -198,12 +212,15 @@ class M0Mw(Magnitude):
         picks = self.getpicks()
         nllocfile = self.getNLLocfile()
         wfdat = self.getwfstream()
+        # get vertical component data only
+        zdat = wfdat.select(component="Z")
+        if len(zdat) == 0:  # check for other components
+            zdat = wfdat.select(component="3")
+
         for key in picks:
-             if picks[key]['P']['weight'] < 4 and picks[key]['P']['w0'] is not None:
+             if picks[key]['P']['weight'] < 4:
                  # select waveform
-                 selwf = wfdat.select(station=key)
+                 selwf = zdat.select(station=key)
-                 # get corresponding height of source spectrum plateau w0
-                 w0 = picks[key]['P']['w0']
                  # get hypocentral distance of station
                  # from NLLoc-location file
                  if len(key) > 4:
@@ -214,14 +231,27 @@ class M0Mw(Magnitude):
                      Ppattern = '%s    ?    ?    ? P' % key
                  nllocline = getPatternLine(nllocfile, Ppattern)
                  delta = float(nllocline.split(None)[21])
-                 # call subfunction
+                 # call subfunction to estimate source spectrum
-                 [Mo, Mw] = calcMoMw(selwf, w0, self.getrho(), self.getvp(), delta)
+                 # and to derive w0 and fc 
-                 # add Mo and Mw to dictionary
+                 [w0, fc] = calcsourcespec(selwf, picks[key]['P']['mpp'], \
+                             self.getiplot(), self.getinvdir())
+
+                 if w0 is not None:
+                     # call subfunction to calculate Mo and Mw
+                     [Mo, Mw] = calcMoMw(selwf, w0, self.getrho(), self.getvp(), \
+                                 delta, self.getinvdir())
+                 else:
+                     Mo = None
+                     Mw = None
+
+                 # add w0, fc, Mo and Mw to dictionary
+                 picks[key]['P']['w0'] = w0
+                 picks[key]['P']['fc'] = fc
                  picks[key]['P']['Mo'] = Mo
                  picks[key]['P']['Mw'] = Mw
                  self.picdic = picks
 
-def calcMoMw(wfstream, w0, rho, vp, delta):
+def calcMoMw(wfstream, w0, rho, vp, delta, inv):
     '''
     Subfunction of run_calcMoMw to calculate individual
     seismic moments and corresponding moment magnitudes.
@@ -245,25 +275,59 @@ def calcMoMw(wfstream, w0, rho, vp, delta):
     return Mo, Mw
 
         
-class w0fc(Magnitude):
-    '''
-    Method to calculate the source spectrum and to derive from that the plateau
-    (usually called omega0) and the corner frequency assuming Aki's omega-square
-    source model. Has to be derived from instrument corrected displacement traces!
-    '''
 
-    def calcsourcespec(self):
+def calcsourcespec(wfstream, onset, iplot, inventory):
+    '''
+    Subfunction to calculate the source spectrum and to derive from that the plateau
+    (usually called omega0) and the corner frequency assuming Aki's omega-square
+    source model. Has to be derived from instrument corrected displacement traces,
+    thus restitution and integration necessary!
+    '''
     print ("Calculating source spectrum ....")
 
-        self.w0 = None # DC-value
+    fc = None
-        self.fc = None # corner frequency
+    w0 = None
+    data = Data()
+    z_copy = wfstream.copy()
 
-        stream = self.getwfstream()
+    [corzdat, restflag] = data.restituteWFData(inventory, z_copy)
-        tr = stream[0]
 
+    if restflag == 1:
+        # integrate to displacment
+        corintzdat = integrate.cumtrapz(corzdat[0], None, corzdat[0].stats.delta)
+        z_copy[0].data = corintzdat
+        tr = z_copy[0]
+        # get window after P pulse for 
+        # calculating source spectrum
+        if tr.stats.sampling_rate <= 100:
+            winzc = tr.stats.sampling_rate
+        elif tr.stats.sampling_rate > 100 and \
+              tr.stats.sampling_rate <= 200:
+             winzc = 0.5 * tr.stats.sampling_rate
+        elif tr.stats.sampling_rate > 200 and \
+              tr.stats.sampling_rate <= 400:
+             winzc = 0.2 * tr.stats.sampling_rate
+        elif tr.stats.sampling_rate > 400:
+             winzc = tr.stats.sampling_rate
+        tstart = UTCDateTime(tr.stats.starttime)
+        tonset = onset.timestamp -tstart.timestamp
+        impickP = tonset * tr.stats.sampling_rate
+        wfzc = tr.data[impickP : impickP + winzc]
         # get time array
         t = np.arange(0, len(tr) * tr.stats.delta, tr.stats.delta)
-        iwin = getsignalwin(t, self.getTo(), self.getpwin())
+        # calculate spectrum using only first cycles of
+        # waveform after P onset!
+        zc = crossings_nonzero_all(wfzc)
+        if np.size(zc) == 0 or len(zc) <= 3:
+            print ("Something is wrong with the waveform, "
+                   "no zero crossings derived!")
+            print ("No calculation of source spectrum possible!")
+            plotflag = 0
+        else:
+            plotflag = 1
+            index = min([3, len(zc) - 1])
+            calcwin = (zc[index] - zc[0]) * z_copy[0].stats.delta
+            iwin = getsignalwin(t, tonset, calcwin)
             xdat = tr.data[iwin]
 
             # fft
@@ -301,30 +365,36 @@ class w0fc(Magnitude):
                    "Determined corner frequency: %f Hz" % (w01, fc1))
         
             # use of conventional fitting 
-        [w02, fc2] = fitSourceModel(F, YY.real, Fcin, self.getiplot())
+            [w02, fc2] = fitSourceModel(F, YY.real, Fcin, iplot)
  
             # get w0 and fc as median 
-        self.w0 = np.median([w01, w02])
+            w0 = np.median([w01, w02])
-        self.fc = np.median([fc1, fc2])
+            fc = np.median([fc1, fc2])
-        print("w0fc: Using w0-value = %e m/Hz and fc = %f Hz" % (self.w0, self.fc))
+            print("w0fc: Using w0-value = %e m/Hz and fc = %f Hz" % (w0, fc))
 
-	if self.getiplot() > 1:
+    if iplot > 1:
         f1 = plt.figure()
         plt.subplot(2,1,1)
         # show displacement in mm
         plt.plot(t, np.multiply(tr, 1000), 'k')
+        if plotflag == 1:
             plt.plot(t[iwin], np.multiply(xdat, 1000), 'g')
-            plt.title('Seismogram and P pulse, station %s' % tr.stats.station)
+            plt.title('Seismogram and P Pulse, Station %s-%s' \
+                        % (tr.stats.station, tr.stats.channel))
+        else:
+            plt.title('Seismogram, Station %s-%s' \
+                        % (tr.stats.station, tr.stats.channel))
         plt.xlabel('Time since %s' % tr.stats.starttime)
         plt.ylabel('Displacement [mm]')
 
+        if plotflag == 1:
             plt.subplot(2,1,2)
             plt.loglog(f, Y.real, 'k')
             plt.loglog(F, YY.real)
             plt.loglog(F, fit, 'g')
-            plt.loglog([self.fc, self.fc], [self.w0/100, self.w0], 'g') 
+            plt.loglog([fc, fc], [w0/100, w0], 'g') 
             plt.title('Source Spectrum from P Pulse, w0=%e m/Hz, fc=%6.2f Hz' \
-                       % (self.w0, self.fc))
+                       % (w0, fc))
             plt.xlabel('Frequency [Hz]')
             plt.ylabel('Amplitude [m/Hz]')
             plt.grid()
@@ -332,6 +402,7 @@ class w0fc(Magnitude):
         raw_input()
         plt.close(f1)
 
+    return w0, fc
      
 
 def synthsourcespec(f, omega0, fcorner):