From 30ee81a39d1d5238fb43dfa982ccf7ac3bb2578b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 23 Sep 2015 16:31:48 +0200 Subject: [PATCH 1/6] New class DCfc of object Magnitude for calculating source spectrum and to derive DC value and corner frequency. --- pylot/core/analysis/magnitude.py | 56 +++++++++++++++++++++++++++++--- 1 file changed, 51 insertions(+), 5 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 072984eb..2af92b37 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -25,7 +25,8 @@ class Magnitude(object): :type: float :param: pwin, pick window [To To+pwin] to get maximum - peak-to-peak amplitude + peak-to-peak amplitude (WApp) or to calculate + source spectrum (DCfc) :type: float :param: iplot, no. of figure window for plotting interims results @@ -40,6 +41,7 @@ class Magnitude(object): self.setpwin(pwin) self.setiplot(iplot) self.calcwapp() + self.calcsourcespec() def getwfstream(self): @@ -68,18 +70,23 @@ class Magnitude(object): def getwapp(self): return self.wapp + + def getw0(self): + return self.w0 + + def getfc(self): + return self.fc def calcwapp(self): self.wapp = None - def calcsourcespec(self): self.sourcespek = None class WApp(Magnitude): ''' Method to derive peak-to-peak amplitude as seen on a Wood-Anderson- - seismograph. Has to be derived from corrected traces! + seismograph. Has to be derived from instrument corrected traces! ''' def calcwapp(self): @@ -110,6 +117,7 @@ class WApp(Magnitude): iwin = getsignalwin(th, self.getTo(), self.getpwin()) self.wapp = np.max(sqH[iwin]) print ("Determined Wood-Anderson peak-to-peak amplitude: %f mm") % self.wapp + if self.getiplot() > 1: stream.plot() f = plt.figure(2) @@ -128,10 +136,48 @@ class WApp(Magnitude): class DCfc(Magnitude): ''' Method to calculate the source spectrum and to derive from that the plateau - (the so-called DC-value) and the corner frequency assuming Aki's omega-square - source model. Has to be derived from corrected traces! + (so-called DC-value) and the corner frequency assuming Aki's omega-square + source model. Has to be derived from instrument corrected displacement traces! ''' def calcsourcespec(self): print ("Calculating source spectrum ....") + self.w0 = None # DC-value + self.fc = None # corner frequency + + stream = self.getwfstream() + tr = stream[0] + + # 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] + + # fft + fny = tr.stats.sampling_rate / 2 + N = 1024 + 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) + + #if self.getiplot() > 1: + iplot=2 + if iplot > 1: + f1 = plt.figure(1) + plt.subplot(2,1,1) + plt.plot(t, np.multiply(tr, 1000), 'k') # show displacement in mm + plt.plot(t[iwin], np.multiply(xdat, 1000), 'g') # show displacement in mm + 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.title('Source Spectrum from P Pulse') + plt.xlabel('Frequency [Hz]') + plt.ylabel('Amplitude [m/Hz]') + plt.show() + raw_input() + plt.close(f1) From c790b4f35322e5f16330d6fbc9766aa264be9f42 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 23 Sep 2015 16:32:18 +0200 Subject: [PATCH 2/6] Implemented new class DCfc. --- pylot/core/pick/autopick.py | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 456023f5..10b4f7d5 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -11,12 +11,13 @@ function conglomerate utils. import matplotlib.pyplot as plt import numpy as np +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 + getSNR, fmpicker, checkPonsets, wadaticheck, crossings_nonzero_all from pylot.core.read.data import Data -from pylot.core.analysis.magnitude import WApp +from pylot.core.analysis.magnitude import WApp, DCfc def autopickevent(data, param): stations = [] @@ -309,6 +310,26 @@ def autopickstation(wfstream, pickparam): else: FM = 'N' + ############################################################## + # 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) + if restflag == 1: + # integrate to displacement + corintzdat = integrate.cumtrapz(corzdat[0], None, corzdat[0].stats.delta) + # class needs stream object => build it + z_copy = zdat.copy() + z_copy[0].data = corintzdat + # calculate source spectrum and get w0 and fc + calcwin = 1 / bpz2[0] # largest detectable period == window length + # around P pulse for calculating source spectrum + specpara = DCfc(z_copy, mpickP, calcwin, iplot) + w0 = specpara.getw0() + fc = specpara.getfc() + print 'autopickstation: P-weight: %d, SNR: %f, SNR[dB]: %f, ' \ 'Polarity: %s' % (Pweight, SNRP, SNRPdB, FM) Sflag = 1 From b391f5e08216513047e14e41d70f983967a159fa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 23 Sep 2015 17:02:27 +0200 Subject: [PATCH 3/6] Modified class DCfc: calculation of fft bins after Bath. --- pylot/core/analysis/magnitude.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 2af92b37..4391d16b 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -156,15 +156,17 @@ class DCfc(Magnitude): # fft fny = tr.stats.sampling_rate / 2 - N = 1024 + 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) + f = np.arange(0, fny, 1/L) - #if self.getiplot() > 1: - iplot=2 - if iplot > 1: + if self.getiplot() > 1: f1 = plt.figure(1) plt.subplot(2,1,1) plt.plot(t, np.multiply(tr, 1000), 'k') # show displacement in mm From 33164d4d1f116ef04d1038a246d7c4362403b473 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Fri, 25 Sep 2015 09:03:59 +0200 Subject: [PATCH 4/6] corrected for some python 3 compatibility issues; added a new delete picks button to the picking window --- icons.qrc | 1 + icons/delete.png | Bin 0 -> 5087 bytes icons_rc.py | 8 ++++---- pylot/core/read/inputs.py | 18 +++++++++--------- pylot/core/util/widgets.py | 11 ++++++++--- 5 files changed, 22 insertions(+), 16 deletions(-) create mode 100755 icons/delete.png diff --git a/icons.qrc b/icons.qrc index 0b815788..efd43234 100644 --- a/icons.qrc +++ b/icons.qrc @@ -3,6 +3,7 @@ icons/pylot.ico icons/pylot.png icons/printer.png + icons/delete.png icons/key_E.png icons/key_N.png icons/key_P.png diff --git a/icons/delete.png b/icons/delete.png new file mode 100755 index 0000000000000000000000000000000000000000..74c37cf9e2a9ac6f0c8f2ddecc302a491dee3f50 GIT binary patch literal 5087 zcmb7H2{@GN-yVjsN65as36acd&u;8XwnC+5%#0;tW@0P}V{1@I>J(Z?Orp67R6li}`WFD9%Y3<80~ z939Y}+$a3oPY?nE@c{g7n?WGH)(~%Bs;{ey362~F#p1~U1ZYecg&PF|nOek9u((hH z6$l^%g^xZH6MLfxGTkm5FFwV8&2?wb@Rr>hT@FzY8K|grZFa5fG`3T z3&eyGNf9P7W@_JYO}OvhKEu?2?;+GsGc~Jkgn+NBJ77Z&Cjfd-gf0$_G6eJtphzR6 zfw93(z?Dn_tSI3C90?%da61Npou5i?oUSU$V@GmN~M^< zVA0Xh&}bBt93BKi7#kbI;7Aw}smqPfjW|l8Vqg$6lw^U^cyBN zfE-CRQ&Z#413#C8rFez}1ydt_QII13G>&T_ECx$~A)s(r`1cstKQjOCoWBO+alc|I zk>SMe3&7)G1R^1fK%z!)IS@aH!L(@bIi zK>fyo|8FJ!!em9H{+;PtQA~(fQjnQij4qxKh>avt)y!wF)&Z_qJOSXA5zBqUZPnZs zW3}u1jSzz)Mn-Vu0fFRjz$*+J5gbAa;sym_V?#)Ue-rwX^}Eo23i~EDA~N7G0r#i0 z@0-yY8}xUf<{>yTDdO+wzw8PmVx!2B;s48Q#7~W3|Elp9%FikTR5I|R&VL8|G5V?O z`(pko{M$l)?Argo;f*)G9Zq8);3)7vN)P;-ioaRsI*DL=k(PdB%>zh;H8#%8cB?`r!gk5aTbWoUPbKq;b|@){(rD`#_R0nr zEak9N-d`@iQA>L&kuas62bph>1g2@V9E>khON3z+XO}nOBepag6|LIo2oiSgxY*s&&Z^(4CfcHsJQxS?XKufzzac0hg4k|GYjZDW zbO4+ck(8O=1O()bQ%lZ5uRFdad;sg<*(;M$DXa5+nJz;^=1&Y)RtBC9y{pB6uqpey z-=|1}_$25zEu$b7HS&j>8tb1d$&RI+-{eM`7JGRO`%}|;dem@DJd3* z#*y3LtaQ#Te<0Ie9!huuv#9V9WMc1MUGTxFL_uYsZV(wLl-{^MTF#fo!LmD_H9jwA zG|qh(Fa$Qba#$KwLP{#zRZL^bbvvFmKWl546Td$4`%NFwt(uymp}FQ;oZ|Hn zLwM`*P+hgkOvyft--1pjo91>9nW9BSX~sg*{50CvKzPDjtT`OWm9td9$O&3D+)8k+3TB9*#+Uu;NAJ&s@y=nbSQOs0| zhn5F$Zm?QAI_jJF+yX^8mJO@RFcuMY>}eOH08CL-%dHGzjOe3>GeMVb|m zqET910B zL?F;lj)EOCe9DpQ8>?9(vRaI!%>ulbp(GZ2h(%}dNS54rogvGDyP&9r2yf=5J$ab>Q<3zT+0neRy%UW#- zk!?451cBT>$s39M=g}_v-tkrL_Zc($U7k-)KIvcpo^Sgm3H~E}H29b*ycJvrAIkW= z48IQ5S*4RU(Zwx8qBI7gl_>!{^etL7a()^3EOv0w}7w433 z%B~js<^)f#*A$$00O^Q`phVQglMod#0gzs1qqV5x z_O8Tu5V*xM1!erZ+y<{0L}QLtA0TqJ!b+)n>o(Yvvf{fG+iC94 zE38E~#P22YILyYf(?=Hi%X8A{EuS||o|o&r_tIl0TjaqojJ2?xtXjXIS#2G&v?8a~ z3vB)(>sbZz?~PnmVb%DM>)(;0tAc!QmqG$l+DqH}x88>EZ`OiZt|qQpYgsYFHa@R)tEL61o=x;yLA=5X-$R#rl;o?56B|B$2lD+1E`&Z18hmcY|C-!fD!<* znSYCi%v-@jPA7ad-lpatB~+PFp&^|EuSg{WWydZ!a~`^O9cGC4)_1R?K3wEjM0C+= zXD=A<28xYq$0rqZZib3nDmR(gE7wC%B-%M;r=Kk}@jMnOeJo$1@k&ix(AcF<_JL2( z_m|Fls^}X(E-WcK%KK2=)2!r#L^%-GpM#>?y&Er&pPPICQVQinVU#iOFGGF$+!UQY zr0=MZ>KN_YF3st``7*SBVX05^I`g)l|MPdR9{Qi|rA(h0o>G#~+hTNUE4@V2 zBma;}p`msL;f(CV1ImGvxY#=}?r^p?piUeo%{wZlM9B4BZ+eCGLlKhig12;k>9UVP z8PNqZnW9y$srp$}J?6Rj`19Qn!O!^CXre}N{epID$TKvrkf{Gr%9KdcxJ5}*%d~qc zm06-oi2Y-5#q4P-qgPA=B{|6~i;E8ym5|EdE%WnAyj^XkmH0qR@I!zjLQ+8BG^kd~ z%cT(^6KnCCU+pj}XvpJkQnicG;k%_}I&V(oqz*ZMBzD_NOA|&v1QOH~gY|0l;*_%Q z?6;ac2ctZG>wMok@2qPvT1ueOoxLX5TuGm*Q!gwAa3K;4&3wx98To)Gv8f?0-THsh|ydxJ(GM#gYIG1?aZ6 zOOOcLQC6GK%yE4i4YO+VS6a$d+1Xn+YPE*#w=3Iv?`G|pUGnfR=|PN}Ek608vzj33 z4LTOdt15fKCD@VKvVxkx*9x@onKOz-Zi=ik%gXK$Fh=5>1Wndbt%pA+m+=XCP|{^s zDVk3yGv)G1T0J}OwAoX2!HshQ4=&=FyiRG>%PYMfU37yFH+t#@z8sFOb}fyc&wF|$ zUDeu{aXi$wtSY*Fx>QDQ3n!T&HDzn0z*JBffJiNEoS^`YHW5)RBh7d57lc#?<>Y7U zl(R|n;h{j6)tshH58l6;%FLRS)T+>yo1P3_JT~3RCI}ud4Bf-iaBg^LYlSGg6KbDZ zBlJ9qA1y9VpXnLv3DlRIp+_B9uu7Nu3cEY;0Q7cWDAC(i+I;!l+Zocp+umE!`>%f% zl<2c7mQfcTkyo^|E+}eH9p0v(r&&;hvGF*m4a-s(J87Z0BfUva+bI}q{v6aI04cjw zk#w#qmzUq^ar9$i^m12B;LPB$pa=e!c{46t8hm(a&}^u~IHGp8{quP*ZpFDN2zHs2oMAr$i_JzV{;%ts~jc)bwr6v7CLwDt4 zZbglKeD6N+=p62rICIfz?qJq?_MY*;6N>>kjlIIIz?Trd^XV3Xi*;L$U+$eBOMqT5 zs#?v`-!p9}TK`8uBwBzLQ}I>D zLV4PzJ%16ijid-&jS5LsEN65%73>lVS5)Jxa7Cp1pQud2Sey`46|*@;Xr9st&obJU zSk@nUt$+62+DaGb$(eRd$9E;ybF7YfpRU&`E}HB}WSd;ot;0a+w>2BR)2o zf4t7jyH6Z@qrZmawR~{orGDJdV#cRYv<2iwm}#^!V+%9CBKV20nEy3QwEcF6>~7;{ zqe5fYsCkEk*phNtbE6k{lU!v1zQA*t$!TUd4Xb1LoWCfA74JIj;3jt;A`TXj!;-}=nwxxuMrR9`OBdib62Lx?o zJsAOovQgARuvorD`XiyWofnYuGbim%?`uxlRDig5%2CX2MY=KQMfu@3p6OeZ4}W@{ z=;O0VXxQQ7>+$W5gZA<_e1xKfb{s&L7EFS7oX~m|)PG#(>>2gUJcP2r)K|wQ&;(}F z$Ib4G2_AcpmRAC^mM8pBn5o}U7r{F&9#nMm>tQ@cpRH2rf0U+K`QUi=?#{f7-S<@b z*ZEa^Jp0r+2DbT!=j}VSkKM@!YRVg09r#{#UxdK)5>AQ;y83zE0$N^z-W_u1Nvhgw z)%K7!U?(MyKCwTaG;S!`pHcaCpB<=Y!(rC8&M6LB+y;{V?aAQI8-Dz8&(_L5R`ytj ztaU==igFIuc*%@CeY1M0V&wLG+!rFa)Z+mtXDbMO27)y>SO$r#)*zwBPned5H`{ex zwx*A?M=p*973?eO1HQz+G9^ZpG&S!NQ24w!68L=Vh)Rj&$LF*6Kl4RKZQf9EK<8zh z=dFl8#;KP-`MA}B2n`(-C@r1kd)7FiyEVP6K)GaUH4kG>v^xK0A3{L*ZhV(##`0A` z;>scWE6TXk5=U97oF%4u(fnim#;~=~@Pf%wpR|>|x*hYB>#A#moMvQ!8J4#? z#i~{fA@nTy!gQqV{(=U_}p*i4f*0IbkR_Yd1Of4?SsK2)m>VVhNJ>7Hf_xm#$B3huM_>9pI7-l}5mf1%Ld++qM_ zx-H*J8F?+eV7;{UVx^Y*K!0D`0USQDEv6%Tj!$7I-aj^R+I@J0wdZKZgC&ix*4{dj=(3{Ty(aJcxUES*>$ip?_yEe6lAj zcIw5kc-j5)y)p)4!7-zZ^jJRc#y?oE-|wEUwtG7pDPZ^5d)C2q>3wLyuIt)yvX9c9 zqsh7IZVGlxGsl3s7RiwzNbdIez2GX7(1<@6!2wEpzVy?Th`)=buy3Xa>^~O_?eW&` zmgvbXZnni8ciVdktN6*KVlk+!d!5|AgTrh8xRB Date: Fri, 25 Sep 2015 15:06:59 +0200 Subject: [PATCH 5/6] some pick plotting issues solved; replotting main window if picks have been changed, otherwise just plot additional/deleted picks --- QtPyLoT.py | 27 +++++++++++++++++---------- pylot/core/util/widgets.py | 30 ++++++++++++++++++++++++++---- 2 files changed, 43 insertions(+), 14 deletions(-) diff --git a/QtPyLoT.py b/QtPyLoT.py index e4904c98..cae007c4 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -386,8 +386,8 @@ class MainWindow(QMainWindow): else: raise DatastructureError('not specified') return self.fnames - except DatastructureError, e: - print e + except DatastructureError as e: + print(e) props = PropertiesDlg(self) if props.exec_() == QDialog.Accepted: return self.getWFFnames() @@ -410,7 +410,7 @@ class MainWindow(QMainWindow): def saveData(self): def getSavePath(e): - print 'warning: {0}'.format(e) + print('warning: {0}'.format(e)) directory = os.path.join(self.getRoot(), self.getEventFileName()) file_filter = "QuakeML file (*.xml);;VELEST observation file format (*.cnv);;NonLinLoc observation file (*.obs)" fname = QFileDialog.getSaveFileName(self, 'Save event data ...', @@ -513,7 +513,7 @@ class MainWindow(QMainWindow): def plotWaveformData(self): zne_text = {'Z': 'vertical', 'N': 'north-south', 'E': 'east-west'} comp = self.getComponent() - title = 'overview: {0} components'.format(zne_text[comp]) + title = 'section: {0} components'.format(zne_text[comp]) wfst = self.getData().getWFData().select(component=comp) self.getPlotWidget().plotWFData(wfdata=wfst, title=title) self.draw() @@ -574,8 +574,8 @@ class MainWindow(QMainWindow): def getFilterOptions(self): try: return self.filteroptions[self.getSeismicPhase()] - except AttributeError, e: - print e + except AttributeError as e: + print(e) return FilterOptions(None, None, None) def getFilters(self): @@ -592,12 +592,12 @@ class MainWindow(QMainWindow): settings = QSettings() if settings.value("filterdefaults", None) is None and not self.getFilters(): - for key, value in FILTERDEFAULTS.iteritems(): + for key, value in FILTERDEFAULTS.items(): self.setFilterOptions(FilterOptions(**value), key) elif settings.value("filterdefaults", None) is not None: for key, value in settings.value("filterdefaults"): self.setFilterOptions(FilterOptions(**value), key) - except Exception, e: + except Exception as e: self.updateStatus('Error ...') emsg = QErrorMessage(self) emsg.showMessage('Error: {0}'.format(e)) @@ -636,8 +636,12 @@ class MainWindow(QMainWindow): if pickDlg.exec_(): self.setDirty(True) self.updateStatus('picks accepted ({0})'.format(station)) - self.addPicks(station, pickDlg.getPicks()) - self.drawPicks(station) + replot = self.addPicks(station, pickDlg.getPicks()) + if replot: + self.plotWaveformData() + self.drawPicks() + else: + self.drawPicks(station) else: self.updateStatus('picks discarded ({0})'.format(station)) @@ -667,6 +671,7 @@ class MainWindow(QMainWindow): def addPicks(self, station, picks): stat_picks = self.getPicksOnStation(station) + rval = False if not stat_picks: stat_picks = picks else: @@ -684,11 +689,13 @@ class MainWindow(QMainWindow): ret = msgBox.exec_() if ret == QMessageBox.Save: stat_picks = picks + rval = True elif ret == QMessageBox.Cancel: pass else: raise Exception('FATAL: Should never occur!') self.getPicks()[station] = stat_picks + return rval def updatePicks(self): evt = self.getData().getEvtData() diff --git a/pylot/core/util/widgets.py b/pylot/core/util/widgets.py index 1e7dd07a..6e680d81 100644 --- a/pylot/core/util/widgets.py +++ b/pylot/core/util/widgets.py @@ -240,7 +240,7 @@ class PickDlg(QDialog): home_icon = QIcon() home_icon.addPixmap(QPixmap(':/icons/zoom_0.png')) del_icon = QIcon() - del_icon.addPixmap(QPixmap(':/icon/delete.png')) + del_icon.addPixmap(QPixmap(':/icons/delete.png')) # create actions self.filterAction = createAction(parent=self, text='Filter', @@ -275,6 +275,8 @@ class PickDlg(QDialog): _dialtoolbar.addAction(self.zoomAction) _dialtoolbar.addSeparator() _dialtoolbar.addAction(self.resetZoomAction) + _dialtoolbar.addSeparator() + _dialtoolbar.addAction(self.resetPicksAction) # layout the innermost widget _innerlayout = QVBoxLayout() @@ -376,7 +378,7 @@ class PickDlg(QDialog): traceIDs = [] for channel in channels: channel = channel.upper() - for traceID, channelID in plotDict.iteritems(): + for traceID, channelID in plotDict.items(): if channelID[1].upper().endswith(channel): traceIDs.append(traceID) return traceIDs @@ -427,6 +429,13 @@ class PickDlg(QDialog): def getPicks(self): return self.picks + def resetPicks(self): + self.picks = {} + + def delPicks(self): + self.resetPicks() + self.resetPlot() + def setIniPick(self, gui_event): trace_number = round(gui_event.ydata) @@ -677,8 +686,21 @@ class PickDlg(QDialog): zoomx=self.getXLims(), zoomy=self.getYLims()) self.setPlotLabels() + self.drawPicks() self.draw() + def resetPlot(self): + self.updateCurrentLimits() + data = self.getWFData().copy() + title = self.getPlotWidget().getAxes().get_title() + self.getPlotWidget().plotWFData(wfdata=data, title=title, + zoomx=self.getXLims(), + zoomy=self.getYLims()) + self.setPlotLabels() + self.drawPicks() + self.draw() + + def setPlotLabels(self): # get channel labels @@ -716,7 +738,7 @@ class PickDlg(QDialog): else: # deal with something that should never happen scale_factor = 1 - print gui_event.button + print(gui_event.button) new_xlim = gui_event.xdata - \ scale_factor * (gui_event.xdata - self.getXLims()) @@ -747,7 +769,7 @@ class PickDlg(QDialog): def apply(self): picks = self.getPicks() for pick in picks: - print pick, picks[pick] + print(pick, picks[pick]) def accept(self): self.apply() From 4bd92565b9585b3f754ffe595deea1cadd3ea466 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Mon, 28 Sep 2015 10:57:48 +0200 Subject: [PATCH 6/6] reformat code (fix indentation and python 3.x issues) --- pylot/core/pick/Picker.py | 366 +++++++++++++++++++------------------- 1 file changed, 183 insertions(+), 183 deletions(-) diff --git a/pylot/core/pick/Picker.py b/pylot/core/pick/Picker.py index 61eade5e..f9ef82d4 100644 --- a/pylot/core/pick/Picker.py +++ b/pylot/core/pick/Picker.py @@ -147,7 +147,7 @@ class AICPicker(AutoPicking): def calcPick(self): - print 'AICPicker: Get initial onset time (pick) from AIC-CF ...' + print('AICPicker: Get initial onset time (pick) from AIC-CF ...') self.Pick = None self.slope = None @@ -155,7 +155,7 @@ class AICPicker(AutoPicking): #find NaN's nn = np.isnan(self.cf) if len(nn) > 1: - self.cf[nn] = 0 + self.cf[nn] = 0 #taper AIC-CF to get rid off side maxima tap = np.hanning(len(self.cf)) aic = tap * self.cf + max(abs(self.cf)) @@ -163,15 +163,15 @@ class AICPicker(AutoPicking): ismooth = int(round(self.Tsmooth / self.dt)) aicsmooth = np.zeros(len(aic)) if len(aic) < ismooth: - print 'AICPicker: Tsmooth larger than CF!' - return + print('AICPicker: Tsmooth larger than CF!') + return else: - for i in range(1, len(aic)): - if i > ismooth: - ii1 = i - ismooth - aicsmooth[i] = aicsmooth[i - 1] + (aic[i] - aic[ii1]) / ismooth - else: - aicsmooth[i] = np.mean(aic[1 : i]) + for i in range(1, len(aic)): + if i > ismooth: + ii1 = i - ismooth + aicsmooth[i] = aicsmooth[i - 1] + (aic[i] - aic[ii1]) / ismooth + else: + aicsmooth[i] = np.mean(aic[1 : i]) #remove offset offset = abs(min(aic) - min(aicsmooth)) aicsmooth = aicsmooth - offset @@ -180,7 +180,7 @@ class AICPicker(AutoPicking): #find NaN's nn = np.isnan(diffcf) if len(nn) > 1: - diffcf[nn] = 0 + diffcf[nn] = 0 #taper CF to get rid off side maxima tap = np.hanning(len(diffcf)) diffcf = tap * diffcf * max(abs(aicsmooth)) @@ -189,104 +189,104 @@ class AICPicker(AutoPicking): #find minimum in AIC-CF front of maximum lpickwindow = int(round(self.PickWindow / self.dt)) for i in range(icfmax - 1, max([icfmax - lpickwindow, 2]), -1): - if aicsmooth[i - 1] >= aicsmooth[i]: - self.Pick = self.Tcf[i] - break + if aicsmooth[i - 1] >= aicsmooth[i]: + self.Pick = self.Tcf[i] + break #if no minimum could be found: #search in 1st derivative of AIC-CF if self.Pick is None: - for i in range(icfmax -1, max([icfmax -lpickwindow, 2]), -1): - if diffcf[i -1] >= diffcf[i]: - self.Pick = self.Tcf[i] - break + for i in range(icfmax -1, max([icfmax -lpickwindow, 2]), -1): + if diffcf[i -1] >= diffcf[i]: + self.Pick = self.Tcf[i] + break # quality assessment using SNR and slope from CF if self.Pick is not None: - # get noise window - inoise = getnoisewin(self.Tcf, self.Pick, self.TSNR[0], self.TSNR[1]) - # check, if these are counts or m/s, important for slope estimation! - # this is quick and dirty, better solution? - if max(self.Data[0].data < 1e-3): - self.Data[0].data = self.Data[0].data * 1000000 - # get signal window - isignal = getsignalwin(self.Tcf, self.Pick, self.TSNR[2]) - # calculate SNR from CF - self.SNR = max(abs(aic[isignal] - np.mean(aic[isignal]))) / max(abs(aic[inoise] \ - - np.mean(aic[inoise]))) - # calculate slope from CF after initial pick - # get slope window - tslope = self.TSNR[3] #slope determination window - islope = np.where((self.Tcf <= min([self.Pick + tslope, len(self.Data[0].data)])) \ - & (self.Tcf >= self.Pick)) - # find maximum within slope determination window - # 'cause slope should be calculated up to first local minimum only! - imax = np.argmax(self.Data[0].data[islope]) - if imax == 0: - print 'AICPicker: Maximum for slope determination right at the beginning of the window!' - print 'Choose longer slope determination window!' - if self.iplot > 1: - p = plt.figure(self.iplot) - x = self.Data[0].data - p1, = plt.plot(self.Tcf, x / max(x), 'k') - p2, = plt.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r') - plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF']) - plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) - plt.yticks([]) - plt.title(self.Data[0].stats.station) - plt.show() - raw_input() - plt.close(p) - return - islope = islope[0][0 :imax] - dataslope = self.Data[0].data[islope] - # calculate slope as polynomal fit of order 1 - xslope = np.arange(0, len(dataslope), 1) - P = np.polyfit(xslope, dataslope, 1) - datafit = np.polyval(P, xslope) - if datafit[0] >= datafit[len(datafit) - 1]: - print 'AICPicker: Negative slope, bad onset skipped!' - return - self.slope = 1 / tslope * (datafit[len(dataslope) - 1] - datafit[0]) + # get noise window + inoise = getnoisewin(self.Tcf, self.Pick, self.TSNR[0], self.TSNR[1]) + # check, if these are counts or m/s, important for slope estimation! + # this is quick and dirty, better solution? + if max(self.Data[0].data < 1e-3): + self.Data[0].data = self.Data[0].data * 1000000 + # get signal window + isignal = getsignalwin(self.Tcf, self.Pick, self.TSNR[2]) + # calculate SNR from CF + self.SNR = max(abs(aic[isignal] - np.mean(aic[isignal]))) / \ + max(abs(aic[inoise] - np.mean(aic[inoise]))) + # calculate slope from CF after initial pick + # get slope window + tslope = self.TSNR[3] #slope determination window + islope = np.where((self.Tcf <= min([self.Pick + tslope, len(self.Data[0].data)])) \ + & (self.Tcf >= self.Pick)) + # find maximum within slope determination window + # 'cause slope should be calculated up to first local minimum only! + imax = np.argmax(self.Data[0].data[islope]) + if imax == 0: + print('AICPicker: Maximum for slope determination right at the beginning of the window!') + print('Choose longer slope determination window!') + if self.iplot > 1: + p = plt.figure(self.iplot) + x = self.Data[0].data + p1, = plt.plot(self.Tcf, x / max(x), 'k') + p2, = plt.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r') + plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF']) + plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) + plt.yticks([]) + plt.title(self.Data[0].stats.station) + plt.show() + raw_input() + plt.close(p) + return + islope = islope[0][0 :imax] + dataslope = self.Data[0].data[islope] + # calculate slope as polynomal fit of order 1 + xslope = np.arange(0, len(dataslope), 1) + P = np.polyfit(xslope, dataslope, 1) + datafit = np.polyval(P, xslope) + if datafit[0] >= datafit[len(datafit) - 1]: + print('AICPicker: Negative slope, bad onset skipped!') + return + self.slope = 1 / tslope * (datafit[len(dataslope) - 1] - datafit[0]) else: - self.SNR = None - self.slope = None + self.SNR = None + self.slope = None if self.iplot > 1: - p = plt.figure(self.iplot) - x = self.Data[0].data - p1, = plt.plot(self.Tcf, x / max(x), 'k') - p2, = plt.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r') - if self.Pick is not None: - p3, = plt.plot([self.Pick, self.Pick], [-0.1 , 0.5], 'b', linewidth=2) - plt.legend([p1, p2, p3], ['(HOS-/AR-) Data', 'Smoothed AIC-CF', 'AIC-Pick']) - else: - plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF']) - plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) - plt.yticks([]) - plt.title(self.Data[0].stats.station) + p = plt.figure(self.iplot) + x = self.Data[0].data + p1, = plt.plot(self.Tcf, x / max(x), 'k') + p2, = plt.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r') + if self.Pick is not None: + p3, = plt.plot([self.Pick, self.Pick], [-0.1 , 0.5], 'b', linewidth=2) + plt.legend([p1, p2, p3], ['(HOS-/AR-) Data', 'Smoothed AIC-CF', 'AIC-Pick']) + else: + plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF']) + plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) + plt.yticks([]) + plt.title(self.Data[0].stats.station) - if self.Pick is not None: - plt.figure(self.iplot + 1) - p11, = plt.plot(self.Tcf, x, 'k') - p12, = plt.plot(self.Tcf[inoise], self.Data[0].data[inoise]) - p13, = plt.plot(self.Tcf[isignal], self.Data[0].data[isignal], 'r') - p14, = plt.plot(self.Tcf[islope], dataslope, 'g--') - p15, = plt.plot(self.Tcf[islope], datafit, 'g', linewidth=2) - plt.legend([p11, p12, p13, p14, p15], ['Data', 'Noise Window', 'Signal Window', 'Slope Window', 'Slope'], \ - loc='best') - plt.title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station, \ - self.SNR, self.slope)) - plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) - plt.ylabel('Counts') - plt.yticks([]) + if self.Pick is not None: + plt.figure(self.iplot + 1) + p11, = plt.plot(self.Tcf, x, 'k') + p12, = plt.plot(self.Tcf[inoise], self.Data[0].data[inoise]) + p13, = plt.plot(self.Tcf[isignal], self.Data[0].data[isignal], 'r') + p14, = plt.plot(self.Tcf[islope], dataslope, 'g--') + p15, = plt.plot(self.Tcf[islope], datafit, 'g', linewidth=2) + plt.legend([p11, p12, p13, p14, p15], ['Data', 'Noise Window', 'Signal Window', 'Slope Window', 'Slope'], + loc='best') + plt.title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station, + self.SNR, self.slope)) + plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) + plt.ylabel('Counts') + plt.yticks([]) - plt.show() - raw_input() - plt.close(p) + plt.show() + raw_input() + plt.close(p) if self.Pick == None: - print 'AICPicker: Could not find minimum, picking window too short?' + print('AICPicker: Could not find minimum, picking window too short?') class PragPicker(AutoPicking): @@ -297,102 +297,102 @@ class PragPicker(AutoPicking): def calcPick(self): if self.getpick1() is not None: - print 'PragPicker: Get most likely pick from HOS- or AR-CF using pragmatic picking algorithm ...' + print('PragPicker: Get most likely pick from HOS- or AR-CF using pragmatic picking algorithm ...') - self.Pick = None - self.SNR = None - self.slope = None - pickflag = 0 - #smooth CF - ismooth = int(round(self.Tsmooth / self.dt)) - cfsmooth = np.zeros(len(self.cf)) - if len(self.cf) < ismooth: - print 'PragPicker: Tsmooth larger than CF!' - return - else: - for i in range(1, len(self.cf)): - if i > ismooth: - ii1 = i - ismooth; - cfsmooth[i] = cfsmooth[i - 1] + (self.cf[i] - self.cf[ii1]) / ismooth - else: - cfsmooth[i] = np.mean(self.cf[1 : i]) + self.Pick = None + self.SNR = None + self.slope = None + pickflag = 0 + #smooth CF + ismooth = int(round(self.Tsmooth / self.dt)) + cfsmooth = np.zeros(len(self.cf)) + if len(self.cf) < ismooth: + print('PragPicker: Tsmooth larger than CF!') + return + else: + for i in range(1, len(self.cf)): + if i > ismooth: + ii1 = i - ismooth; + cfsmooth[i] = cfsmooth[i - 1] + (self.cf[i] - self.cf[ii1]) / ismooth + else: + cfsmooth[i] = np.mean(self.cf[1 : i]) - #select picking window - #which is centered around tpick1 - ipick = np.where((self.Tcf >= self.getpick1() - self.PickWindow / 2) \ - & (self.Tcf <= self.getpick1() + self.PickWindow / 2)) - cfipick = self.cf[ipick] - np.mean(self.cf[ipick]) - Tcfpick = self.Tcf[ipick] - cfsmoothipick = cfsmooth[ipick]- np.mean(self.cf[ipick]) - ipick1 = np.argmin(abs(self.Tcf - self.getpick1())) - cfpick1 = 2 * self.cf[ipick1] + #select picking window + #which is centered around tpick1 + ipick = np.where((self.Tcf >= self.getpick1() - self.PickWindow / 2) \ + & (self.Tcf <= self.getpick1() + self.PickWindow / 2)) + cfipick = self.cf[ipick] - np.mean(self.cf[ipick]) + Tcfpick = self.Tcf[ipick] + cfsmoothipick = cfsmooth[ipick]- np.mean(self.cf[ipick]) + ipick1 = np.argmin(abs(self.Tcf - self.getpick1())) + cfpick1 = 2 * self.cf[ipick1] - #check trend of CF, i.e. differences of CF and adjust aus regarding this trend - #prominent trend: decrease aus - #flat: use given aus - cfdiff = np.diff(cfipick); - i0diff = np.where(cfdiff > 0) - cfdiff = cfdiff[i0diff] - minaus = min(cfdiff * (1 + self.aus)); - aus1 = max([minaus, self.aus]); + #check trend of CF, i.e. differences of CF and adjust aus regarding this trend + #prominent trend: decrease aus + #flat: use given aus + cfdiff = np.diff(cfipick) + i0diff = np.where(cfdiff > 0) + cfdiff = cfdiff[i0diff] + minaus = min(cfdiff * (1 + self.aus)) + aus1 = max([minaus, self.aus]) - #at first we look to the right until the end of the pick window is reached - flagpick_r = 0 - flagpick_l = 0 - cfpick_r = 0 - cfpick_l = 0 - lpickwindow = int(round(self.PickWindow / self.dt)) - for i in range(max(np.insert(ipick, 0, 2)), min([ipick1 + lpickwindow + 1, len(self.cf) - 1])): - if self.cf[i + 1] > self.cf[i] and self.cf[i - 1] >= self.cf[i]: - if cfsmooth[i - 1] * (1 + aus1) >= cfsmooth[i]: - if cfpick1 >= self.cf[i]: - pick_r = self.Tcf[i] - self.Pick = pick_r - flagpick_l = 1 - cfpick_r = self.cf[i] - break + #at first we look to the right until the end of the pick window is reached + flagpick_r = 0 + flagpick_l = 0 + cfpick_r = 0 + cfpick_l = 0 + lpickwindow = int(round(self.PickWindow / self.dt)) + for i in range(max(np.insert(ipick, 0, 2)), min([ipick1 + lpickwindow + 1, len(self.cf) - 1])): + if self.cf[i + 1] > self.cf[i] and self.cf[i - 1] >= self.cf[i]: + if cfsmooth[i - 1] * (1 + aus1) >= cfsmooth[i]: + if cfpick1 >= self.cf[i]: + pick_r = self.Tcf[i] + self.Pick = pick_r + flagpick_l = 1 + cfpick_r = self.cf[i] + break - # now we look to the left - for i in range(ipick1, max([ipick1 - lpickwindow + 1, 2]), -1): - if self.cf[i + 1] > self.cf[i] and self.cf[i - 1] >= self.cf[i]: - if cfsmooth[i - 1] * (1 + aus1) >= cfsmooth[i]: - if cfpick1 >= self.cf[i]: - pick_l = self.Tcf[i] - self.Pick = pick_l - flagpick_r = 1 - cfpick_l = self.cf[i] - break + # now we look to the left + for i in range(ipick1, max([ipick1 - lpickwindow + 1, 2]), -1): + if self.cf[i + 1] > self.cf[i] and self.cf[i - 1] >= self.cf[i]: + if cfsmooth[i - 1] * (1 + aus1) >= cfsmooth[i]: + if cfpick1 >= self.cf[i]: + pick_l = self.Tcf[i] + self.Pick = pick_l + flagpick_r = 1 + cfpick_l = self.cf[i] + break - # now decide which pick: left or right? - if flagpick_l > 0 and flagpick_r > 0 and cfpick_l <= cfpick_r: - self.Pick = pick_l - pickflag = 1 - elif flagpick_l > 0 and flagpick_r > 0 and cfpick_l >= cfpick_r: - self.Pick = pick_r - pickflag = 1 - elif flagpick_l == 0 and flagpick_r > 0 and cfpick_l >= cfpick_r: + # now decide which pick: left or right? + if flagpick_l > 0 and flagpick_r > 0 and cfpick_l <= cfpick_r: self.Pick = pick_l pickflag = 1 - else: - print 'PragPicker: Could not find reliable onset!' - self.Pick = None - pickflag = 0 + elif flagpick_l > 0 and flagpick_r > 0 and cfpick_l >= cfpick_r: + self.Pick = pick_r + pickflag = 1 + elif flagpick_l == 0 and flagpick_r > 0 and cfpick_l >= cfpick_r: + self.Pick = pick_l + pickflag = 1 + else: + print('PragPicker: Could not find reliable onset!') + self.Pick = None + pickflag = 0 - if self.getiplot() > 1: - p = plt.figure(self.getiplot()) - p1, = plt.plot(Tcfpick,cfipick, 'k') - p2, = plt.plot(Tcfpick,cfsmoothipick, 'r') - if pickflag > 0: - p3, = plt.plot([self.Pick, self.Pick], [min(cfipick), max(cfipick)], 'b', linewidth=2) - plt.legend([p1, p2, p3], ['CF', 'Smoothed CF', 'Pick']) - plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) - plt.yticks([]) - plt.title(self.Data[0].stats.station) - plt.show() - raw_input() - plt.close(p) + if self.getiplot() > 1: + p = plt.figure(self.getiplot()) + p1, = plt.plot(Tcfpick,cfipick, 'k') + p2, = plt.plot(Tcfpick,cfsmoothipick, 'r') + if pickflag > 0: + p3, = plt.plot([self.Pick, self.Pick], [min(cfipick), max(cfipick)], 'b', linewidth=2) + plt.legend([p1, p2, p3], ['CF', 'Smoothed CF', 'Pick']) + plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) + plt.yticks([]) + plt.title(self.Data[0].stats.station) + plt.show() + raw_input() + plt.close(p) else: - print 'PragPicker: No initial onset time given! Check input!' - self.Pick = None - return + print('PragPicker: No initial onset time given! Check input!') + self.Pick = None + return