From a854982400e8408015945aedb31ecce521c5e948 Mon Sep 17 00:00:00 2001 From: xsong Date: Mon, 2 May 2016 00:51:51 -0400 Subject: [PATCH 1/5] Add xcorr_full_polarity function to support mccc without changing polarity --- src/pysmo/aimbat/xcorr.f90 | 43 ++++++++++++++++++++++++++++++++++++++ src/pysmo/aimbat/xcorr.py | 30 ++++++++++++++++++++++++++ 2 files changed, 73 insertions(+) mode change 100644 => 100755 src/pysmo/aimbat/xcorr.py diff --git a/src/pysmo/aimbat/xcorr.f90 b/src/pysmo/aimbat/xcorr.f90 index fc377e7b..326687d1 100644 --- a/src/pysmo/aimbat/xcorr.f90 +++ b/src/pysmo/aimbat/xcorr.f90 @@ -84,6 +84,49 @@ end subroutine xcorr_full !!!---------------------------------------------------------------------!!! +!!!---------------------------------------------------------------------!!! +subroutine xcorr_full_polarity(x, y, n, shift, delay, ccmax, ccpol) +! Cross-correlation of array x and y. +! It does not correct the polarity, meaning it does only calculate the maximum value of cross correlation, not including minimum, which is different from xcorr_full +! Return time shift, correlation coefficient, and polarity at maximum correlation. +implicit none +real(8) :: crosscorr +integer :: n, shift, delay, ccpol +real(8), dimension(0:n-1) :: x, y +real(8) :: cc, ccmax, ccmin +integer :: k, kmin +!f2py intent(in) :: x, y, n, shift +!f2py intent(out) :: delay, ccmax, ccpol +!f2py intent(hide):: c, k, kmin, ccmin + +shift = 1 +ccmax = 0 +ccmin = 0 +kmin = 0 +do k = -n+1,n-1,shift + cc = crosscorr(x,y,n,k) + if (cc.gt.ccmax) then + ccmax = cc + delay = k + endif + if (cc.lt.ccmin) then + ccmin = cc + kmin = k + endif +enddo +! if (ccmax.gt.-ccmin) then +! ccpol = 1 +! else +! ccmax = -ccmin +! delay = kmin +! ccpol = -1 +ccpol = 1 +ccmax = ccmax/sqrt( dot_product(x,x) * dot_product(y,y) ) + +end subroutine xcorr_full_polarity +!!!---------------------------------------------------------------------!!! + + !!!---------------------------------------------------------------------!!! subroutine xcorr_fast(x, y, n, shift, delay, ccmax, ccpol) ! Fast cross-correlation using 1 level of coarse shift. diff --git a/src/pysmo/aimbat/xcorr.py b/src/pysmo/aimbat/xcorr.py old mode 100644 new mode 100755 index f0a9656b..7289e01a --- a/src/pysmo/aimbat/xcorr.py +++ b/src/pysmo/aimbat/xcorr.py @@ -50,6 +50,28 @@ def _xcorr(x, y, cmode='full'): delay -= (len(y)-1)/2 return delay, ccmax, ccpol +def _xcorr_polarity(x, y, cmode='full'): + """ + Cross-correlation of two 1-D arrays using 'full' or 'same' mode. + c[k] = sum_i x[i]*y[i+k] + Full mode: indices of array c --> j=0:nx+ny-2 <-- k=ny-1:-nx+1:-1 (j=ny-1-k) + It does not correct the polarity, meaning it does only calculate the maximum + value of cross correlaton, not including minimum, which is different from + _xcorr + Return time shift of y relative to x at maximum correlation and polarity. + """ + cc = correlate(x, y, cmode) + imax = argmax(cc) + imin = argmin(cc) + ccmax = cc[imax] + ccmin = cc[imin] + ccpol = 1 + delay = len(y)-1-imax + ccmax = ccmax / sqrt(dot(x,x)*dot(y,y)) + if cmode == 'same': + delay -= (len(y)-1)/2 + return delay, ccmax, ccpol + def xcorr_full(x, y, shift=1): """ Cross-correlation of two 1-D arrays using 'full' mode. @@ -57,6 +79,14 @@ def xcorr_full(x, y, shift=1): """ return _xcorr(x, y, 'full') +def xcorr_full_polarity(x, y, shift=1): + """ + Cross-correlation of two 1-D arrays using 'full' mode. + Argument shift=1 is here only in order to make the same number of arguments for all xcorr functions. + Not correct the polarity + """ + return _xcorr_polarity(x, y, 'full') + def xcorr_same(x, y): """ Cross-correlation of two 1-D arrays using 'same' mode. From a8bb7b02456edca022abde85836a45bbf4f33282 Mon Sep 17 00:00:00 2001 From: xsong Date: Mon, 2 May 2016 01:01:57 -0400 Subject: [PATCH 2/5] Show mccc results for every trace in the graphic window Show different colors for different mccc thresholds Write inversion matrix to output file --- src/pysmo/aimbat/algmccc.py | 20 +++++++++++--------- src/pysmo/aimbat/pickphase.py | 26 +++++++++++++++++++++++--- src/pysmo/aimbat/qualctrl.py | 4 ++-- src/pysmo/aimbat/ttconfig.py | 5 +++++ src/pysmo/aimbat/ttdefaults.conf | 4 ++-- 5 files changed, 43 insertions(+), 16 deletions(-) diff --git a/src/pysmo/aimbat/algmccc.py b/src/pysmo/aimbat/algmccc.py index 5761eefb..3f8972aa 100755 --- a/src/pysmo/aimbat/algmccc.py +++ b/src/pysmo/aimbat/algmccc.py @@ -26,7 +26,7 @@ from time import strftime, tzname from optparse import OptionParser from ttconfig import MCConfig -from sacpickle import loadData, saveData, SacDataHdrs, taper, taperWindow, windowIndex, windowData +from sacpickle import loadData, saveData, SacDataHdrs, taper, taperWindow, windowIndex, windowData from qualsort import initQual, seleSeis def getOptions(): @@ -279,7 +279,7 @@ def WriteFileWithDelay(mcpara, solist, solution, outvar, outcc, t0_times, delay_ line1 = 'station, mccc delay, std, cc coeff, cc std, pol , t0_times , delay_times\n' ofile.write( line0 ) ofile.write( line1 ) - fmt = ' {0:<9s} {1:9.4f} {2:9.4f} {3:>9.4f} {4:>9.4f} {5:4d} {6: 10000 or abs(dt) < threshd[0]: + colors[j][k] = colsel + elif threshd[0] <= abs(dt) < threshd[1]: + colors[j][k] = colthd[0] + elif threshd[1] <= abs(dt) < threshd[2]: + colors[j][k] = colthd[1] + if threshd[2] <= abs(dt) < 10000: + colors[j][k] = colthd[2] else: colors[j][k] = coldel + # plot pps = [] for j in range(2): diff --git a/src/pysmo/aimbat/qualctrl.py b/src/pysmo/aimbat/qualctrl.py index 0977308c..a70f1226 100755 --- a/src/pysmo/aimbat/qualctrl.py +++ b/src/pysmo/aimbat/qualctrl.py @@ -167,7 +167,7 @@ def setLabels(self): trans = transforms.blended_transform_factory(axpp.transAxes, axpp.transData) font = FontProperties() font.set_family('monospace') - axpp.text(1.025, 0, ' '*8+'qual= CCC/SNR/COH', transform=trans, va='center', + axpp.text(1.025, 0, ' '*8+'qual= MCCC/CCC/SNR/COH', transform=trans, va='center', color='k', fontproperties=font) def plotStack(self): @@ -1093,7 +1093,7 @@ def mccc(self, event): mcpara.evline = evline mcpara.mcname = mcname mcpara.kevnm = gsac.kevnm - solution, solist_LonLat, delay_times = mccc(gsac, mcpara) + solution, solist_LonLat, delay_times, itmean = mccc(gsac, mcpara) self.gsac.solist_LonLat = solist_LonLat self.gsac.delay_times = delay_times diff --git a/src/pysmo/aimbat/ttconfig.py b/src/pysmo/aimbat/ttconfig.py index 2f3801aa..9b2a06c2 100644 --- a/src/pysmo/aimbat/ttconfig.py +++ b/src/pysmo/aimbat/ttconfig.py @@ -80,6 +80,8 @@ def __init__(self): self.srate = config.getfloat('sacplot', 'srate') self.tapertype = config.get('signal', 'tapertype') self.taperwidth = config.getfloat('signal', 'taperwidth') + self.thresholds = [ float(val) for val in config.get('sacplot', 'thresholds').split() ] + self.colorthresholds = config.get('sacplot', 'colorthresholds').split() @@ -109,6 +111,9 @@ def __init__(self): self.pickstyles = config.get('sacplot', 'pickstyles').split() self.minspan = config.getint('sacplot', 'minspan') self.fstack = config.get('iccs', 'fstack') + self.thresholds = [ float(val) for val in config.get('sacplot', 'thresholds').split() ] + self.colorthresholds = config.get('sacplot', 'colorthresholds').split() + class CCConfig: diff --git a/src/pysmo/aimbat/ttdefaults.conf b/src/pysmo/aimbat/ttdefaults.conf index bfe36eeb..91695f6c 100644 --- a/src/pysmo/aimbat/ttdefaults.conf +++ b/src/pysmo/aimbat/ttdefaults.conf @@ -47,7 +47,7 @@ qweights = 0.3333 0.3333 0.3333 srate = -1 ### cross-correlation xcorr_modu = xcorrf90 -xcorr_func = xcorr_fast +xcorr_func = xcorr_full_polarity shift = 10 maxiter = 10 ### convergence criterion @@ -67,7 +67,7 @@ srate = -1 ### output file: "$evdate".mc if mc ofilename = mc xcorr_modu = xcorrf90 -xcorr_func = xcorr_faster +xcorr_func = xcorr_full_polarity shift = 10 extraweight = 1000. ### lsqr solution without weighting; From 4faf3a714b715732922d15710e6f4c7b3813ff4b Mon Sep 17 00:00:00 2001 From: xsong Date: Fri, 29 Jul 2016 01:36:14 -0400 Subject: [PATCH 3/5] Fix missing configures in ttdefaults.conf --- src/pysmo/aimbat/ttdefaults.conf | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/pysmo/aimbat/ttdefaults.conf b/src/pysmo/aimbat/ttdefaults.conf index 91695f6c..3e476835 100644 --- a/src/pysmo/aimbat/ttdefaults.conf +++ b/src/pysmo/aimbat/ttdefaults.conf @@ -12,8 +12,8 @@ alphatwfill = 0.2 alphatwsele = 0.6 # plot picks t0-t5 npick = 6 -pickcolors = kmrcgyb -pickstyles = -- : +pickcolors = rkmcgyb +pickstyles = - : ### Fig size and rect for seismogram [x0, y0, dx, dy] figsize = 8 10 @@ -24,6 +24,10 @@ minspan = 5 ### sample rate. read first file if srate<0 srate = -1 +### thresholds for relative traveltime after mccc +thresholds = 3.5 4.0 5.0 +colorthresholds = green magenta red + ############################################## [sachdrs] ### time window @@ -49,12 +53,12 @@ srate = -1 xcorr_modu = xcorrf90 xcorr_func = xcorr_full_polarity shift = 10 -maxiter = 10 +maxiter = 20 ### convergence criterion ### type: coef for correlation coefficient, resi for residual ### epsi: small number convtype = coef -convepsi = 0.001 +convepsi = 0.0001 ### stack weight: coef for correlation coefficient, otherwise even weight stackwgt = coef From 5813ca29d26671431d72026f52eb6455e890369e Mon Sep 17 00:00:00 2001 From: xsong Date: Tue, 18 Apr 2017 21:42:25 -0400 Subject: [PATCH 4/5] Fix the following two bugs about MCCC results shown in GUI: - The MCCC shown on right end, the end of each seismogram, should be direct MCCC results, without substracting the t2 ( the picking arrival traveltime ) - The colors for different traveltime thresholds should be based on direct MCCC results, without substracting the t2 ( the picking arrival traveltime ) --- src/pysmo/aimbat/algmccc.py | 1 + src/pysmo/aimbat/pickphase.py | 6 +++--- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/pysmo/aimbat/algmccc.py b/src/pysmo/aimbat/algmccc.py index c569985b..d5de3718 100755 --- a/src/pysmo/aimbat/algmccc.py +++ b/src/pysmo/aimbat/algmccc.py @@ -372,6 +372,7 @@ def corrite(solist, mcpara, reftimes, solution, outvar, outcc, invdata): itmean = mean(reftimes) for i in range(nsta): wt = itmean + solution[i,0] + solist[i].sethdr('user5', solution[i,0]) solist[i].sethdr(wpick, wt) t0_times = [sacdh.gethdr('t0') for sacdh in solist] delay_times = [(solution[i,0]-(t0_times[i]-itmean)) for i in xrange(nsta)] diff --git a/src/pysmo/aimbat/pickphase.py b/src/pysmo/aimbat/pickphase.py index 91aaed88..03cbcee4 100755 --- a/src/pysmo/aimbat/pickphase.py +++ b/src/pysmo/aimbat/pickphase.py @@ -204,7 +204,7 @@ def labelStation(self): cc = sacdh.gethdr(hdrcc) sn = sacdh.gethdr(hdrsn) co = sacdh.gethdr(hdrco) - mc = sacdh.gethdr('user3') + mc = sacdh.gethdr('user5') slab += 'qual={0:5.3f}/{1:4.2f}/{2:.1f}/{3:4.2f}'.format(mc, cc, sn, co) trans = transforms.blended_transform_factory(axpp.transAxes, axpp.transData) font = FontProperties() @@ -464,8 +464,8 @@ def plotWave(self): colors = [[None,] * npsel , [None,] * npdel] for j in range(2): for k in range(nsede[j]): - plists[j][k].users[3] = plists[j][k].thdrs[3] - self.itmean - dt = plists[j][k].users[3] + #plists[j][k].users[3] = plists[j][k].thdrs[3] - self.itmean + dt = plists[j][k].users[5] if plists[j][k].selected: if abs(dt) > 10000 or abs(dt) < threshd[0]: colors[j][k] = colsel From 6cc98f17675560a2ab1ed23a3cae9687e49940e8 Mon Sep 17 00:00:00 2001 From: xsong Date: Mon, 28 Oct 2019 17:47:03 -0400 Subject: [PATCH 5/5] Add the usage of new features in the README.md --- README.md | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/README.md b/README.md index 56550051..283d1a0c 100644 --- a/README.md +++ b/README.md @@ -51,5 +51,23 @@ Running Unit Tests python run_unit_tests.py +The modification by Xin Song +------------------ +* Colour alarm for traces with very large traveltimes + +When the measured relative traveltime after MCCC exceeds the traveltime thresholds, the corresponding traces will show different colours. This feature is to avoid some high risky traces that might have very large traveltimes. It can have 3 traveltime thresholds maximum. + +This feature is set in the `ttconfig.py` in either `aimbat/src/pysmo/aimbat/` or either local directory + +For example in `ttconfig.py`: +``` +thresholds = 2 2.5 3 +colorthresholds = green magenta red +``` +It means the colour of trace will become green, magenta, or red if the traveltime of this trace is over 2s, 2.5s, or 3s after MCCC. + +* Perform the MCCC without changing the polarity of traces by setting `xcorr_func` in the `ttconfig.py`: + +`xcorr_func = xcorr_full_polarity`