Merge branch 'develop' of ariadne.geophysik.ruhr-uni-bochum.de:/data/git/pylot into develop

Conflicts:
	pylot/core/pick/autopick.py

autoPyLoT.py

@@ -88,9 +88,10 @@ def autoPyLoT(inputfile):
             maxnumit = 3 # maximum number of iterations for re-picking
         else:
             locflag = 0
-            print ("                 !!!              ")
+            print("                 !!!              ")
-            print ("!!No location routine available, autoPyLoT is running in non-location mode!!")
+            print("!!No location routine available, autoPyLoT is running in non-location mode!!")
-            print ("                 !!!              ")
+            print("!!No source parameter estimation possible!!")
+            print("                 !!!              ")
 
 
         # multiple event processing
@@ -132,6 +133,18 @@ def autoPyLoT(inputfile):
 
                     if len(badpicks) == 0:
                         print("autoPyLoT: No bad onsets found, thus no iterative picking necessary!") 
+                        # get NLLoc-location file
+                        locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout)
+                        if len(glob.glob(locsearch)) > 0:
+                            # get latest NLLoc-location file if several are available
+                            nllocfile = max(glob.glob(locsearch), key=os.path.getctime) 
+                            # calculating seismic moment Mo and moment magnitude Mw
+                            finalpicks = M0Mw(wfdat, None, None,  parameter.getParam('iplot'), \
+                                          nllocfile, picks, parameter.getParam('rho'), \
+                                          parameter.getParam('vp'), parameter.getParam('invdir'))
+                        else:
+                            print("autoPyLoT: No NLLoc-location file available!")
+                            print("No source parameter estimation possible!")
                     else:
                         # get theoretical P-onset times from NLLoc-location file
                         locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout)
@@ -151,6 +164,9 @@ def autoPyLoT(inputfile):
                                 # locate the event
                                 locate(nlloccall, ctrfile)
                                 print("autoPyLoT: Iteration No. %d finished." % nlloccounter)
+                                # get updated NLLoc-location file
+                                nllocfile = max(glob.glob(locsearch), key=os.path.getctime) 
+                                # check for bad picks
                                 badpicks = []
                                 for key in picks:
                                      if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4:
@@ -159,23 +175,19 @@ def autoPyLoT(inputfile):
                                        len(badpicks)))
                                 if len(badpicks) == 0:
                                     print("autoPyLoT: No more bad onsets found, stop iterative picking!")
-                                    break
+                                    nlloccounter = maxnumit
-                            # calculating seismic moment Mo and corresponding moment
+
-                            # magnitude Mw after Hanks and Kanamori (1979) from reliable
+                            # calculating seismic moment Mo and moment magnitude Mw
-                            # picks/waveforms 
+                            finalpicks = M0Mw(wfdat, None, None,  parameter.getParam('iplot'), \
-                            for key in picks:
+                                          nllocfile, picks, parameter.getParam('rho'), \
-                                 if picks[key]['P']['weight'] < 4 and picks[key]['P']['w0'] is not None:
+                                          parameter.getParam('vp'), parameter.getParam('invdir'))
-                                     selwf = wfdat.select(station=key)
-                                     w0 = picks[key]['P']['w0']
-                                     sourcepara = M0Mw(selwf, None, None,  None, w0, 5, \
-                                     parameter.getParam('rho'), parameter.getParam('vp'))
                         else:
                             print("autoPyLoT: No NLLoc-location file available! Stop iteration!")
                 ##########################################################
                 # write phase files for various location routines
                 # HYPO71
                 hypo71file = '%s/%s/autoPyLoT_HYPO71.pha' % (datapath, evID)
-                writephases(picks, 'HYPO71', hypo71file)
+                writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file)
 
                 endsplash = '''------------------------------------------\n'
                                -----Finished event %s!-----\n' 
@@ -220,6 +232,18 @@ def autoPyLoT(inputfile):
 
                 if len(badpicks) == 0:
                     print("autoPyLoT: No bad onsets found, thus no iterative picking necessary!") 
+                    # get NLLoc-location file
+                    locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout)
+                    if len(glob.glob(locsearch)) > 0:
+                        # get latest NLLOc-location file if several are available
+                        nllocfile = max(glob.glob(locsearch), key=os.path.getctime) 
+                        # calculating seismic moment Mo and moment magnitude Mw
+                        finalpicks = M0Mw(wfdat, None, None,  parameter.getParam('iplot'), \
+                                      nllocfile, picks, parameter.getParam('rho'), \
+                                      parameter.getParam('vp'), parameter.getParam('invdir'))
+                    else:
+                        print("autoPyLoT: No NLLoc-location file available!")
+                        print("No source parameter estimation possible!")
                 else:
                     # get theoretical P-onset times from NLLoc-location file
                     locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout)
@@ -239,6 +263,9 @@ def autoPyLoT(inputfile):
                             # locate the event
                             locate(nlloccall, ctrfile)
                             print("autoPyLoT: Iteration No. %d finished." % nlloccounter)
+                            # get updated NLLoc-location file
+                            nllocfile = max(glob.glob(locsearch), key=os.path.getctime) 
+                            # check for bad picks
                             badpicks = []
                             for key in picks:
                                  if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4:
@@ -247,25 +274,19 @@ def autoPyLoT(inputfile):
                                    len(badpicks)))
                             if len(badpicks) == 0:
                                 print("autoPyLoT: No more bad onsets found, stop iterative picking!")
-                                break
+                                nlloccounter = maxnumit
-                        
-                        # calculating seismic moment Mo and corresponding moment
-                        # magnitude Mw after Hanks and Kanamori (1979) from reliable
-                        # picks/waveforms 
-                        for key in picks:
-                             if picks[key]['P']['weight'] < 4 and picks[key]['P']['w0'] is not None:
-                                 selwf = wfdat.select(station=key)
-                                 w0 = picks[key]['P']['w0']
-                                 sourcepara = M0Mw(selwf, None, None,  None, w0, 5, \
-                                 parameter.getParam('rho'), parameter.getParam('vp'))
                         
+                        # calculating seismic moment Mo and moment magnitude Mw
+                        finalpicks = M0Mw(wfdat, None, None,  parameter.getParam('iplot'), \
+                                      nllocfile, picks, parameter.getParam('rho'), \
+                                      parameter.getParam('vp'), parameter.getParam('invdir'))
                     else:
                         print("autoPyLoT: No NLLoc-location file available! Stop iteration!")
             ##########################################################
             # write phase files for various location routines
             # HYPO71
             hypo71file = '%s/%s/autoPyLoT_HYPO71.pha' % (datapath, parameter.getParam('eventID'))
-            writephases(picks, 'HYPO71', hypo71file)
+            writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file)
            
             endsplash = '''------------------------------------------\n'
                            -----Finished event %s!-----\n' 

pylot/core/analysis/magnitude.py

@@ -8,17 +8,22 @@ 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):
     '''
     Superclass for calculating Wood-Anderson peak-to-peak
-    amplitudes, local magnitudes and moment magnitudes.
+    amplitudes, local magnitudes, seismic moments
+    and moment magnitudes.
     '''
 
-    def __init__(self, wfstream, To, pwin, iplot, w0=None, delta=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
@@ -28,12 +33,27 @@ 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
         :type:  integer
 
+        :param: NLLocfile, name and full path to NLLoc-location file
+                needed when calling class MoMw
+        :type:  string
+
+        :param: picks, dictionary containing picking results
+        :type:  dictionary
+
+        :param: rho [kg/m³], rock density, parameter from autoPyLoT.in
+        :type:  integer
+
+        :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)
@@ -42,13 +62,14 @@ class Magnitude(object):
         self.setTo(To)
         self.setpwin(pwin)
         self.setiplot(iplot)
-        self.setw0(w0)
+        self.setNLLocfile(NLLocfile)
         self.setrho(rho)
-        self.setdelta(delta)
+        self.setpicks(picks)
         self.setvp(vp)
+        self.setinvdir(invdir)
         self.calcwapp()
         self.calcsourcespec()
-        self.calcMoMw()
+        self.run_calcMoMw()
 
 
     def getwfstream(self):
@@ -75,11 +96,11 @@ class Magnitude(object):
     def setiplot(self, iplot):
         self.iplot = iplot
 
-    def setw0(self, w0):
+    def setNLLocfile(self, NLLocfile):
-        self.w0 = w0
+        self.NLLocfile = NLLocfile
 
-    def getw0(self):
+    def getNLLocfile(self):
-        return self.w0
+        return self.NLLocfile
 
     def setrho(self, rho):
         self.rho = rho
@@ -93,11 +114,11 @@ class Magnitude(object):
     def getvp(self):
         return self.vp
 
-    def setdelta(self, delta):
+    def setpicks(self, picks):
-        self.delta = delta
+        self.picks = picks
 
-    def getdelta(self):
+    def getpicks(self):
-        return self.delta
+        return self.picks
 
     def getwapp(self):
         return self.wapp
@@ -108,11 +129,14 @@ class Magnitude(object):
     def getfc(self):
         return self.fc
 
-    def getMo(self):
+    def setinvdir(self, invdir):
-        return self.Mo
+        self.invdir = invdir
 
-    def getMw(self):
+    def getinvdir(self):
-        return self.Mw
+        return self.invdir
+
+    def getpicdic(self):
+        return self.picdic
 
     def calcwapp(self):
         self.wapp = None
@@ -120,9 +144,8 @@ class Magnitude(object):
     def calcsourcespec(self):
         self.sourcespek = None
 
-    def calcMoMw(self):
+    def run_calcMoMw(self):
-        self.Mo = None
+        self.pickdic = None
-        self.Mw = None
 
 class WApp(Magnitude):
     '''
@@ -177,48 +200,134 @@ class WApp(Magnitude):
 class M0Mw(Magnitude):
     '''
     Method to calculate seismic moment Mo and moment magnitude Mw.
-    Uses class w0fc for calculating plateau wo and corner frequency 
+    Requires results of class w0fc for calculating plateau w0 
-    fc of source spectrum, respectively.
+    and corner frequency fc of source spectrum, respectively. Uses 
+    subfunction calcMoMw.py. Returns modified dictionary of picks including 
+    Dc-value, corner frequency fc, seismic moment Mo and 
+    corresponding moment magntiude Mw.
     '''
 
-    def calcMoMw(self):
+    def run_calcMoMw(self):
 
-        stream = self.getwfstream()
+        picks = self.getpicks()
-        tr = stream[0]
+        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")
 
-        print("Calculating seismic moment Mo and moment magnitude Mw for station %s ..." \
+        for key in picks:
+             if picks[key]['P']['weight'] < 4:
+                 # select waveform
+                 selwf = zdat.select(station=key)
+                 # get hypocentral distance of station
+                 # from NLLoc-location file
+                 if len(key) > 4:
+                     Ppattern = '%s  ?    ?    ? P' % key
+                 elif len(key) == 4:
+                     Ppattern = '%s   ?    ?    ? P' % key
+                 elif len(key) < 4:
+                     Ppattern = '%s    ?    ?    ? P' % key
+                 nllocline = getPatternLine(nllocfile, Ppattern)
+                 delta = float(nllocline.split(None)[21])
+                 # call subfunction to estimate source spectrum
+                 # and to derive w0 and fc 
+                 [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, inv):
+    '''
+    Subfunction of run_calcMoMw to calculate individual
+    seismic moments and corresponding moment magnitudes.
+    '''
+
+    tr = wfstream[0]
+
+    print("calcMoMw: Calculating seismic moment Mo and moment magnitude Mw for station %s ..." \
            % tr.stats.station)
 
     # additional common parameters for calculating Mo
     rP = 2 / np.sqrt(15) # average radiation pattern of P waves (Aki & Richards, 1980)
     freesurf = 2.0       # free surface correction, assuming vertical incidence 
 
-        self.Mo = (self.getw0() * 4 * np.pi * self.getrho() * np.power(self.getvp(), 3) * \
+    Mo = w0 * 4 * np.pi * rho * np.power(vp, 3) * delta / (rP * freesurf) 
-                  self.getdelta()) / (rP * freesurf) 
 
-        self.Mw = 2/3 * np.log10(self.Mo) - 6
+    Mw = np.log10(Mo * 1e07) * 2 / 3 - 10.7 #after Hanks & Kanamori (1979), defined for [dyn*cm]!
-        print("MoMw: Calculated seismic moment Mo = %e Nm => Mw = %3.1f " % (self.Mo, self.Mw))
+
+    print("calcMoMw: Calculated seismic moment Mo = %e Nm => Mw = %3.1f " % (Mo, Mw))
+
+    return Mo, Mw
 
         
-class w0fc(Magnitude):
+
+def calcsourcespec(wfstream, onset, iplot, inventory):
     '''
-    Method to calculate the source spectrum and to derive from that the plateau
+    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!
+    source model. Has to be derived from instrument corrected displacement traces,
+    thus restitution and integration necessary!
     '''
-
-    def calcsourcespec(self):
     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
@@ -256,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()
@@ -287,6 +402,7 @@ class w0fc(Magnitude):
         raw_input()
         plt.close(f1)
 
+    return w0, fc
      
 
 def synthsourcespec(f, omega0, fcorner):

pylot/core/pick/autopick.py

@@ -15,10 +15,10 @@ from scipy import integrate
 from pylot.core.pick.Picker import AICPicker, PragPicker
 from pylot.core.pick.CharFuns import HOScf, AICcf, ARZcf, ARHcf, AR3Ccf
 from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker,\
-    getSNR, fmpicker, checkPonsets, wadaticheck, crossings_nonzero_all
+    getSNR, fmpicker, checkPonsets, wadaticheck
 from pylot.core.util.utils import getPatternLine
 from pylot.core.read.data import Data
-from pylot.core.analysis.magnitude import WApp, w0fc
+from pylot.core.analysis.magnitude import WApp
 
 def autopickevent(data, param):
     stations = []
@@ -138,8 +138,6 @@ def autopickstation(wfstream, pickparam, verbose=False):
     Sflag = 0
     Pmarker = []
     Ao = None     # Wood-Anderson peak-to-peak amplitude
-    w0 = None     # plateau of source spectrum
-    fc = None     # corner frequancy of source spectrum
 
     # split components
     zdat = wfstream.select(component="Z")
@@ -317,38 +315,6 @@ def autopickstation(wfstream, pickparam, verbose=False):
                 else:
                     FM = 'N'
 
-                ##############################################################
-                # get DC value and corner frequency (fc) of source spectrum
-                # from P pulse
-                # initialize Data object
-                data = Data()
-                z_copy = zdat.copy()
-                [corzdat, restflag] = data.restituteWFData(invdir, z_copy)
-                if restflag == 1:
-                    # integrate to displacement
-                    corintzdat = integrate.cumtrapz(corzdat[0], None, corzdat[0].stats.delta)
-                    z_copy[0].data = corintzdat
-                    # largest detectable period == window length
-                    # after P pulse for calculating source spectrum
-                    winzc = (1 / bpz2[0]) * z_copy[0].stats.sampling_rate
-                    impickP = mpickP * z_copy[0].stats.sampling_rate
-                    wfzc = z_copy[0].data[impickP : impickP + winzc]
-                    # 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:
-                        msg = "Something is wrong with the waveform, " \
-                              "no zero crossings derived!\nCannot " \
-                              "calculate source spectrum!"
-                        if verbose: print(msg)
-                    else:
-                        index = min([3, len(zc) - 1])
-                        calcwin = (zc[index] - zc[0]) * z_copy[0].stats.delta
-                        # calculate source spectrum and get w0 and fc
-                        specpara = w0fc(z_copy, mpickP, calcwin, iplot)
-                        w0 = specpara.getw0()
-                        fc = specpara.getfc()
-
                 msg = "autopickstation: P-weight: {0}, " \
                       "SNR: {1}, SNR[dB]: {2}, Polarity: {3}".format(Pweight,
                                                                      SNRP,
@@ -808,8 +774,7 @@ def autopickstation(wfstream, pickparam, verbose=False):
     # for P phase
     phase = 'P'
     phasepick = {'lpp': lpickP, 'epp': epickP, 'mpp': mpickP, 'spe': Perror,
-                 'snr': SNRP, 'snrdb': SNRPdB, 'weight': Pweight, 'fm': FM,
+                 'snr': SNRP, 'snrdb': SNRPdB, 'weight': Pweight, 'fm': FM}
-                 'w0': w0, 'fc': fc}
     picks = {phase: phasepick}
     # add P marker
     picks[phase]['marked'] = Pmarker