diff --git a/autoPyLoT.py b/autoPyLoT.py index 54003edb..eb205dbd 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -88,9 +88,10 @@ def autoPyLoT(inputfile): maxnumit = 3 # maximum number of iterations for re-picking else: locflag = 0 - print (" !!! ") - print ("!!No location routine available, autoPyLoT is running in non-location mode!!") - print (" !!! ") + print(" !!! ") + print("!!No location routine available, autoPyLoT is running in non-location mode!!") + 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 - # 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')) + nlloccounter = maxnumit + + # 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, 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 - - # 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')) + nlloccounter = maxnumit + # 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' diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index ff1f580e..d0e895c4 100644 --- a/pylot/core/analysis/magnitude.py +++ b/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 pylot.core.pick.utils import getsignalwin +from obspy.core import Stream, UTCDateTime +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 + :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): - self.w0 = w0 + def setNLLocfile(self, NLLocfile): + self.NLLocfile = NLLocfile - def getw0(self): - return self.w0 + def getNLLocfile(self): + 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): - self.delta = delta + def setpicks(self, picks): + self.picks = picks - def getdelta(self): - return self.delta + def getpicks(self): + 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): - return self.Mo + def setinvdir(self, invdir): + self.invdir = invdir - def getMw(self): - return self.Mw + def getinvdir(self): + 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): - self.Mo = None - self.Mw = None + def run_calcMoMw(self): + self.pickdic = None class WApp(Magnitude): ''' @@ -177,117 +200,210 @@ 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 - fc of source spectrum, respectively. + 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 + Dc-value, corner frequency fc, seismic moment Mo and + corresponding moment magntiude Mw. ''' - def calcMoMw(self): + def run_calcMoMw(self): - stream = self.getwfstream() - tr = stream[0] + 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") - print("Calculating seismic moment Mo and moment magnitude Mw for station %s ..." \ - % tr.stats.station) + 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()) - # 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 + 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 - self.Mo = (self.getw0() * 4 * np.pi * self.getrho() * np.power(self.getvp(), 3) * \ - self.getdelta()) / (rP * freesurf) + # 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 - self.Mw = 2/3 * np.log10(self.Mo) - 6 - print("MoMw: Calculated seismic moment Mo = %e Nm => Mw = %3.1f " % (self.Mo, self.Mw)) +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 + + Mo = w0 * 4 * np.pi * rho * np.power(vp, 3) * delta / (rP * freesurf) + + Mw = np.log10(Mo * 1e07) * 2 / 3 - 10.7 #after Hanks & Kanamori (1979), defined for [dyn*cm]! + + 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! ''' + print ("Calculating source spectrum ....") - def calcsourcespec(self): - print ("Calculating source spectrum ....") + fc = None + w0 = None + data = Data() + z_copy = wfstream.copy() - self.w0 = None # DC-value - self.fc = None # corner frequency - - stream = self.getwfstream() - tr = stream[0] + [corzdat, restflag] = data.restituteWFData(inventory, z_copy) + 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()) - xdat = tr.data[iwin] + # 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 - fny = tr.stats.sampling_rate / 2 - l = len(xdat) / tr.stats.sampling_rate - n = tr.stats.sampling_rate * l # number of fft bins after Bath - # find next power of 2 of data length - m = pow(2, np.ceil(np.log(len(xdat)) / np.log(2))) - N = int(np.power(m, 2)) - y = tr.stats.delta * np.fft.fft(xdat, N) - Y = abs(y[: N/2]) - L = (N - 1) / tr.stats.sampling_rate - f = np.arange(0, fny, 1/L) + # fft + fny = tr.stats.sampling_rate / 2 + l = len(xdat) / tr.stats.sampling_rate + n = tr.stats.sampling_rate * l # number of fft bins after Bath + # find next power of 2 of data length + m = pow(2, np.ceil(np.log(len(xdat)) / np.log(2))) + N = int(np.power(m, 2)) + y = tr.stats.delta * np.fft.fft(xdat, N) + Y = abs(y[: N/2]) + 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 - w0in = np.mean(YY[0:100]) - # initial guess of corner frequency - # where spectral level reached 50% of flat level - iin = np.where(YY >= 0.5 * w0in) - Fcin = F[iin[0][np.size(iin) - 1]] + # 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 + w0in = np.mean(YY[0:100]) + # initial guess of corner frequency + # where spectral level reached 50% of flat level + iin = np.where(YY >= 0.5 * w0in) + Fcin = F[iin[0][np.size(iin) - 1]] - # use of implicit scipy otimization function - fit = synthsourcespec(F, w0in, Fcin) - [optspecfit, pcov] = curve_fit(synthsourcespec, F, YY.real, [w0in, Fcin]) - w01 = optspecfit[0] - fc1 = optspecfit[1] - print ("w0fc: Determined w0-value: %e m/Hz, \n" - "Determined corner frequency: %f Hz" % (w01, fc1)) + # use of implicit scipy otimization function + fit = synthsourcespec(F, w0in, Fcin) + [optspecfit, pcov] = curve_fit(synthsourcespec, F, YY.real, [w0in, Fcin]) + w01 = optspecfit[0] + fc1 = optspecfit[1] + print ("w0fc: Determined w0-value: %e m/Hz, \n" + "Determined corner frequency: %f Hz" % (w01, fc1)) - # use of conventional fitting - [w02, fc2] = fitSourceModel(F, YY.real, Fcin, self.getiplot()) + # use of conventional fitting + [w02, fc2] = fitSourceModel(F, YY.real, Fcin, iplot) - # get w0 and fc as median - self.w0 = np.median([w01, w02]) - self.fc = np.median([fc1, fc2]) - print("w0fc: Using w0-value = %e m/Hz and fc = %f Hz" % (self.w0, self.fc)) + # get w0 and fc as median + w0 = np.median([w01, w02]) + fc = np.median([fc1, fc2]) + print("w0fc: Using w0-value = %e m/Hz and fc = %f Hz" % (w0, fc)) - if self.getiplot() > 1: - f1 = plt.figure() - plt.subplot(2,1,1) - # show displacement in mm - plt.plot(t, np.multiply(tr, 1000), 'k') + 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.xlabel('Time since %s' % tr.stats.starttime) - plt.ylabel('Displacement [mm]') + 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() - plt.show() - raw_input() - plt.close(f1) - + plt.show() + raw_input() + plt.close(f1) + return w0, fc + def synthsourcespec(f, omega0, fcorner): ''' diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 366a36e3..b6c1a11b 100755 --- a/pylot/core/pick/autopick.py +++ b/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, - 'w0': w0, 'fc': fc} + 'snr': SNRP, 'snrdb': SNRPdB, 'weight': Pweight, 'fm': FM} picks = {phase: phasepick} # add P marker picks[phase]['marked'] = Pmarker