[cig-commits] r16010 - seismo/3D/ADJOINT_TOMO/flexwin/scripts
liuqy at geodynamics.org
liuqy at geodynamics.org
Thu Nov 19 14:01:00 PST 2009
Author: liuqy
Date: 2009-11-19 14:01:00 -0800 (Thu, 19 Nov 2009)
New Revision: 16010
Modified:
seismo/3D/ADJOINT_TOMO/flexwin/scripts/process_ds.py
seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_out.py
Log:
process_ds.py ---
Choose data of the channel with largest SNR for flexwin input file
write_flexwin_out.py --
Allow merging of manual pick log file
Modified: seismo/3D/ADJOINT_TOMO/flexwin/scripts/process_ds.py
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/process_ds.py 2009-11-19 20:19:15 UTC (rev 16009)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/process_ds.py 2009-11-19 22:01:00 UTC (rev 16010)
@@ -1,13 +1,11 @@
#!/usr/bin/python
-# find matching data and synthetics from data and syn directories
+# find matching data and synthetics from data and syn directories, process them
# and write input to flexwin package
-# similar to Carl's process_data_and_syn.pl and prepare_input
-# works mainly for global dataset
-# Qinya Liu, Oct 2009, UofT
+# similar to Carl's process_data_and_syn.pl
+# Qinya Liu, Oct 2009
# remaining question: do we need -h for process_syn?
-
# back up your raw data and synthetics before processing
# channels may become a command-line option in the future
#
@@ -19,6 +17,24 @@
import os,sys,glob,getopt
+##########
+def calc_snr(dfile):
+ asc_file=dfile+'.a'
+ error=os.system("sac2ascii.pl "+dfile+" >/dev/null")
+ if (error != 0):
+ print "Error converting sac file"; exit()
+ t1= str(float(os.popen("saclst t1 f "+dfile).readline().split()[1])-5)
+ mm_nl=os.popen("awk '$1< "+ t1 + " {print $0}' " +asc_file+ " | minmax -C").readline().split()
+ if (len(mm_nl) != 4):
+ print "Error computing noise level"; exit()
+ max_nl=max(abs(float(mm_nl[2])),abs(float(mm_nl[3])))
+ mm_sn=os.popen("minmax -C "+asc_file).readline().split()
+ max_sn=max(abs(float(mm_sn[2])),abs(float(mm_sn[3])))
+ data_snr=max_sn/max_nl
+ os.system("rm -f ascfile")
+ return(data_snr)
+###########
+
usage= 'Usage: process_ds.py\n -d ddir,dext -s sdir,sext (required)\n -p cmt,sta -t tmin,tmax -m -w -D (optional)\n\n -m merge sac files from the same (sta,net,cmp,hole)\n -w generate flexwin input\n -D allow simultaneous processing of derivative synthetics\n'
preprocess=False; bandpass=False; merge=False; windowing=False; der_syn=False
@@ -70,15 +86,20 @@
merge_dir=datadir+'/merge_dir/'
if (not os.path.isdir(merge_dir)):
os.mkdir(merge_dir)
+
+# clean tmp files
+os.system('rm -f data.tmp '+datadir+'/*.c '+syndir+'/*.c '+syndir+'/*.c??')
#channel name
chan='LH'; sps='1.0';
+resp_dir='/data2/Datalib/global_resp/'
#padding zeros
bb=-40
if (bandpass):
bb=-float(tmin)*3.0
######### obtain data/syn list##################
+os.system('rm -f data.tmp '+datadir+'/*.c '+syndir+'/*.c '+syndir+'/*.c??')
nsta=0
data_list=[]; syn_list=[]; all_syn_list=[]; hole_list=[]; comp_list=[]
@@ -90,8 +111,14 @@
[sta,net,cmp,hole]=tmp
else:
[sta,net,cmp]=tmp; hole=''
+ if (list(cmp)[2] == '1'):
+ scmp=cmp[0:2]+'E'
+ elif (list(cmp)[2] == '2'):
+ scmp=cmp[0:2]+'N'
+ else:
+ scmp=cmp
dat=glob.glob(datadir+'/*'+net+'.'+sta+'.'+hole+'.'+cmp+'*'+dataext)
- syn=glob.glob(syndir+'/'+sta+'.'+net+'.'+cmp+synext)
+ syn=glob.glob(syndir+'/'+sta+'.'+net+'.'+scmp+synext)
if len(syn) != 1:
print 'no matching syn. for ', dat[0]
else:
@@ -118,10 +145,11 @@
all_syn_list.append(syn_list[i]+' '+' '.join(glob.glob(syn_list[i]+'.???')))
else:
all_syn_list=syn_list
-print "**** Total number of corresponding traces ", nsta, " *****"
+print "**** Total number of matching traces ", nsta, " *****"
############ preprocessing data and synthetics ################
if preprocess:
+
print '**** pre-processing data and synthetics ****'
error1=os.system('process_data_new.pl -m '+cmt+' -p -s '+sps+' '+' '.join(data_list)+'>/dev/null')
error2=os.system('process_syn_new.pl -S -m '+cmt+' -a '+st+' -p -s '+sps+' '+' '.join(all_syn_list)+'>/dev/null')
@@ -151,7 +179,7 @@
print 'check dt for',data, ' and ', syn; exit()
else:
# b=min(float(db),float(sb),bb); e=max(float(de),float(se))
- b=max(float(db),float(sb),bb); e=min(float(de),float(se))
+ b=max(float(db),float(sb),bb); e=max(min(float(de),float(se)),b+30)
sac_input=sac_input + 'cut '+str(b)+' '+str(e)+'\n' \
+'r '+data+' '+allsyn+'\n cut off\n w ' \
+ data_list[i]+'.c '+hole_ext.join(allsyn.split(' ')+[''])+'\n'
@@ -161,33 +189,41 @@
rdata_list=rdata_list+new_data_list[i]+' '
rsyn_list=rsyn_list+new_syn_list[i]+' '
if der_syn:
- rsyn_list=rsyn_list+new_syn_list[i]+'.??? '
+ rsyn_list=rsyn_list+' '+new_syndir+'/'+os.path.basename(syn)+'.???.c'+hole_list[i]+' '
sac_input=sac_input+'quit\n'
-# print 'sac <<EOF\n'+sac_input+'EOF\n'; exit()
if (os.system('sac <<EOF > /dev/null\n'+sac_input+'EOF\n') != 0):
print 'error cutting data and synthetics'; exit()
# further band-pass filtering data and all syn
print '**** band-passing data and synthetics ****'
- error1=os.system('process_data_new.pl -i -d '+new_datadir+' -t '+tmin+'/'+tmax+' '+cdata_list+'>/dev/null')
+# error1=os.system('process_data_new.pl -R '+resp_dir+' -d '+new_datadir+' -t '+tmin+'/'+tmax+' '+cdata_list+'>/dev/null ')
+ error1=os.system('process_data_new.pl -i -d '+new_datadir+' -t '+tmin+'/'+tmax+' '+cdata_list+'>/dev/null ')
error2=os.system('process_syn_new.pl -S -d '+new_syndir+' -t'+tmin+'/'+tmax+' '+csyn_list+'>/dev/null')
if (error1 != 0 or error2 != 0):
print 'Error bandpass filtering data and syn', error1, error2; exit()
# rotating LHE and LH1 components
print '**** rotating processed data and synthetics ****'
- error1=os.system('rotate.pl -d '+rdata_list+' >/dev/null')
+ error1=os.system('rotate.pl -d '+rdata_list+' >/dev/null ')
error2=os.system('rotate.pl '+rsyn_list+' >/dev/null')
if (error1 != 0 or error2 != 0):
print 'Error rotating raw data and syn', error1, error2; exit()
+ # correct derivative extension
+ if (der_syn):
+ for ext in ['Mrr','Mtt','Mpp','Mrt','Mrp','Mtp','dep','lat','lon']:
+ for file in glob.glob(new_syndir+'/*'+ext+'.c')+glob.glob(new_syndir+'/*'+ext+'.c??'):
+ file_list=os.path.basename(file).split('.')
+ os.rename(file,new_syndir+'/'+'.'.join(file_list[:-2]+[file_list[-1],file_list[-2]]))
# ############### flexwin input file ##########################
if windowing:
if (not bandpass):
print "use -t tmin,tmax before -w (flexwin input)"; exit()
print '**** writing flexwin input file ****'
- string=''
+ string=[]
+ j=-1
+ flexwin_dict={}
for i in range(0,nsta):
# flexwin input file
[sta,net,cmp]=os.popen('saclst kstnm knetwk kcmpnm f '+new_data_list[i]).readline().split()[1:]
@@ -198,16 +234,38 @@
new_data=glob.glob(new_datadir+'/*'+net+'.'+sta+'.'+hole_list[i]+'.'+cmp+'*'+dataext+'.c')
new_syn=glob.glob(new_syndir+'/'+sta+'.'+net+'.'+cmp+synext+'.c'+hole_list[i])
if (len(new_data) == 1 and len(new_syn) == 1):
- string=string+new_data[0]+'\n'+new_syn[0]+'\n'+'MEASURE/'+sta+'.'+net+'.'+cmp+'\n'
+ string.append(new_data[0]+'\n'+new_syn[0]+'\n'+'MEASURE/'+sta+'.'+net+'.'+cmp+'\n')
+ j=j+1
+ key=sta+'.'+net+'.'+cmp
+ if (key in flexwin_dict.keys()):
+ print 'multiple files for key ' + key + '===='
+ # compare data quality to decide data from which hole
+ data1=string[flexwin_dict[key]].split('\n')[0]; data2=string[j].split('\n')[0]
+ snr1=calc_snr(data1); snr2=calc_snr(data2)
+ if (snr1 < snr2):
+ flexwin_dict[key]=j
+ print 'Keeping '+data2+' instead of '+data1+', determined by SNR', snr2, snr1
+ else:
+ print 'Keeping '+data1+' over '+data2+', determined by SNR', snr1, snr2
+ else:
+ flexwin_dict[key]=j
else:
print "No exact matching for ", sta, net, cmp, hole_list[i]
+
+ nkeys=len(flexwin_dict.keys())
+ outstring=str(nkeys)+'\n'
+ for key in sorted(flexwin_dict,key=flexwin_dict.__getitem__):
+ outstring=outstring+string[flexwin_dict[key]]
+
f=open('INPUT_flexwin','w')
- f.write(str(nsta)+'\n')
- f.write(string.rstrip())
+ f.write(outstring.rstrip('\n'))
f.close()
+
if (not os.path.isdir('MEASURE')):
os.mkdir('MEASURE')
os.system('rm -f data.tmp '+datadir+'/*.c '+syndir+'/*.c '+syndir+'/*.c??')
print '**** Done ****'
+
+
Modified: seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_out.py
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_out.py 2009-11-19 20:19:15 UTC (rev 16009)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_out.py 2009-11-19 22:01:00 UTC (rev 16010)
@@ -1,23 +1,103 @@
#!/usr/bin/python
# similar to Carl's prepare_meas_all.pl script
# Qinya Liu, Oct 2009, UofT
+
+# adding capability to apply a manual pick file, Nov 2009
+# manual pick file format
+# input_flexwin --- line 1
+# - KIEV.IU.LHZ 1 --- delete window 1 pick for KIEV.IU.LHZ
+# + UCH.KN.LHR 400 725 --- add window pick between [400,725] for UCH.KN.LHR
+# possibly can be used for automation after manual picking in a GUI interface
+
import os,sys,glob
-if (len(sys.argv) != 3):
- print "write_flexwin_out.py measure_dir out_filename"; exit()
+if (len(sys.argv) != 3 and len(sys.argv) != 4):
+ print "write_flexwin_out.py measure_dir out_filename [manual-pick-file]"; exit()
-dir=sys.argv[1]; out=sys.argv[2]
+dir=sys.argv[1]; outfile=sys.argv[2]
if (not os.path.isdir(dir)):
print 'check if '+dir+' exists or not'; exit()
+manual_pick=False
+if (len(sys.argv) == 4):
+ manual_pick=True
+ manual_file=sys.argv[3]
+ if (not os.path.isfile(manual_file)):
+ print 'Check if '+manual_file+' exists or not'; exit()
+ manual_lines=open(manual_file,'r').readlines()
+ flexwin_input = manual_lines[0].rstrip('\n');
+ print '*** Adding manual pick info from '+ manual_file+ ' ***'
+
output=''
files=glob.glob(dir+'/*mt*');
-nfiles=len(files);
+nfiles=0
+
for file in files:
- output=output+''.join(os.popen('cat '+file).readlines())
+ if (not manual_pick):
+ nfiles=nfiles+1
+ output=output+''.join(os.popen('cat '+file).readlines())
+
+ else: # manual pick file
+ basefile='.'.join(os.path.basename(file).split('.')[0:3])
+# print 'processing '+basefile+' ...'
+ out=os.popen('grep '+basefile+' '+manual_file).readlines()
+ if (len(out) == 0):
+ output=output+''.join(os.popen('cat '+file).readlines())
+ nfiles=nfiles+1
+ elif (len(out) == 1):
+ out2=out[0].rstrip('\n').split(' ') # manual pick info
+ mtfile=os.popen('cat '+file).readlines()
+ nwin=int(mtfile[2])
+ if (out2[0] == '-'): # delete window(s)
+ nwin_subtract=len(out2)-2
+ print 'Deleting windows from '+file, nwin, nwin_subtract
+ if (nwin > nwin_subtract): # delete some windows
+ nfiles=nfiles+1
+ output=output+mtfile[0]+mtfile[1]+' '+str(nwin-nwin_subtract)+'\n'
+ for i in range(0,nwin):
+ if not (str(i+1) in out2[2:]):
+ output=output+mtfile[i+3]
+ elif (out2[0] == '+'): # add windows
+ nwin_add=(len(out2)-2)/2
+ nfiles=nfiles+1
+ print 'Adding windows for '+file, nwin, nwin_add
+ output=output+mtfile[0]+mtfile[1]+' '+str(nwin+nwin_add)+'\n'
+ output=output+''.join(mtfile[3:])
+ for i in range(0,nwin_add):
+ if (float(out2[i*2+2]) < float(out2[i*2+3])):
+ output=output+' '+out2[i*2+2]+' '+out2[i*2+3]+'\n'
+ else:
+ print 'End time < start time in line: '+ out[0]; exit()
+ else:
+ print 'Format error in line: '+out[0]; exit()
+ else:
+ print 'Multiple lines for '+basefile; exit()
+# now add files that did not have mt_input
+if manual_pick:
+ print ' *** Adding missing data and syn pairs *** '
+ f=open(outfile, 'w')
+ f.write(output)
+ f.close()
+ for i in range(1,len(manual_lines)):
+ out2=manual_lines[i].split(' ')
+ if (os.system('cat '+outfile+' | grep '+out2[1]+'>/dev/null') != 0 and out2[0] == '+'): # not an existing window
+ out3=os.popen('grep -n "'+dir+'.*'+out2[1]+'" '+flexwin_input).readlines()
+ if (len(out3) != 1):
+ print 'Error grep '+out2[1]+' from '+flexwin_input+' : '+synfile; exit()
+ else:
+ line_num=int(out3[0].split(':')[0])
+ datafile=os.popen("awk 'NR=="+str(line_num-2)+" {print $0}' "+flexwin_input).readline()
+ synfile=os.popen("awk 'NR=="+str(line_num-1)+" {print $0}' "+flexwin_input).readline()
+ print 'Adding windows to '+datafile+' and '+synfile
+ nfiles=nfiles+1
+ nwin_add=(len(out2)-2)/2
+ output=output+datafile+synfile+str(nwin_add)+'\n'
+ for i in range(0,nwin_add):
+ output=output+' '+out2[i*2+2]+' '+out2[i*2+3]+'\n'
+
output=str(nfiles)+'\n'+output
-f=open(out, 'w')
+f=open(outfile, 'w')
f.write(output)
f.close()
More information about the CIG-COMMITS
mailing list