[cig-commits] r19435 - seismo/3D/ADJOINT_TOMO/flexwin/scripts
liuqy at geodynamics.org
liuqy at geodynamics.org
Tue Jan 24 09:30:36 PST 2012
Author: liuqy
Date: 2012-01-24 09:30:35 -0800 (Tue, 24 Jan 2012)
New Revision: 19435
Added:
seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_manual_pick.py
Modified:
seismo/3D/ADJOINT_TOMO/flexwin/scripts/extract_event_windowing_stats.sh
seismo/3D/ADJOINT_TOMO/flexwin/scripts/process_ds.py
seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_out.py
seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_shift_weight.py
Log:
Update scripts to accomodate more general cases.
Modified: seismo/3D/ADJOINT_TOMO/flexwin/scripts/extract_event_windowing_stats.sh
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/extract_event_windowing_stats.sh 2012-01-24 17:28:40 UTC (rev 19434)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/extract_event_windowing_stats.sh 2012-01-24 17:30:35 UTC (rev 19435)
@@ -110,6 +110,7 @@
# set output filename
out=${basename}/event_winstats.eps
+outpdf=${basename}/event_winstats.pdf
short_basename=`echo $basename | awk -F"/" '{print $NF}'`
#evlo=`echo $evlo | awk '{printf "%.0f", $1}'`
@@ -143,10 +144,11 @@
mv t1 $out
echo $out
-
+ps2pdf $out $outpdf
# plot window recordsection
out2=${basename}/event_recordsection.eps
+out2pdf=${basename}/event_recordsection.pdf
gmtset PAPER_MEDIA a4+ MEASURE_UNIT inch
region=$min_time/$max_time/$min_ddg/$max_ddg
@@ -168,6 +170,7 @@
mv t1 $out2
echo $out2
+ps2pdf $out2 $out2pdf
rm -f $t0 $t1 $t2z $t2r $t2t $t3
Modified: seismo/3D/ADJOINT_TOMO/flexwin/scripts/process_ds.py
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/process_ds.py 2012-01-24 17:28:40 UTC (rev 19434)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/process_ds.py 2012-01-24 17:30:35 UTC (rev 19435)
@@ -9,9 +9,11 @@
# back up your raw data and synthetics before processing
# channels may become a command-line option in the future
# when number of matching files is large, we may face issues of
-# 'argument list too long' (especially when derivatives are includede)
+# 'argument list too long' (especially when derivatives are included)
+# a. data/syn_tmin_tmax are wiped out before processing
#
-# need process_data_new.pl, process_syn_new.pl, pad_zeros.pl and rotate.pl as well as sac_merge
+# need process_data_new.pl, process_syn_new.pl, pad_zeros.pl and rotate.pl
+# as well as saclst sacch (?), and sac_merge
# data_list,syn_list -- input data/syn
# all_syn_list -- input syn+derivative syn (.???)
@@ -24,11 +26,11 @@
asc_file=dfile+'.a'
error=os.system("sac2ascii.pl "+dfile+" >/dev/null")
if (error != 0):
- print "Error converting sac file"; exit()
+ sys.exit("Error converting sac file")
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()
+ sys.exit("Error computing noise level")
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])))
@@ -37,80 +39,105 @@
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'
+usage= 'Usage: process_ds.py\n -d ddir,dext,chan -s sdir,sext,chan (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 -r (for regional data) -M measure_dir\n'
preprocess=False; bandpass=False; merge=False; windowing=False; der_syn=False
+regional_data=False; model_name='global'; eps=0.001; pzdir=''; rot_reg='-d'; chan='LH'; synchan=chan; meas_dir='MEASURE'
# command-line options
try:
- opts,arg=getopt.getopt(sys.argv[1:],'d:s:p:t:mwD')
+ opts,arg=getopt.getopt(sys.argv[1:],'d:s:p:t:mwDrM:')
except getopt.GetoptError:
- print usage; exit()
+ sys.exit(usage)
if (len(opts) == 0):
- print usage; exit()
+ sys.exit(usage)
for o,a in opts:
if o == '-d':
- (datadir,dataext)=a.split(',')
+ (datadir,dataext,chan)=a.split(',')
if o == '-s':
- (syndir,synext)=a.split(',')
+ aa=a.split(',')
+ syndir=aa[0]
+ if (len(aa) >= 2):
+ synext=aa[1]
+ if (len(aa) == 3):
+ synchan=aa[2]
if (synext != ''):
synext='.'+synext
+ if (len(aa) == 2):
+ synchan=chan
if o == '-p':
preprocess=True
(cmt,st)=a.split(',')
if (cmt == '' or st== '' or not os.path.isfile(cmt) or not os.path.isfile(st)):
- print 'check cmt and sta file'; exit()
+ sys.exit('check cmt and sta file')
if o == '-t':
bandpass=True
(tmin,tmax)=a.split(',')
- junk=float(tmin); junk=float(tmax)
+ junk=float(tmin); junk=float(tmax) #?????
ext='_'+tmin+'_'+tmax
new_datadir=datadir+ext; new_syndir=syndir+ext
if (not os.path.isdir(new_datadir)):
os.mkdir(new_datadir)
- else:
- os.system('rm -rf '+new_datadir+'/*')
if (not os.path.isdir(new_syndir)):
os.mkdir(new_syndir)
- else:
- os.system('rm -rf '+new_syndir+'/*')
+ os.system('rm -rf '+new_datadir+'/*')
+ os.system('rm -rf '+new_syndir+'/*')
if o == '-m':
merge=True
if o == '-w':
windowing=True
if o == '-D':
der_syn=True
+ if o == '-r':
+ regional_data=True
+ model_name='socal'
+ pzdir='/opt/seismo/datalib/CALIFORNIA_PZS'; rot_reg=''
+ if o == '-M':
+ meas_dir=a
if (not os.path.isdir(datadir) or not os.path.isdir(syndir)):
- print datadir, ' ', syndir, 'exist or not'; exit()
+ sys.exit(datadir+' '+syndir+' exist or not')
if merge:
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??')
-
+if (windowing and not bandpass):
+ sys.exit('pass bands -t has to be set before windowing (-w)')
+if (windowing):
+ meas_dir=meas_dir+ext
+ input_flexwin='INPUT_FLEXWIN'+ext
+#if (preprocess and bandpass): # only clean up directories if preprocessing is done
+# os.system('rm -rf '+new_datadir+'/*')
+# os.system('rm -rf '+new_syndir+'/*')
#channel name
-chan='LH'; sps='1.0';
-resp_dir='/data2/Datalib/global_resp/'
+if chan == 'LH':
+ sps='20'
+else:
+ sps='1'
+ resp_dir='/data2/Datalib/global_resp/'
+
#padding zeros
-bb=-40
+if (regional_data):
+ bb=-40
+else:
+ bb=-100
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
+bb2=bb+1.0 # to make sure cut range is inside [b,e]
+
+######### obtain data/syn list (sort by station) ##################
+nfile=0; outfile='out_ds.tmp'; os.system('rm -f '+outfile)
data_list=[]; syn_list=[]; all_syn_list=[]; hole_list=[]; comp_list=[]
-print '**** Data/syn list for the same (sta,net,LH?,khole) ****'
+print '**** Data/syn list for the same (sta,net,'+chan+'?,khole) ****'
for dd in os.popen('saclst kstnm knetwk kcmpnm khole f '+datadir+'/*'+chan+'[ZEN12]*'+dataext+"| awk '{print $2,$3,$4,$5}' | sort -k 1 | uniq").readlines():
tmp=dd.split()
if (len(tmp) > 3):
[sta,net,cmp,hole]=tmp
+ if abs(float(hole)+12345)<0.1:
+ hole=''
else:
[sta,net,cmp]=tmp; hole=''
if (list(cmp)[2] == '1'):
@@ -119,126 +146,132 @@
scmp=cmp[0:2]+'N'
else:
scmp=cmp
- dat=glob.glob(datadir+'/*'+net+'.'+sta+'.'+hole+'.'+cmp+'*'+dataext)
- syn=glob.glob(syndir+'/'+sta+'.'+net+'.'+scmp+synext)
+ if regional_data: # the naming convention is subject to change
+ dat=glob.glob(datadir+'/*'+net+'.'+sta+'.'+cmp+'*'+dataext)
+ else:
+ dat=glob.glob(datadir+'/*'+net+'.'+sta+'.'+hole+'.'+cmp+'*'+dataext)
+ syn=glob.glob(syndir+'/'+sta+'.'+net+'.'+synchan+list(scmp)[2]+synext)
if len(syn) != 1:
- print 'no matching syn. for ', dat[0]
+ print 'no matching syn. for ', dat[0], ':', syndir+'/'+sta+'.'+net+'.'+synchan+list(scmp)[2]+synext
else:
if len(dat) > 1:
print 'multiple files for merging ===', dd.rstrip()
if merge:
dat.sort()
- os.system('process_data_new.pl -m '+cmt+' '+' '.join(dat)+'>/dev/null')
- error=os.system('sac_merge '+' '.join(dat)+' merge_tmp.sac')
- if error != 0:
+ os.system('process_data_new.pl -m '+cmt+' '+' '.join(dat)+'>>'+outfile)
+ if os.system('sac_merge '+' '.join(dat)+' merge_tmp.sac') != 0:
print 'cannot merge sac files ', dat
os.system('mv -f '+' '.join(dat[1:])+' '+merge_dir)
else:
print 'successfully merging sac files ', dat
os.system('mv -f '+' '.join(dat)+' '+merge_dir+';mv -f merge_tmp.sac '+dat[0]);
else:
- print 'sac files need to be merged or resolved', ' '.join(dat), ' use -m'; exit()
- data_list.append(dat[0]); syn_list.append(syn[0])
+ sys.exit('sac files need to be merged or resolved', ' '.join(dat), ' use -m');
+ data_list.append(dat[0]); syn_list.append(syn[0]) # ENZ components
hole_list.append(hole); comp_list.append(cmp)
- nsta=nsta+1
+ nfile=nfile+1
if der_syn:
for i in range(0,len(syn_list)):
all_syn_list.append(syn_list[i]+' '+' '.join(glob.glob(syn_list[i]+'.???')))
else:
all_syn_list=syn_list
-print "**** Total number of matching traces ", nsta, " *****"
+print "**** Total number of matching traces ", nfile, " *****"
############ 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')
- if (error1 != 0 or error2 != 0):
- print 'Error pre-processing data and syn'; exit()
- error1=os.system('pad_zeros.pl -s -l '+str(bb)+' '+' '.join(all_syn_list)+'> /dev/null')
- if (error1 != 0):
- print 'Error padding zeros to syn'; exit()
-
+ data_cmd='process_data_new.pl -m '+cmt+' -p --model='+model_name
+ print data_cmd+' data_files'
+ error=os.system(data_cmd+' '+' '.join(data_list)+'>>'+outfile)
+ if error != 0:
+ sys.exit('Error pre-processing data')
+
+ # why can't we run this separately, because of the limitation on arg length?
+ f=open('syn.command','w')
+ syn_cmd='process_syn_new.pl -S -m '+cmt+' -a '+st+' -p --model='+model_name
+ print syn_cmd+' all_syn_files'
+ f.write(syn_cmd+' '+' '.join(all_syn_list)+' >> '+outfile+'\n')
+ pad_cmd='pad_zeros.pl -s -l '+str(bb)
+ print pad_cmd+' all_syn_files'
+ f.write(pad_cmd+' '+' '.join(all_syn_list)+' >>'+outfile+'\n')
+ f.close()
+ if os.system('chmod a+x syn.command; syn.command') != 0:
+ sys.exit('Error pre-processing syn and padding synthetics')
+
########### cutting and band-passing data and synthetics #############
-new_data_list=[]; new_syn_list=[]
-cdata_list=''; csyn_list=''; rdata_list=''; rsyn_list=''
-
+cdata_list=[];
if bandpass:
print '**** cutting data and synthetics ****'
- sac_input='cuterr fillz\n'
- for i in range(0,nsta):
- data=data_list[i]; syn=syn_list[i]; allsyn=all_syn_list[i].rstrip()
- new_data_list.append(new_datadir+'/'+os.path.basename(data)+'.c')
- new_syn_list.append(new_syndir+'/'+os.path.basename(syn)+'.c'+hole_list[i])
- hole_ext='.c'+hole_list[i]+' '
+ sac_input='cuterr fillz\n' # this should be redundant
+ for i in range(0,nfile):
+ data=data_list[i]; syn=syn_list[i]; allsyn=all_syn_list[i].rstrip()
[db,de,ddt]=os.popen('saclst b e delta f '+data).readline().split()[1:]
[sb,se,sdt]=os.popen('saclst b e delta f '+syn).readline().split()[1:]
- if (float(ddt)-float(sdt)) > 0.001:
- print 'check dt for',data, ' and ', syn; exit()
+
+ if (float(ddt)-float(sdt)) > eps:
+ sys.exit('check dt for '+data+' and '+syn)
else:
-# b=min(float(db),float(sb),bb); e=max(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'
- cdata_list=cdata_list+data_list[i]+'.c '
- csyn_list=csyn_list+' '+hole_ext.join(allsyn.split(' ')+[''])
- if (list(comp_list[i])[2] == 'E' or list(comp_list[i])[2] == '1'):
- rdata_list=rdata_list+new_data_list[i]+' '
- rsyn_list=rsyn_list+new_syn_list[i]+' '
- if der_syn:
- rsyn_list=rsyn_list+' '+new_syndir+'/'+os.path.basename(syn)+'.???.c'+hole_list[i]+' '
- sac_input=sac_input+'quit\n'
- f = open('sac.cmd', "w"); f.write('#!/bin/bash\nsac<<EOF\n')
- f.write(''.join(sac_input)); f.write('EOF\n'); f.close()
- if (os.system('chmod a+rx sac.cmd; sac.cmd') != 0):
- print 'error cutting data and synthetics '+str(nerror); exit()
+ b=max(float(db),float(sb),bb2); e=max(min(float(de),float(se)),b+30)
+ if (os.system('sacch t3 '+str(b)+' t4 '+str(e)+' f '+data+' '+allsyn+'>>'+outfile)!=0):
+ sys.exit('error setting start and end time to t3 and t4')
+ cdata_list.append(new_datadir+'/'+os.path.basename(data));
# further band-pass filtering data and all syn
print '**** band-passing data and synthetics ****'
-# 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()
+ data_cmd= 'process_data_new.pl -l t3/t4 -d '+new_datadir+' -t '+tmin+'/'+tmax+' -i '+pzdir+' -f -s '+sps+' --model='+model_name
+ print data_cmd+' data_files'
+ error=os.system(data_cmd+' '+' '.join(data_list)+'>>'+outfile)
+ if (error != 0):
+ sys.exit('Error bandpass filtering data '+str(error1))
+
+ f=open('syn.command','w')
+ syn_cmd='process_syn_new.pl -S -l t3/t4 -d '+new_syndir+' -t '+tmin+'/'+tmax+' -f -s '+sps+' --model='+model_name
+ print syn_cmd+' all_syn_files'
+ f.write(syn_cmd+' '+' '.join(all_syn_list)+'>>'+ outfile+'\n')
+ f.close()
+ if os.system('chmod a+x syn.command; syn.command') != 0:
+ sys.exit('Error filtering syn ')
# rotating LHE and LH1 components
print '**** rotating processed data and synthetics ****'
- error1=os.system('rotate.pl -d '+rdata_list+' >/dev/null ')
- error2=os.system('rotate.pl '+rsyn_list+' >/dev/null')
+ rot_cmd='rotate.pl '+rot_reg+' '+new_datadir+'/*'+chan+'[E1]*'+dataext
+ print rot_cmd
+ error1=os.system(rot_cmd+' >>'+outfile)
+
+ rot_cmd='rotate.pl '+new_syndir+'/*'+synchan+'[E1]'+synext+'*'# '+new_syndir+'/*'+synchan+'[E1]'+synext+'.???'
+ print rot_cmd
+ error2=os.system(rot_cmd+' >>'+outfile)
+
if (error1 != 0 or error2 != 0):
- print 'Error rotating raw data and syn', error1, error2; exit()
+ sys.exit('Error rotating raw data and syn '+str(error1)+' '+str(error2))
- # 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()
+ sys.exit("make sure -t tmin,tmax is used before -w (flexwin input)")
print '**** writing flexwin input file ****'
string=[]
j=-1
flexwin_dict={}
- for i in range(0,nsta):
+ for i in range(0,nfile):
# flexwin input file
- [sta,net,cmp]=os.popen('saclst kstnm knetwk kcmpnm f '+new_data_list[i]).readline().split()[1:]
+ [sta,net,cmp]=os.popen('saclst kstnm knetwk kcmpnm f '+cdata_list[i]).readline().split()[1:]
if (list(cmp)[2] == 'E' or list(cmp)[2] == '1'):
cmp=''.join(list(cmp[0:2])+['T'])
elif (list(cmp)[2] == 'N' or list(cmp)[2] == '2'):
cmp=''.join(list(cmp[0:2])+['R'])
- 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 (not regional_data):
+ new_data=glob.glob(new_datadir+'/*'+net+'.'+sta+'.'+hole_list[i]+'.'+cmp+'*'+dataext)
+ else:
+ new_data=glob.glob(new_datadir+'/*'+net+'.'+sta+'.'+cmp+'*'+dataext)
+ new_syn=glob.glob(new_syndir+'/'+sta+'.'+net+'.'+cmp+synext)
if (len(new_data) == 1 and len(new_syn) == 1):
- string.append(new_data[0]+'\n'+new_syn[0]+'\n'+'MEASURE/'+sta+'.'+net+'.'+cmp+'\n')
+ string.append(new_data[0]+'\n'+new_syn[0]+'\n'+meas_dir+'/'+sta+'.'+net+'.'+cmp+'\n')
j=j+1
key=sta+'.'+net+'.'+cmp
if (key in flexwin_dict.keys()):
@@ -255,21 +288,21 @@
flexwin_dict[key]=j
else:
print "No exact matching for ", sta, net, cmp, hole_list[i]
+# check the beginning and end times
-
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(outstring.rstrip('\n'))
+ f=open(input_flexwin,'w')
+ f.write(outstring)
f.close()
- if (not os.path.isdir('MEASURE')):
- os.mkdir('MEASURE')
+ if (not os.path.isdir(meas_dir)):
+ os.mkdir(meas_dir)
-os.system('rm -f data.tmp sac.cmd '+datadir+'/*.c '+syndir+'/*.c '+syndir+'/*.c??')
+#os.system('rm -f sac.cmd syn.command')
print '**** Done ****'
Added: seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_manual_pick.py
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_manual_pick.py (rev 0)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_manual_pick.py 2012-01-24 17:30:35 UTC (rev 19435)
@@ -0,0 +1,105 @@
+#!/usr/bin/env 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 and len(sys.argv) != 4):
+ sys.exit("write_flexwin_out.py measure_dirs out_filename [manual-pick-file]")
+
+dir=sys.argv[1]; outfile=sys.argv[2]
+if (not os.path.isdir(dir)):
+ sys.exit('check if '+dir+' exists or not')
+
+manual_pick=False
+if (len(sys.argv) == 4):
+ manual_pick=True
+ manual_file=sys.argv[3]
+ if (not os.path.isfile(manual_file)):
+ sys.exit('Check if '+manual_file+' exists or not')
+ 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*')
+if len(files) <= 0:
+ sys.exit('No such file as '+dir+'/*mt*')
+nfiles=0
+
+for file in files:
+ 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(outfile, 'w')
+f.write(output)
+f.close()
+
Property changes on: seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_manual_pick.py
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_out.py
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_out.py 2012-01-24 17:28:40 UTC (rev 19434)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_out.py 2012-01-24 17:30:35 UTC (rev 19435)
@@ -2,104 +2,30 @@
# 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 and len(sys.argv) != 4):
- sys.exit("write_flexwin_out.py measure_dir out_filename [manual-pick-file]")
+if (len(sys.argv) < 3):
+ sys.exit("write_flexwin_out.py measure_dirs out_filename ")
-dir=sys.argv[1]; outfile=sys.argv[2]
-if (not os.path.isdir(dir)):
- sys.exit('check if '+dir+' exists or not')
+dirs=sys.argv[1:-1]; outfile=sys.argv[-1]
+f=open(outfile, 'w')
-manual_pick=False
-if (len(sys.argv) == 4):
- manual_pick=True
- manual_file=sys.argv[3]
- if (not os.path.isfile(manual_file)):
- sys.exit('Check if '+manual_file+' exists or not')
- 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*')
-if len(files) <= 0:
- sys.exit('No such file as '+dir+'/*mt*')
nfiles=0
+output=''
+for dir in dirs:
+ if (not os.path.isdir(dir)):
+ sys.exit('check if '+dir+' exists or not')
+
+ files=glob.glob(dir+'/*mt*')
+ if len(files) <= 0:
+ sys.exit('No such file as '+dir+'/*mt*')
-for file in files:
- if (not manual_pick):
+ for file in files:
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(outfile, 'w')
f.write(output)
f.close()
Modified: seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_shift_weight.py
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_shift_weight.py 2012-01-24 17:28:40 UTC (rev 19434)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_shift_weight.py 2012-01-24 17:30:35 UTC (rev 19435)
@@ -2,18 +2,22 @@
# note: geometrical spreading factor needs to be significantly improved,
# right now only scale by 1/dist**exp
from numpy import *
-import os,sys,commands,glob,getopt
+import os,sys,commands,glob,getopt,re
# command-line options
usage="""
write_flexwin_shift_weight.py [-C]
[-c Z/R/T -d B/R/L -a az_weight -r corr_weight -t ]
- -s MEASURE/ [-o adj-dir]
+ -D data_dir -s measure_dir [-o adj-dir]
flexwin_summary_file new_flexwin_summary_file
+ -c Z/R/T component weights; -d B/R/L body/rayleigh/love wave weights
+ -a azimuthal weight -r correlation weight -t apply time shift
+ -C socal (regional) weight scheme in which -d pnl/rayleigh/love waves
+ -D and -s gives data and measure dir prefixes
"""
try:
- opts,args=getopt.getopt(sys.argv[1:],'c:d:a:s:r:o:tC')
+ opts,args=getopt.getopt(sys.argv[1:],'c:d:a:D:s:r:o:tC')
except getopt.GetoptError:
sys.exit(usage)
if len(args) != 2:
@@ -22,7 +26,7 @@
# default weight values
comp_z_weight=1.; comp_r_weight=1.; comp_t_weight=1.
pnl_dist_weigth=0.; rayleigh_dist_weight=0.; love_dist_weight=0.
-az_exp_weigth=0.; corr_exp_weight = 0.
+az_exp_weigth=0.; corr_exp_weight = 0.; data_dir='data'
adj_dir='.'; write_tshift=False; write_adj_dir=False; socal=False
#opts
@@ -30,12 +34,12 @@
if o == '-c':
tmp=a.split('/')
if len(tmp)!=3:
- sys.exit('-c Z/R/T')
+ sys.exit('-c Z/R/T '+a)
comp_z_weight=float(tmp[0]); comp_r_weight=float(tmp[1]); comp_t_weight=float(tmp[2])
if o == '-d':
tmp=a.split('/')
if len(tmp)!=3:
- sys.exit('-c P/R/L')
+ sys.exit('-d P/R/L '+a)
pnl_dist_weight =float(tmp[0]); rayleigh_dist_weight=float(tmp[1]); love_dist_weight=float(tmp[2])
if o == '-a':
az_exp_weight=float(a)
@@ -43,8 +47,8 @@
corr_exp_weight=float(a)
if o == '-s':
measure_dir=a
- if not os.path.isdir(measure_dir):
- sys.exit('Check if '+measure_dir+' exist or not')
+# if not os.path.isdir(measure_dir): # measure_dir is just a prefix
+# sys.exit('Check if '+measure_dir+' exist or not')
if o == '-o':
write_adj_dir=True
adj_dir=a.rstrip('/')
@@ -54,6 +58,8 @@
write_tshift=True
if o == '-C':
socal=True
+ if o == '-D':
+ data_dir=a # data_dir is a prefix
infile=args[0]; outfile=args[1]
@@ -71,12 +77,12 @@
ref_dist=60 # degree
vr=3.1; vl=4.1
naz=10; daz=360/naz;
-eps=1.0e-4; pi=3.1416926
+eps=1.0e-2; pi=3.1416926;
naz_array=ones((naz))
-alldata=commands.getoutput('grep data '+infile).split('\n')
+alldata=commands.getoutput('grep '+data_dir+' '+infile).split('\n')
if (len(alldata) == 0):
- sys.exit('Right now I assume data directories contain the name "data"')
+ sys.exit('No selected files of '+data_dir+' in '+infile)
output=commands.getoutput('saclst az f '+' '.join(alldata)+"|awk '{print $2}'").split('\n')
nfile_az=len(output)
for i in range(0,nfile_az):
@@ -105,13 +111,14 @@
# read data and synthetics, check headers
[status1,out1]=commands.getstatusoutput('saclst kstnm knetwk kcmpnm az dist gcarc f '+datfile)
- [status2,out2]=commands.getstatusoutput('saclst kstnm knetwk kcmpnm az dist f '+synfile)
+ [status2,out2]=commands.getstatusoutput('saclst kstnm knetwk kcmpnm az dist gcarc f '+synfile)
if status1 != 0 or status2 != 0:
sys.exit('Error taking kcmpnm az dist from '+datfile+' and '+synfile)
[sta,net,comp,az,dist,gcarc]=out1.split()[1:7]
az=float(az); dist=float(dist); gcarc=float(gcarc)
- [sta2,net2,comp2,az2,dist2]=out2.split()[1:6]; az2=float(az2); dist2=float(dist2)
- if (sta != sta2 or net != net2 or comp != comp2 or abs(az-az2) >eps or abs(dist-dist2) > eps):
+ [sta2,net2,comp2,az2,dist2,gcarc2]=out2.split()[1:7]
+ az2=float(az2); dist2=float(dist2); gcarc2=float(gcarc2)
+ if (sta != sta2 or net != net2 or comp != comp2 or abs(az-az2) >eps or abs(dist-dist2)/dist > eps):
# sys.exit('Check if kstnm knetwk kcmpnm az dist are the same for '+datfile+' and '+synfile)
print 'Unmatching sta, net, comp, az, dist for =====\n',out1, '\n', out2
@@ -125,7 +132,7 @@
else:
sys.exit('Only deal with Z/R/T components right now')
- filename=sta+'.'+net+'.'+comp+'.win.qual'
+ filename=re.sub(data_dir,measure_dir,os.path.dirname(datfile))+'/'+ sta+'.'+net+'.'+comp+'.win.qual'
if write_adj_dir:
g.write(adj_dir+'/'+sta+'.'+net+'.'+comp+'.adj.sac\n')
@@ -158,7 +165,7 @@
# azimuth
nn=3+j;
- [istat,out]=commands.getstatusoutput("awk 'NF == 6 {print $2,$3,$4,$5}' "+measure_dir+'/'+filename)
+ [istat,out]=commands.getstatusoutput("awk 'NF == 6 {print $2,$3,$4,$5}' "+filename)
if (istat != 0 or len(out)== 0):
sys.exit('Error with '+"awk 'NF == 6 {print $2,$3,$4,$5}' "+measure_dir+'/'+filename)
More information about the CIG-COMMITS
mailing list