[cig-commits] r15901 - in seismo/3D/ADJOINT_TOMO/flexwin: . scripts
liuqy at geodynamics.org
liuqy at geodynamics.org
Fri Oct 30 14:17:57 PDT 2009
Author: liuqy
Date: 2009-10-30 14:17:57 -0700 (Fri, 30 Oct 2009)
New Revision: 15901
Added:
seismo/3D/ADJOINT_TOMO/flexwin/scripts/plot_flexwin.py
seismo/3D/ADJOINT_TOMO/flexwin/scripts/process_ds.py
seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_out.py
Modified:
seismo/3D/ADJOINT_TOMO/flexwin/seismo_subs.f90
Log:
Fixing sta,net,cmp output problem in seismo_subs.f90(no @ anymore)
Add scripts for processing data/syn, generating input file, graphics output, and compact window file output (mainly for global dataset)
Added: seismo/3D/ADJOINT_TOMO/flexwin/scripts/plot_flexwin.py
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/plot_flexwin.py (rev 0)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/plot_flexwin.py 2009-10-30 21:17:57 UTC (rev 15901)
@@ -0,0 +1,33 @@
+#!/usr/bin/python
+# Similar to Carl's plot_windows_all.pl script
+# Qinya Liu, Oct 2009, UofT
+# plots are sort by distance and then components
+
+import sys,os
+
+if (len(sys.argv) != 2):
+ print "plot_flexwin.py measure_dir(MEASURE)"; exit()
+
+dir=sys.argv[1]
+if not os.path.isdir(dir):
+ print "no such dir "+dir; exit()
+script1='scripts/plot_seismos_gmt.sh'
+script2='scripts/extract_event_windowing_stats.sh'
+if not os.path.isfile(script1) or not os.path.isfile(script2):
+ print 'no '+script1+' or '+script2; exit()
+
+if (os.system(script2+' '+dir+' > /dev/null') != 0):
+ print 'Error executing '+ script2; exit()
+
+ps=dir+'/event_winstats.pdf'+' '+dir+'/event_recordsection.pdf'
+
+for basename in os.popen('grep DDG '+dir+"/*.info | awk '{print $1,$4}'| sort -k 2 -g | awk '{print $1}' | sed 's/\.info:#//g'").readlines():
+ input = basename.rstrip()
+ if (os.system(script1+' '+input) != 0):
+ print "Error plotting individual seismograms"; exit()
+ ps = ps + ' '+input+'.seis.pdf'
+
+if (os.system('pdcat -r '+ps+' flexwin_seis.pdf') != 0):
+ print "Error concatenate all seismogram plots"; exit()
+
+
Property changes on: seismo/3D/ADJOINT_TOMO/flexwin/scripts/plot_flexwin.py
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/flexwin/scripts/process_ds.py
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/process_ds.py (rev 0)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/process_ds.py 2009-10-30 21:17:57 UTC (rev 15901)
@@ -0,0 +1,213 @@
+#!/usr/bin/python
+
+# find matching data and synthetics from data and syn directories
+# 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
+
+# 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
+#
+# need process_data_new.pl, process_syn_new.pl, pad_zeros.pl and rotate.pl as well as sac_merge
+
+# data_list,syn_list -- input data/syn
+# all_syn_list -- input syn+derivative syn (.???)
+# new_data_list,new_syn_list -- new_datadir/data.c new_datadir/syn.c(hole)
+
+import os,sys,glob,getopt
+
+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
+
+# command-line options
+try:
+ opts,arg=getopt.getopt(sys.argv[1:],'d:s:p:t:mwD')
+except getopt.GetoptError:
+ print usage; exit()
+if (len(opts) == 0):
+ print usage; exit()
+
+for o,a in opts:
+ if o == '-d':
+ (datadir,dataext)=a.split(',')
+ if o == '-s':
+ (syndir,synext)=a.split(',')
+ if (synext != ''):
+ synext='.'+synext
+ 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()
+ if o == '-t':
+ bandpass=True
+ (tmin,tmax)=a.split(',')
+ 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+'/*')
+ if o == '-m':
+ merge=True
+ if o == '-w':
+ windowing=True
+ if o == '-D':
+ der_syn=True
+
+if (not os.path.isdir(datadir) or not os.path.isdir(syndir)):
+ print datadir, ' ', syndir, 'exist or not'; exit()
+if merge:
+ merge_dir=datadir+'/merge_dir/'
+ if (not os.path.isdir(merge_dir)):
+ os.mkdir(merge_dir)
+
+#channel name
+chan='LH'; sps='1.0';
+#padding zeros
+bb=-40
+if (bandpass):
+ bb=-float(tmin)*3.0
+
+######### obtain data/syn list##################
+nsta=0
+data_list=[]; syn_list=[]; all_syn_list=[]; hole_list=[]; comp_list=[]
+
+print '**** Data/syn list for the same (sta,net,LH?,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
+ else:
+ [sta,net,cmp]=tmp; hole=''
+ dat=glob.glob(datadir+'/*'+net+'.'+sta+'.'+hole+'.'+cmp+'*'+dataext)
+ syn=glob.glob(syndir+'/'+sta+'.'+net+'.'+cmp+synext)
+ if len(syn) != 1:
+ print 'no matching syn. for ', dat[0]
+ 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:
+ 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])
+ hole_list.append(hole); comp_list.append(cmp)
+ nsta=nsta+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 corresponding 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')
+ 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()
+
+########### cutting and band-passing data and synthetics #############
+new_data_list=[]; new_syn_list=[]
+cdata_list=''; csyn_list=''; rdata_list=''; rsyn_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]+' '
+
+ [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()
+ 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))
+ 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_syn_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')
+ 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')
+ 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()
+
+
+# ############### flexwin input file ##########################
+if windowing:
+ if (not bandpass):
+ print "use -t tmin,tmax before -w (flexwin input)"; exit()
+ print '**** writing flexwin input file ****'
+ string=''
+ 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:]
+ 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 (len(new_data) == 1 and len(new_syn) == 1):
+ string=string+new_data[0]+'\n'+new_syn[0]+'\n'+'MEASURE/'+sta+'.'+net+'.'+cmp+'\n'
+ else:
+ print "No exact matching for ", sta, net, cmp, hole_list[i]
+
+ f=open('INPUT_flexwin','w')
+ f.write(str(nsta)+'\n')
+ f.write(string.rstrip())
+ 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 ****'
Property changes on: seismo/3D/ADJOINT_TOMO/flexwin/scripts/process_ds.py
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_out.py
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_out.py (rev 0)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_out.py 2009-10-30 21:17:57 UTC (rev 15901)
@@ -0,0 +1,23 @@
+#!/usr/bin/python
+# similar to Carl's prepare_meas_all.pl script
+# Qinya Liu, Oct 2009, UofT
+import os,sys,glob
+
+if (len(sys.argv) != 3):
+ print "write_flexwin_out.py measure_dir out_filename"; exit()
+
+dir=sys.argv[1]; out=sys.argv[2]
+if (not os.path.isdir(dir)):
+ print 'check if '+dir+' exists or not'; exit()
+
+output=''
+files=glob.glob(dir+'/*mt*');
+nfiles=len(files);
+for file in files:
+ output=output+''.join(os.popen('cat '+file).readlines())
+
+output=str(nfiles)+'\n'+output
+f=open(out, 'w')
+f.write(output)
+f.close()
+
Property changes on: seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_out.py
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/ADJOINT_TOMO/flexwin/seismo_subs.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/seismo_subs.f90 2009-10-30 19:55:27 UTC (rev 15900)
+++ seismo/3D/ADJOINT_TOMO/flexwin/seismo_subs.f90 2009-10-30 21:17:57 UTC (rev 15901)
@@ -327,8 +327,9 @@
! write(*,*)'number of points =',npts
! DEBUG
- if (DEBUG) write(*,*) 'DEBUG : b, dt, npts ', b, dt, npts
+ if (DEBUG) write(*,'(a,f10.1,f10.4,i10)') ' DEBUG : b, dt, npts ', b, dt, npts
+
! read event and station header parameters from observation file
call getfhv('evla', evla, nerr)
call getfhv('evlo', evlo, nerr)
@@ -344,8 +345,8 @@
call getfhv('t2', S_pick, nerr)
endif
- ! CHT: why does this output as: BZN ^@BHR ^@AZ ^@
- if (DEBUG) write(*,*) kstnm, kcmpnm, knetwk
+ ! LQY fixed -- CHT: why does this output as: BZN ^@BHR ^@AZ ^@
+ if (DEBUG) write(*,*) 'DEBUG : sta,net,comp = ', kstnm(1:7), kcmpnm(1:7), knetwk(1:7)
! calculate distances and azimuths
call distaz(evla,evlo,stla,stlo,azimuth,backazimuth,dist_deg,dist_km)
More information about the CIG-COMMITS
mailing list