[cig-commits] r18466 - in seismo/3D/ADJOINT_TOMO/flexwin: . scripts test_data
liuqy at geodynamics.org
liuqy at geodynamics.org
Thu May 26 09:13:29 PDT 2011
Author: liuqy
Date: 2011-05-26 09:13:29 -0700 (Thu, 26 May 2011)
New Revision: 18466
Added:
seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_shift_weight.py
seismo/3D/ADJOINT_TOMO/flexwin/test_data/PAR_FILE
Removed:
seismo/3D/ADJOINT_TOMO/flexwin/test_data/PAR_FILE
Modified:
seismo/3D/ADJOINT_TOMO/flexwin/PAR_FILE
seismo/3D/ADJOINT_TOMO/flexwin/make_gfortran
Log:
revert RUN_BANDPASS to be true in Par_file
cp Par_file to test_data/ directory instead of making it a link.
add make test to make_gfortran
Modified: seismo/3D/ADJOINT_TOMO/flexwin/PAR_FILE
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/PAR_FILE 2011-05-25 20:13:34 UTC (rev 18465)
+++ seismo/3D/ADJOINT_TOMO/flexwin/PAR_FILE 2011-05-26 16:13:29 UTC (rev 18466)
@@ -20,7 +20,7 @@
# -------------------------------------------------------------
# period min/max for filtering
-RUN_BANDPASS = .false.
+RUN_BANDPASS = .true.
WIN_MIN_PERIOD = 50.0
WIN_MAX_PERIOD = 150.0
Modified: seismo/3D/ADJOINT_TOMO/flexwin/make_gfortran
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/make_gfortran 2011-05-25 20:13:34 UTC (rev 18465)
+++ seismo/3D/ADJOINT_TOMO/flexwin/make_gfortran 2011-05-26 16:13:29 UTC (rev 18466)
@@ -33,6 +33,9 @@
flexwin:flexwin.f90 copy_float.o ${LIB2}
$(FC) ${CFLAGS} flexwin.f90 -o flexwin copy_float.o ${LIB2} -L${TAULIBDIR} -L${SACLIBDIR} ${LIBS}
+test: flexwin.f90 ${LIB2}
+ cd test_data; ./flexwin < input.test; ./plot_seismos_gmt.sh MEASURE/ABKT.II.LHZ
+
clean:
rm -f ${LIB2} ${PROGS} *.o *.mod
Added: seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_shift_weight.py
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_shift_weight.py (rev 0)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_shift_weight.py 2011-05-26 16:13:29 UTC (rev 18466)
@@ -0,0 +1,189 @@
+#!/usr/bin/env python
+# 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
+
+# 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]
+ flexwin_summary_file new_flexwin_summary_file
+"""
+
+try:
+ opts,args=getopt.getopt(sys.argv[1:],'c:d:a:s:r:o:tC')
+except getopt.GetoptError:
+ sys.exit(usage)
+if len(args) != 2:
+ sys.exit(usage)
+
+# 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.
+adj_dir='.'; write_tshift=False; write_adj_dir=False; socal=False
+
+#opts
+for o,a in opts:
+ if o == '-c':
+ tmp=a.split('/')
+ if len(tmp)!=3:
+ sys.exit('-c Z/R/T')
+ 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')
+ 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)
+ if o == '-r':
+ 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 o == '-o':
+ write_adj_dir=True
+ adj_dir=a.rstrip('/')
+ if not os.path.isdir(adj_dir):
+ os.mkdir(adj_dir)
+ if o == '-t':
+ write_tshift=True
+ if o == '-C':
+ socal=True
+
+infile=args[0]; outfile=args[1]
+
+if not os.path.isfile(infile):
+ sys.exit('Check if input flexwin summary file '+infile+' exist or not')
+
+print 'Component linear weight = ', comp_z_weight, comp_r_weight, comp_t_weight
+print 'Distance exp weight = ', pnl_dist_weight, rayleigh_dist_weight, love_dist_weight
+print 'Azimuth exp weight = ', az_exp_weight
+print 'correlation exp weight = ', corr_exp_weight
+
+if socal:
+ ref_dist=100 # km
+else:
+ ref_dist=60 # degree
+ vr=3.1; vl=4.1
+naz=10; daz=360/naz;
+eps=1.0e-4; pi=3.1416926
+naz_array=ones((naz))
+
+alldata=commands.getoutput('grep data '+infile).split('\n')
+if (len(alldata) == 0):
+ sys.exit('Right now I assume data directories contain the name "data"')
+output=commands.getoutput('saclst az f '+' '.join(alldata)+"|awk '{print $2}'").split('\n')
+nfile_az=len(output)
+for i in range(0,nfile_az):
+ az=float(output[i])
+ #azimuth counts
+ k = int(divmod(az,daz)[0])
+ if (k < 0 or k >= naz):
+ sys.exit('Error binning azimuth')
+ naz_array[k] = naz_array[k] + 1
+
+print 'azimuth = '+str(naz_array)
+
+f=open(infile,'r')
+nfile=int(f.readline())
+print 'nfile = '+str(nfile)
+if nfile != nfile_az:
+ sys.exit('Check nfile and nfile_az')
+g=open(outfile,'w')
+g.write(str(nfile)+'\n')
+
+for i in range(0,nfile):
+ datfile=f.readline().rstrip('\n')
+ synfile=f.readline().rstrip('\n')
+ g.write(datfile+'\n')
+ g.write(synfile+'\n')
+
+ # 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)
+ 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):
+# 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
+
+ # componennt
+ if comp[2] == 'Z':
+ comp_weight=comp_z_weight
+ elif comp[2] == 'R':
+ comp_weight=comp_r_weight
+ elif comp[2] == 'T':
+ comp_weight=comp_t_weight
+ else:
+ sys.exit('Only deal with Z/R/T components right now')
+
+ filename=sta+'.'+net+'.'+comp+'.win.qual'
+ if write_adj_dir:
+ g.write(adj_dir+'/'+sta+'.'+net+'.'+comp+'.adj.sac\n')
+
+ nwin=int(f.readline())
+ g.write(str(nwin)+'\n')
+
+ for j in range(0,nwin):
+ (tstart,tend)=f.readline().split()
+
+ # distance
+ if socal:
+ if comp[2] == 'T':
+ dist_exp_weight=love_dist_weight
+ elif nwin > 1 and j == 0:
+ dist_exp_weight=pnl_dist_weight
+ else:
+ dist_exp_weight=rayleigh_dist_weight
+ dist_weight=(dist/ref_dist)**dist_exp_weight
+ else:
+ if comp[2] == 'T':
+ if float(tend) < dist/vl:
+ dist_weight=(sin(gcarc*pi/180)/sin(ref_dist*pi/180)) ** pnl_dist_weight
+ else:
+ dist_weight=(sin(gcarc*pi/180)/sin(ref_dist*pi/180)) ** love_dist_weight
+ else:
+ if float(tend) < dist/vr:
+ dist_weight=(sin(gcarc*pi/180)/sin(ref_dist*pi/180)) ** pnl_dist_weight
+ else:
+ dist_weight=(sin(gcarc*pi/180)/sin(ref_dist*pi/180)) ** rayleigh_dist_weight
+
+ # azimuth
+ nn=3+j;
+ [istat,out]=commands.getstatusoutput("awk 'NF == 6 {print $2,$3,$4,$5}' "+measure_dir+'/'+filename)
+ if (istat != 0 or len(out)== 0):
+ sys.exit('Error with '+"awk 'NF == 6 {print $2,$3,$4,$5}' "+measure_dir+'/'+filename)
+
+ set=False
+ for line in out.split('\n'):
+ [t1,t2,tt,cc]=line.split()
+ if abs(float(t1)-float(tstart)) < eps and abs(float(t2)-float(tend)) < eps:
+ tshift=tt; corr=cc; set=True; break
+ if not set:
+ sys.exit('can not set tshift for '+filename)
+
+ # correlation
+ k = int(divmod(az,daz)[0])
+ corr_weight = float(corr)**corr_exp_weight
+
+ weight=comp_weight*dist_weight*corr_weight/(naz_array[k]**az_exp_weight)
+# print k, weight, comp_weight, dist_weight, corr_weight, naz_array[k]**az_exp_weight
+
+ if write_tshift:
+ g.write('%s %s %s %f\n' % (tstart, tend, tshift, weight))
+ else:
+ g.write('%s %s %s %f\n' % (tstart, tend, '0.0', weight))
+
+# print weight
+
+
+f.close()
+g.close()
Property changes on: seismo/3D/ADJOINT_TOMO/flexwin/scripts/write_flexwin_shift_weight.py
___________________________________________________________________
Name: svn:executable
+ *
Deleted: seismo/3D/ADJOINT_TOMO/flexwin/test_data/PAR_FILE
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/test_data/PAR_FILE 2011-05-25 20:13:34 UTC (rev 18465)
+++ seismo/3D/ADJOINT_TOMO/flexwin/test_data/PAR_FILE 2011-05-26 16:13:29 UTC (rev 18466)
@@ -1 +0,0 @@
-link ../PAR_FILE
\ No newline at end of file
Added: seismo/3D/ADJOINT_TOMO/flexwin/test_data/PAR_FILE
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/test_data/PAR_FILE (rev 0)
+++ seismo/3D/ADJOINT_TOMO/flexwin/test_data/PAR_FILE 2011-05-26 16:13:29 UTC (rev 18466)
@@ -0,0 +1,72 @@
+# -------------------------------------------------------------
+#
+# This is the parameter file for FLEXWIN. It is based on the
+# same syntax as the Par_file for SPECFEM. Variable names are
+# put first, values are placed after the 34th column.
+#
+# Comment lines and blank lines are significant. If you
+# change the layout of this file or add/remove parameters
+# you must also modify the user_variables module and the
+# read_parameter_file subroutine at the start of seismo_subs.f90.
+#
+# -------------------------------------------------------------
+
+# -------------------------------------------------------------
+# boolean parameters
+DEBUG = .true.
+MAKE_SEISMO_PLOTS = .true.
+MAKE_WINDOW_FILES = .false.
+BODY_WAVE_ONLY = .false.
+
+# -------------------------------------------------------------
+# period min/max for filtering
+RUN_BANDPASS = .true.
+WIN_MIN_PERIOD = 50.0
+WIN_MAX_PERIOD = 150.0
+
+# -------------------------------------------------------------
+# E(t) water level
+STALTA_BASE = 0.08
+
+# -------------------------------------------------------------
+# maximum allowable time shift from reference TSHIFT
+TSHIFT_BASE = 15.0
+TSHIFT_REFERENCE = 0.0
+
+# -------------------------------------------------------------
+# maximum allowable amplitude measurement relative to reference DLNA
+DLNA_BASE = 1.0
+DLNA_REFERENCE = 0.0
+
+# -------------------------------------------------------------
+# limit on CC for window acceptance
+CC_BASE = 0.85
+
+# -------------------------------------------------------------
+# boolean switch for check_data_quality
+DATA_QUALITY = .true.
+
+# if DATA_QUALITY = .true. and if two different measurements of
+# signal-to-noise ratios exceeds these two base levels,
+# then the data time series (and syn) is kept
+SNR_INTEGRATE_BASE = 3.5
+SNR_MAX_BASE = 3.0
+
+# -------------------------------------------------------------
+# limit on signal to noise ratio in a particular window.
+WINDOW_S2N_BASE = 1.5
+
+# -------------------------------------------------------------
+# Fine tuning constants
+C_0 (internal minima) = 0.7
+C_1 (small windows) = 4.0
+C_2 (prominence) = 0.0
+C_3a (separation height) = 1.0
+C_3b (separation time) = 2.0
+C_4a (curtail on left) = 3.0
+C_4b (curtail on right) = 10.0
+
+WEIGHT_SPACE_COVERAGE = 1.0
+WEIGHT_AVERAGE_CC = 1.0
+WEIGHT_N_WINDOWS = 1.0
+
More information about the CIG-COMMITS
mailing list