[cig-commits] r16350 - seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process
bozdag at geodynamics.org
bozdag at geodynamics.org
Fri Feb 26 13:47:31 PST 2010
Author: bozdag
Date: 2010-02-26 13:47:31 -0800 (Fri, 26 Feb 2010)
New Revision: 16350
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_data.pl
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_syn.pl
Log:
Updated the data processing scripts in UTILS/seis_process
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_data.pl
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_data.pl 2010-02-26 21:04:48 UTC (rev 16349)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_data.pl 2010-02-26 21:47:31 UTC (rev 16350)
@@ -7,15 +7,20 @@
sub Usage{
print STDERR <<END;
- process_data_new.pl -m CMTFILE -l Start/End -t Ts/Tl -i -p -x Ext -d OutDir
- -s Sps -P n/p -T Taper_width -c -y t1/v1/...
+ process_data_new.pl -m CMTFILE -l Start/End -t Ts/Tl -f -i -R resp_dir
+ -p -s Sps -P n/p -T Taper_width -c -y t1/v1/...
+ -x Ext -d OutDir
data_sac_files
with
-m -- use CMTSOLUTION file to set event info
+ -a STAFILE -- add station information from STAFILE, use '-a none' to get
+ information from Vala''s station file.
-l Start/End -- start and end of record from o
-t Ts/Tl -- shortest and longest period
+ -f -- do not apply butterworth filter before deconvolution
-i -- transfer record to displacement
+ -R resp_dir -- transfer record using RESP file in resp_dir
-p -- pick event p and s first arrival into sac headers (t1 and t2)
-x Ext -- extension of new file name
-d Dir -- directory to put the output files (default .)
@@ -28,24 +33,32 @@
data_sac_files --- names of the data sac files to be processed
Notice:
- 0. Make sure you have sac, saclst, phtimes in PATH before execution
- 1. We require that polezero files are in the same directory as the sac
- data files which is generally satisfied. All the needed info are
+ 1. We require that polezero files in the same directory as the sac
+ data files which is generally satisfied. We require that resp dir
+ be specified even if it is current(.). All the needed info is
taken from sac headers
2. Origin time is set to PDE + time_shift (given by the CMTSOLUTION)
3. The displacement outputs after -i option are in the unit of meters
4. For BH? components, please set -s 20, otherwise interpolation to
1 sample/second will be performed
- Qinya Liu, May 2007, Caltech
-
+ NOTE: Please make sure that SAC, saclst and IASP91 packages are installed properly on
+ your system, and that all the environment variables are set properly before
+ running the script.
+
+ Qinya Liu, Originally written in Oct 2002; updated in Feb 2010
END
exit(1);
}
if (@ARGV == 0) { Usage(); }
-if (!getopts('m:l:t:ipx:d:s:P:T:y:c')) {die('Check input arguments\n');}
+if (!getopts('m:a:l:t:fiR:px:d:s:P:T:y:c')) {die('Check input arguments\n');}
+
+if ($opt_m and not -f $opt_m) {die("Check if file $opt_m exists\n");}
+if ($opt_a and not -f $opt_a) {$opt_a="/opt/seismo-util/data/STATIONS_new";}
+if ($opt_R and not -d $opt_R) {die("Check if dir $opt_R exists\n");}
+
if ($opt_t) {($tmin, $tmax) = split(/\//,$opt_t);
$f1 = 1./$tmax; $f2=1./$tmin;}
if ($opt_d) {
@@ -60,14 +73,15 @@
if (!$opt_T) {$opt_T = 0.05;}
if ($opt_l) {($lmin,$lmax) = split(/\//,$opt_l);} else {$lmin = 0; $lmax = 3600;}
-$saclst="/opt/seismo-util/source/saclst/saclst";
-$phtimes="/opt/seismo-util/bin/phtimes";
-if (! -e $saclst) {die(" No $saclst file\n");}
-if (! -e $phtimes) {die("No $phtimes file\n");}
+if ($opt_f and not ($opt_R or $opt_i)) {die("-f option goes together with -i or -R\n");}
+$saclst="saclst";
+$phtimes="phtimes";
+
$eps=1e-5; $undef=-12345.0; $cundef="-12345";
foreach $file (@ARGV) {
+ $no_resp=0;
if (! -f $file) {die("No such file to be processed!!\n"); }
print "Processing data file $file \n";
if (not $opt_d) {$outfile = $file.$ext;}
@@ -82,8 +96,7 @@
print SAC "echo on\n";
print SAC "r $outfile\n";
- # add event information from CMTSOLUTION FILE
- if ($opt_m) {
+ if ($opt_m) { # add event information from CMTSOLUTION FILE
print " Adding event info from $opt_m\n";
($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift,undef,
$elat,$elon,$edep,undef) = get_cmt($opt_m);
@@ -100,7 +113,29 @@
print SAC "echo on\n";
print SAC "r $outfile\n";}
- if ($opt_l){ # cut record and rtrend and rmean
+ if ($opt_a) { # add station information
+ print "Add station information\n";
+ (undef,undef,undef,undef,undef,undef,$net,$sta,undef,$comp)=split(/\./,$file);
+ print SAC "ch kstnm $sta knetwk $net kcmpnm $comp\n ";
+ ($sta_name,$sta_net,$sta_lat,$sta_lon,$sta_ele,$sta_bur)=split(" ",`egrep '$sta +$net' $opt_a`);
+ if (not defined $sta_name) {
+ print("No such station $sta+$net in $opt_a file\n");
+ ($sta_name,$sta_net,$sta_lat,$sta_lon,$sta_ele,$sta_bur)=split(" ",`egrep '$sta' $opt_a`);
+ if (not defined $sta_name) {
+ print " No such station as $sta in the $opt_a file\n";
+ $sta_text .="$sta ";}}
+ print SAC "ch stla $sta_lat stlo $sta_lon stel $sta_ele stdp $sta_bur\nwh\n";
+ print SAC "w $outfile\nquit\n";
+ close(SAC);
+ if ($opt_c) {open(SAC,"|sac ") || die("Can't open sac\n");}
+ else {open(SAC,"|sac > /dev/null") || die("Can't open sac\n");}
+ print SAC "echo on\n";
+ print SAC "r $outfile\n";}
+
+ if ($opt_s) {print SAC "interp delta $dt\n";
+ print SAC "w over\n";}
+
+ if ($opt_l){ # cut record
print " Cut the record from o+$lmin to o+$lmax\n";
(undef,$tmp_o)=split(" ",`$saclst o f $outfile`);
if (abs($tmp_o-$undef) < $eps) {die("O has not been set to cut record\n");}
@@ -114,28 +149,39 @@
print SAC "echo on\n";
print SAC "r $outfile\n";}
- if ($opt_s) {
- print SAC "interp delta $dt\n";}
- if ($opt_t) {# filter records
+ if ($opt_t and not $opt_f) {# filter records
print " Filter record at periods $tmin and $tmax\n";
print SAC "rtrend\n rmean\n taper width $opt_T\n";
printf SAC ("bp n %4d p $pass cor %12.6f %12.6f\n",$poles,$f1,$f2);
printf SAC " rtrend\n rmean\n taper width $opt_T\n";}
- if ($opt_i) { # transfer instrument response
+ if ($opt_i or $opt_R) { # transfer instrument response
if (! $opt_t) {die(" Specify bandpass range by -t\n");}
$f0=$f1*0.8;
$f3=$f2*1.2;
+ print SAC "rtrend\n rmean\n taper width $opt_T\n";
(undef,$network,$sta,$comp,$khole)=split(" ",`$saclst knetwk kstnm kcmpnm khole f $outfile`);
if ($network eq $cundef or $sta eq $cundef or $comp eq $cundef or $khole eq $cundef) {
die("No network station name or component name or khole defined in $outfile\n");}
- $pzfile=`ls -1 $filedir/SAC_PZs_${network}_${sta}_${comp}_${khole}* | head -1`;
- chomp($pzfile);
- if (!-f $pzfile) {die(" No pzfile $pzfile \n");}
- print " Transfer instrument response using pz file $pzfile\n";
- printf SAC ("trans from polezero s %20s to none freq %12.6f%12.6f%12.6f%12.6f \n",
- $pzfile,$f0,$f1,$f2,$f3);}
+ if ($opt_i) { # assume pz file in the same dir
+ $pzfile=`ls -1 $filedir/SAC_PZs_${network}_${sta}_${comp}_${khole}* | head -1`;
+ chomp($pzfile);
+ if (!-f $pzfile) {print " **** No pzfile $pzfile **** \n"; $no_resp=1;}
+ else {print " Transfer instrument response using pz file $pzfile\n";
+ printf SAC ("trans from polezero s %20s to none freq %12.6f%12.6f%12.6f%12.6f \n",
+ $pzfile,$f0,$f1,$f2,$f3);}
+ } else { # assume resp file in opt_R
+ @respfiles=glob("$opt_R/${network}.${sta}.${khole}.${comp}*.RESP");
+ $respfile = find_resp_file($outfile, at respfiles);
+ if (! -f $respfile or @respfiles == 0) {print " **** No respfile $respfiles for $network, $sta, $khole, $comp ****\n"; $no_resp=1;}
+ else {
+ print " Transfer instrument response using response file $respfile\n";
+ printf SAC ("trans from evalresp fname $respfile to none freq %12.6f%12.6f%12.6f%12.6f \n",
+ $f0,$f1,$f2,$f3);
+ printf SAC "mul 1e-9\n"; }}
+ printf SAC " rtrend\n rmean\n taper width $opt_T\n";
+ }
if ($opt_p) { # add p and s arrival info
print " Add first P and S arrival information\n";
@@ -166,6 +212,9 @@
print SAC "w $outfile\n";
print SAC "echo off\nquit\n";
close(SAC);
+ if ($no_resp and defined $opt_d) {
+ print "Deleting $outfile\n";
+ system("rm -f $outfile");}
}
print " Done! \n";
@@ -214,3 +263,21 @@
close(CMT);
return ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift,$hdur,$elat,$elon,$edep,$Mrr,$Mtt,$Mpp,$Mrt,$Mrp,$Mtp);
}
+
+sub find_resp_file {
+ my($file, at resp)=@_;
+ # using sac own saclst
+ my(undef,$year,$jday,$hr,$min,$sec,$msec)=split(" ",`saclst nzyear nzjday nzhour nzmin nzsec nzmsec f $file`);
+ $f_date=sprintf("%04d.%03d.%02d.%02d.%02d.%04d",$year,$jday,$hr,$min,$sec,$msec);
+# print "$f_date\n";
+ my($resp_file)="";
+ foreach $resp (@resp) {
+ ($rs_date) = split(" ", `grep 'Start date' $resp | awk '{print \$4}' | tr , . | tr : .`);
+ ($re_date) = split(" ", `grep 'End date' $resp | awk '{print \$4}' | tr , . | tr : .`);
+
+# print "$rs_date; $re_date\n";
+ if ($rs_date le $f_date and $re_date ge $f_date) {
+ $resp_file=$resp; last;}
+ }
+ return($resp_file);
+}
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_syn.pl
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_syn.pl 2010-02-26 21:04:48 UTC (rev 16349)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_syn.pl 2010-02-26 21:47:31 UTC (rev 16350)
@@ -8,7 +8,7 @@
print STDERR <<END;
Usage: process_syn_new.pl
- -S -m CMTFILE -h -o offset -lStart/End -tTmin/Tmax -P n/p
+ -S -m CMTFILE -h -o offset -lStart/End -tTmin/Tmax -f -P n/p
-i Dir -A amp_factor -p -a STAFILE -s sps -y t0/V0/t1/V1...
-c -d OutDir -x Ext synthetics_files
where
@@ -18,6 +18,7 @@
-o moves the synthetics back by 'offset' from centroid time.
-l lmin/lmax start/end cut of trace
-t Tmin/Tmax -- shortest/longest period of bandpass filter
+ -f -- apply transfer filter instead of butterworth filter for -t
-P -- number of poles and passes in butterworth filter(default 4/2)
-i Dir -- convolve synthetics with instrument response in Dir
-A amp_factor -- multiple all files with an amplification factor
@@ -32,10 +33,12 @@
names of files -- name of syn files to be processed
- Make sure you have sac, saclst, phtimes, asc2sac, ascii2sac.csh
- in PATH before execution
+ NOTE: Please make sure that SAC, saclst and IASP91 packages are installed properly on
+ your system, and that all the environment variables are set properly before
+ running the script.
- Qinya Liu, Caltech, May 2007
+
+ Qinya Liu, originally written in Oct 2002, most recently updated in Feb 2010.
END
exit(1);
@@ -43,7 +46,7 @@
@ARGV > 1 or Usage();
-if (!getopts('Sm:ho:l:t:i:pP:a:cd:x:vs:y:A:')) {die(" check input arguments\n");}
+if (!getopts('Sm:ho:l:t:fi:pP:a:cd:x:vs:y:A:')) {die(" check input arguments\n");}
if ($opt_t) {($tmin, $tmax) = split(/\//,$opt_t);
$f1 = 1./$tmax; $f2=1./$tmin;}
@@ -54,21 +57,24 @@
if(not defined $pass or $pass<1){$pass=2;}}
if ($opt_l) {($lmin,$lmax) = split(/\//,$opt_l);}
else {$lmin = 0; $lmax = 3600;}
-if ($opt_a and not -f $opt_a) {$opt_a="/opt/seismo-util/data/STATIONS_new";}
+if ($opt_a and not -f $opt_a) {$opt_a="/opt/seismo/data/STATIONS_new";}
if ($opt_o and not $opt_m) {die("Specify centroid time first\n");}
if ($opt_d and not -d $opt_d) {die("No such directory as $opt_d\n");}
if ($opt_i and not -d $opt_i) {die("No such directory as $opt_i\n");}
if ($opt_A) {if ($opt_A !~ /^\d/) {die("-A option should be numbers\n");}}
+if ($opt_f) {
+ if (not $opt_t) {die("-t tmin/tmax has to be specified for -f\n");}
+ else {$f0=$f1*0.8; $f3=$f2*1.2;}}
$undef = -12345;
$eps = 0.1;
-$saclst="/opt/seismo-util/source/saclst/saclst";
-$phtimes="/opt/seismo-util/bin/phtimes";
-$asc2sac="/opt/seismo-util/bin/ascii2sac.csh";
-if (! -e $saclst) {die(" No $saclst file\n");}
-if (! -e $phtimes) {die("No $phtimes file\n");}
-if (! -e $asc2sac) {die("No $asc2sac file\n");}
+$saclst="saclst";
+$phtimes="phtimes";
+$asc2sac="ascii2sac.csh";
+#if (! -e $saclst) {die(" No $saclst file\n");}
+#if (! -e $phtimes) {die("No $phtimes file\n");}
+#if (! -e $asc2sac) {die("No $asc2sac file\n");}
$min_hdur=1.0;
@@ -115,9 +121,13 @@
if ($opt_c) {open(SAC,"|sac");}else {open(SAC,"|sac > /dev/null");}
print SAC "echo on\n r $outfile\n";}
- if ($opt_h and $hdur > $min_hdur) { # convolve source time function
- system("/opt/seismo-util/bin/convolve_stf g $hdur $outfile");
- system("mv ${outfile}.conv ${outfile}");}
+ if ($opt_h and $hdur > $min_hdur) { # convolve source time function
+ print SAC "quit\n";
+ close(SAC);
+ system("/opt/seismo-utils/bin/convolve_stf g $hdur $outfile");
+ system("mv ${outfile}.conv ${outfile}");
+ if ($opt_c) {open(SAC,"|sac");}else {open(SAC,"|sac > /dev/null");}
+ print SAC "echo on\n r $outfile\n";}
if ($opt_a) { # add station information
@@ -136,7 +146,13 @@
$sta_text .="$sta ";}}
print SAC "ch stla $sta_lat stlo $sta_lon stel $sta_ele stdp $sta_bur\nwh\n";}
- if ($opt_l) {
+ if ($opt_s) {print SAC "interp delta $dt\n";}
+ print SAC "w over\nquit\n";
+ close(SAC);
+ if ($opt_c) { open(SAC,"|sac");}else {open(SAC,"|sac > /dev/null");}
+ print SAC "echo on\n r $outfile\n ";
+
+ if ($opt_l) { #cut record
print SAC "setbb begin ( max &1,b ( &1,o + $lmin ) ) \n";
print SAC "setbb end ( min &1,e ( &1,o + $lmax ) ) \n";
print SAC "cut %begin% %end% \n";
@@ -145,18 +161,15 @@
if ($opt_A) {print SAC "mul $opt_A\n";}
- if ($opt_s) {print SAC "interp delta $dt\n";}
- print SAC "w over\nquit\n";
- close(SAC);
- if ($opt_c) { open(SAC,"|sac");}else {open(SAC,"|sac > /dev/null");}
- print SAC "echo on\n r $outfile\n ";
-
- if ($opt_t){ # filter record
+ if ($opt_t) { # filter record
print "Filtering records...\n";
print SAC "rtrend\n rmean\n taper\n";
- printf SAC ("bp n $poles p $pass co %10.5f %10.5f\n",$f1,$f2);
+ if ($opt_f) {
+ printf SAC ("trans from none to none freq %12.6f %12.6f %12.6f %12.6f\n",$f0,$f1,$f2,$f3);}
+ else {
+ printf SAC ("bp n $poles p $pass co %10.5f %10.5f\n",$f1,$f2);}
print SAC "rtrend\n rmean\n taper\n";}
-
+
if ($opt_i) {# convolve with instrument response
print "Convolving instrument response...\n";
$pzfile="SAC_PZs_${net}_${sta}_${comp}_";
@@ -247,4 +260,3 @@
return ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift,$hdur,$elat,$elon,$edep,$Mrr,$Mtt,$Mpp,$Mrt,$Mrp,$Mtp);
}
-
More information about the CIG-COMMITS
mailing list