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

pylot/core/analysis/magnitude.py

@@ -9,6 +9,7 @@ import matplotlib.pyplot as plt
 import numpy as np
 from obspy.core import Stream
 from pylot.core.pick.utils import getsignalwin
+from scipy.optimize import curve_fit
 
 class Magnitude(object):
     '''
@@ -166,20 +167,68 @@ class DCfc(Magnitude):
         L = (N - 1) / tr.stats.sampling_rate
         f = np.arange(0, fny, 1/L)
 
+        # remove zero-frequency and frequencies above 
+        # corner frequency of seismometer (assumed
+        # to be 100 Hz)
+        fi = np.where((f >= 1) & (f < 100))
+        F = f[fi]
+        YY = Y[fi]
+        # get plateau (DC value) and corner frequency
+        # initial guess of plateau
+        DCin = np.mean(YY[0:100])
+        # initial guess of corner frequency
+        # where spectral level reached 50% of flat level
+        iin = np.where(YY >= 0.5 * DCin)
+        Fcin = F[iin[0][np.size(iin) - 1]]
+        fit = synthsourcespec(F, DCin, Fcin)
+        [optspecfit, pcov] = curve_fit(synthsourcespec, F, YY.real, [DCin, Fcin])
+        self.w0 = optspecfit[0]
+        self.fc = optspecfit[1]
+        print ("DCfc: Determined DC-value: %e m/Hz, \n" \
+               "Determined corner frequency: %f Hz" % (self.w0, self.fc)) 
+        
+        
         if self.getiplot() > 1:
-            f1 = plt.figure(1)
+            f1 = plt.figure()
             plt.subplot(2,1,1)
-            plt.plot(t, np.multiply(tr, 1000), 'k') # show displacement in mm
+            # show displacement in mm
-            plt.plot(t[iwin], np.multiply(xdat, 1000), 'g') # show displacement in mm
+            plt.plot(t, np.multiply(tr, 1000), 'k') 
+            plt.plot(t[iwin], np.multiply(xdat, 1000), 'g') 
             plt.title('Seismogram and P pulse, station %s' % tr.stats.station)
             plt.xlabel('Time since %s' % tr.stats.starttime)
             plt.ylabel('Displacement [mm]')
 
             plt.subplot(2,1,2)
-            plt.semilogy(f, Y.real)
+            plt.loglog(f, Y.real, 'k')
-            plt.title('Source Spectrum from P Pulse')
+            plt.loglog(F, YY.real)
+            plt.loglog(F, fit, 'g')
+            plt.title('Source Spectrum from P Pulse, DC=%e m/Hz, fc=%4.1f Hz' \
+                       % (self.w0, self.fc))
             plt.xlabel('Frequency [Hz]')
             plt.ylabel('Amplitude [m/Hz]')
+            plt.grid()
             plt.show()
             raw_input()
             plt.close(f1)
+
+
+def synthsourcespec(f, omega0, fcorner):
+    '''
+    Calculates synthetic source spectrum from given plateau and corner 
+    frequency assuming Akis omega-square model.  
+
+    :param: f, frequencies
+    :type:  array
+
+    :param: omega0, DC-value (plateau) of source spectrum
+    :type:  float
+
+    :param: fcorner, corner frequency of source spectrum
+    :type:  float
+    '''
+        
+    #ssp = omega0 / (pow(2, (1 + f / fcorner)))
+    ssp = omega0 / (1 + pow(2, (f / fcorner)))
+
+    return ssp
+

pylot/core/pick/autopick.py

@@ -313,7 +313,6 @@ def autopickstation(wfstream, pickparam):
                 ##############################################################
                 # get DC value (w0) and corner frequency (fc) of source spectrum
                 # from P pulse
-                # restitute streams
                 # initialize Data object
                 data = Data()
                 [corzdat, restflag] = data.restituteWFData(invdir, zdat)
@@ -323,33 +322,45 @@ def autopickstation(wfstream, pickparam):
                         # class needs stream object => build it
                         z_copy = zdat.copy()
                         z_copy[0].data = corintzdat
-                	# calculate source spectrum and get w0 and fc
+                        # largest detectable period == window length
-                        calcwin = 1 / bpz2[0] # largest detectable period == window length 
+                        # after P pulse for calculating source spectrum
-                                              # around P pulse for calculating source spectrum
+                        winzc = (1 / bpz2[0]) * z_copy[0].stats.sampling_rate  
-                	specpara = DCfc(z_copy, mpickP, calcwin, iplot)
+                        impickP = mpickP * z_copy[0].stats.sampling_rate
-                        w0 = specpara.getw0()
+                        wfzc = z_copy[0].data[impickP : impickP + winzc]
-                        fc = specpara.getfc()
+                        # calculate spectrum using only first cycles of
+                        # waveform after P onset!
+                        zc = crossings_nonzero_all(wfzc)
+                        if np.size(zc) == 0:
+                            print ("Something is wrong with the waveform, " \
+                                   "no zero crossings derived!")
+                            print ("Cannot calculate source spectrum!")
+                        else:
+                            calcwin = (zc[3] - zc[0]) * z_copy[0].stats.delta
+                	    # calculate source spectrum and get w0 and fc
+                	    specpara = DCfc(z_copy, mpickP, calcwin, iplot)
+                            w0 = specpara.getw0()
+                            fc = specpara.getfc()
 
-                print 'autopickstation: P-weight: %d, SNR: %f, SNR[dB]: %f, ' \
+                print ("autopickstation: P-weight: %d, SNR: %f, SNR[dB]: %f, " \
-                      'Polarity: %s' % (Pweight, SNRP, SNRPdB, FM)
+                       "Polarity: %s" % (Pweight, SNRP, SNRPdB, FM))
                 Sflag = 1
 
         else:
-            print 'Bad initial (AIC) P-pick, skipping this onset!'
+            print ("Bad initial (AIC) P-pick, skipping this onset!")
             print 'AIC-SNR=', aicpick.getSNR(), 'AIC-Slope=', aicpick.getSlope(), 'counts/s'
             print '(min. AIC-SNR=', minAICPSNR, ', min. AIC-Slope=', minAICPslope, 'counts/s)'
             Sflag = 0
 
     else:
-        print 'autopickstation: No vertical component data available!, ' \
+        print ("autopickstation: No vertical component data available!, " \
-              'Skipping station!'
+               "Skipping station!")
 
     if edat is not None and ndat is not None and len(edat) > 0 and len(
             ndat) > 0 and Pweight < 4:
-        print 'Go on picking S onset ...'
+        print ("Go on picking S onset ...")
-        print '##################################################'
+        print ("##################################################")
-        print 'Working on S onset of station %s' % edat[0].stats.station
+        print ("Working on S onset of station %s" % edat[0].stats.station)
-        print 'Filtering horizontal traces ...'
+        print ("Filtering horizontal traces ...")
 
         # determine time window for calculating CF after P onset
         cuttimesh = [round(max([mpickP + sstart, 0])),