[cig-commits] r22274 - seismo/3D/GRD_CMT3D/scripts
liuqy at geodynamics.org
liuqy at geodynamics.org
Fri Jun 14 14:46:13 PDT 2013
Author: liuqy
Date: 2013-06-14 14:46:13 -0700 (Fri, 14 Jun 2013)
New Revision: 22274
Added:
seismo/3D/GRD_CMT3D/scripts/check_ds_exist.bash
seismo/3D/GRD_CMT3D/scripts/diff_meas.bash
seismo/3D/GRD_CMT3D/scripts/mech_table3.pl
Modified:
seismo/3D/GRD_CMT3D/scripts/cmt.py
seismo/3D/GRD_CMT3D/scripts/cp_data.bash
seismo/3D/GRD_CMT3D/scripts/cp_syn.bash
seismo/3D/GRD_CMT3D/scripts/mech_table.pl
seismo/3D/GRD_CMT3D/scripts/mech_table2.pl
Log:
update some scripts.
Added: seismo/3D/GRD_CMT3D/scripts/check_ds_exist.bash
===================================================================
--- seismo/3D/GRD_CMT3D/scripts/check_ds_exist.bash (rev 0)
+++ seismo/3D/GRD_CMT3D/scripts/check_ds_exist.bash 2013-06-14 21:46:13 UTC (rev 22274)
@@ -0,0 +1,16 @@
+#!/bin/bash
+# this script checks if all the files listed in MEASURE*ALL file exist or not
+if [ $# != 1 ]; then
+ echo "check_ds_exist.bash event-list"
+ exit
+fi
+
+for evid in `cat $1 `; do
+cd $evid
+for file in ` cat MEASURE*ALL | grep '^[ds]' `; do
+ if [ ! -f $file ]; then
+ echo $file
+ fi
+done
+cd ..
+done
Property changes on: seismo/3D/GRD_CMT3D/scripts/check_ds_exist.bash
___________________________________________________________________
Added: svn:executable
+ *
Modified: seismo/3D/GRD_CMT3D/scripts/cmt.py
===================================================================
--- seismo/3D/GRD_CMT3D/scripts/cmt.py 2013-06-14 21:33:52 UTC (rev 22273)
+++ seismo/3D/GRD_CMT3D/scripts/cmt.py 2013-06-14 21:46:13 UTC (rev 22274)
@@ -287,6 +287,7 @@
# ====================================
# extension name array is just a name holder, the deriatives may differ from the names they are called.
+# dlocation mean different things for utm = 0, and utm in [1,60], and dlocation should be set to the same value as ddepth for utm=-1
def gen_cmt_der(cmt,npar=9,dmoment=1e22,ddepth=1,dlocation=0.05,utm=0,nevent=1,ext=['Mrr','Mtt','Mpp','Mrt','Mrp','Mtp','dep','lon','lat']):
# ddepth and dlocation are in terms of kms
Modified: seismo/3D/GRD_CMT3D/scripts/cp_data.bash
===================================================================
--- seismo/3D/GRD_CMT3D/scripts/cp_data.bash 2013-06-14 21:33:52 UTC (rev 22273)
+++ seismo/3D/GRD_CMT3D/scripts/cp_data.bash 2013-06-14 21:46:13 UTC (rev 22274)
@@ -1,4 +1,5 @@
#!/bin/bash
+# this is an old program used to copy tar'ed derivative synthetics from caltech machines to local machines.
if [ $# != 1 ]; then
echo "Usage: cp_data.bash evid-file"; exit
@@ -9,9 +10,10 @@
for evid in `cat $1`; do
mkdir -p $evid/data $evid/backup; cd $evid/data
echo " == $evid ==="
- ssh wagholi.gps.caltech.edu "cd $dir/$evid; tar cjvf ~/data_${evid}.tar.bz *.sac >/dev/null "
+ ssh -l lqy wagholi.gps.caltech.edu "cd $dir/$evid; tar cjvf ~/data_${evid}.tar.bz *.sac >/dev/null "
scp lqy at wagholi.gps.caltech.edu:data_${evid}.tar.bz .
- ssh wagholi.gps.caltech.edu "rm -f data_${evid}.tar.bz"
- tar xjvf data_${evid}.tar.bz; mv data_${evid}.tar.bz ../backup
+# ssh wagholi.gps.caltech.edu "rm -f data_${evid}.tar.bz"
+ mv ../../data/data_${evid}.tar.bz .
+ tar xjvf data_${evid}.tar.bz > /dev/null; mv data_${evid}.tar.bz ../backup
cd ../..
done
Modified: seismo/3D/GRD_CMT3D/scripts/cp_syn.bash
===================================================================
--- seismo/3D/GRD_CMT3D/scripts/cp_syn.bash 2013-06-14 21:33:52 UTC (rev 22273)
+++ seismo/3D/GRD_CMT3D/scripts/cp_syn.bash 2013-06-14 21:46:13 UTC (rev 22274)
@@ -1,5 +1,7 @@
#!/bin/bash
-
+# this is an old script used to copy and untar derivative syns,
+# copy measurement files, add up syns, change file extension,
+# copy inversion pars, figure out the dlat, dlon from CMTs.
if [ $# != 1 ]; then
echo "Usage: cp_syn.bash evid-file"; exit
fi
Added: seismo/3D/GRD_CMT3D/scripts/diff_meas.bash
===================================================================
--- seismo/3D/GRD_CMT3D/scripts/diff_meas.bash (rev 0)
+++ seismo/3D/GRD_CMT3D/scripts/diff_meas.bash 2013-06-14 21:46:13 UTC (rev 22274)
@@ -0,0 +1,11 @@
+# diff measurement files
+for freq in 2 3 6; do
+for evid in `cat event_list`; do
+ res1=`head -n 1 $evid/MEASUREMENT_WINDOWS_${evid}_T00${freq}_T030_m12 `
+ res2=`head -n 1 event_info/MEASUREMENT_WINDOWS_${evid}_T00${freq}_T030_m12`
+ if [ $res1 -ne $res2 ]; then
+ echo $evid -- $freq --- $res1 --- $res2
+ fi
+ done
+done
+
Property changes on: seismo/3D/GRD_CMT3D/scripts/diff_meas.bash
___________________________________________________________________
Added: svn:executable
+ *
Modified: seismo/3D/GRD_CMT3D/scripts/mech_table.pl
===================================================================
--- seismo/3D/GRD_CMT3D/scripts/mech_table.pl 2013-06-14 21:33:52 UTC (rev 22273)
+++ seismo/3D/GRD_CMT3D/scripts/mech_table.pl 2013-06-14 21:46:13 UTC (rev 22274)
@@ -4,9 +4,9 @@
use GMT_PLOT_SC;
use CMT_TOOLS;
use POSIX;
+
-
-# parameter files
+# set up different parameter settings
$inv_file="INVERSION.PAR.SAVE";
$outdir="mech"; $psfile="mech.ps";
@@ -50,7 +50,8 @@
print BASH "cp -f $cmt $outdir\n";
if (-f "${cmt}_GRD") {print BASH "cp -f ${cmt}_GRD $outdir\n";}
-for ($k=2;$k<$nbeach;$k++) { # cmt3d runs
+# cmt3d runs for different INVERSION.PAR
+for ($k=2;$k<$nbeach;$k++) {
print BASH "echo Running cmt3d_flexwin for INVERSION.$k ...\n";
open(INV,">INVERSION.$k");
print INV "${cmt}\n${cmt_new}\n$npars[$k]\n$delta${flex}.true.\n$wdatas[$k]\n$scorrs[$k]\n$cons[$k]\n${wns}";
@@ -72,7 +73,7 @@
$var[$k]=$tmp[8]; print "$k -- VR = $var[$k]\n";}
#die("here\n");
-
+# plot results
open(BASH,">mech_plot.bash");
print BASH "gmtset BASEMAP_TYPE plain ANOT_FONT_SIZE 9 HEADER_FONT_SIZE 10 MEASURE_UNIT inch PAPER_MEDIA letter TICK_LENGTH 0.1c\n";
@@ -87,9 +88,11 @@
if (-f "$outdir/$cmt[$k]") {
plot_psmeca_raw(\*BASH,$psfile,"-JX -R -Sm1.0 -B3/3wesn","1.5 2","$outdir/$cmt[$k]");
@output = `cmtsol2faultpar.pl $outdir/$cmt[$k]`;
- (undef,$tmp1) = split(" ", $output[4]);
- (undef,$tmp2) = split(" ", $output[6]);
+ (undef,$tmp1) = split(" ", $output[5]);
+ (undef,$tmp2) = split(" ", $output[7]);
+
($s,$d,$r,$p,$t) = split(/\//,$tmp2);
+
$tex="1 1 9 0 4 CM $names[$k]\n";
$tex.="0.2 0.8 9 0 4 LM Mw/dep/eps/= $tmp1\n0.2 0.6 9 0 4 LM S/D/R= $s/$d/$r";
if ($k>=2) {$tex.=sprintf("\n0.2 0.4 9 0 4 LM Var=%6.2f%",$var[$k]);}
Modified: seismo/3D/GRD_CMT3D/scripts/mech_table2.pl
===================================================================
--- seismo/3D/GRD_CMT3D/scripts/mech_table2.pl 2013-06-14 21:33:52 UTC (rev 22273)
+++ seismo/3D/GRD_CMT3D/scripts/mech_table2.pl 2013-06-14 21:46:13 UTC (rev 22274)
@@ -4,6 +4,7 @@
use GMT_PLOT_SC;
use CMT_TOOLS;
use POSIX;
+use List::Util qw[max];
## ok, this version includes both 2, 3 and 6 seconds measurements
# parameter files
@@ -65,6 +66,7 @@
print INV "${cmt}\n${cmt_new}\n$npars[$k]\n$delta$meas[$k]\n.true.\n$wdatas[$k]\n$scorrs[$k]\n$cons[$k]\n";
if ($k==2) {print INV ".true.\n";} else { print INV "${wns}\n";}
close(INV);
+
print BASH "# --- mech $k ----\n";
print BASH "cp -f INVERSION.$k INVERSION.PAR\n";
print BASH "cmt3d_flexwin > cmt3d.stdout\n";
@@ -76,13 +78,21 @@
print BASH "mv -f INVERSION.$k $outdir\n";
}
close(BASH);
+
system("chmod a+x mech_table.bash; mech_table.bash");
for ($k=2;$k<$nbeach;$k++) {
(@tmp)=split(" ",`grep Variance $outdir/cmt3d_stdout.$k`);
$var[$k]=$tmp[8]; print "$k -- VR = $var[$k]\n";}
+
#die("here\n");
+# get depth information to display
+($odep)=split(" ",`cmtsol2faultpar.pl $cmt | grep mw | cut -d = -f 2 | cut -d / -f 2`);
+($ndep)=split(" ",`cmtsol2faultpar.pl $outdir/$cmt[2] | grep mw | cut -d = -f 2 | cut -d / -f 2`);
+($dep1)=split(" ",`cmtsol2faultpar.pl $outdir/$cmt[6] | grep mw | cut -d = -f 2 | cut -d / -f 2`);
+($dep2)=split(" ",`cmtsol2faultpar.pl $outdir/$cmt[7] | grep mw | cut -d = -f 2 | cut -d / -f 2`);
+($dep3)=split(" ",`cmtsol2faultpar.pl $outdir/$cmt[8] | grep mw | cut -d = -f 2 | cut -d / -f 2`);
+$ddep=sprintf("%5.2f",max(abs($dep1-$ndep),abs($dep2-$ndep),abs($dep3-$ndep)));
-
open(BASH,">mech_plot.bash");
print BASH "gmtset BASEMAP_TYPE plain ANOT_FONT_SIZE 9 HEADER_FONT_SIZE 10 MEASURE_UNIT inch PAPER_MEDIA letter TICK_LENGTH 0.1c\n";
@@ -97,8 +107,8 @@
if (-f "$outdir/$cmt[$k]") {
plot_psmeca_raw(\*BASH,$psfile,"-JX -R -Sm1.0 ","1.5 2","$outdir/$cmt[$k]");
@output = `cmtsol2faultpar.pl $outdir/$cmt[$k]`;
- (undef,$tmp1) = split(" ", $output[4]);
- (undef,$tmp2) = split(" ", $output[6]);
+ (undef,$tmp1) = split(" ", $output[5]);
+ (undef,$tmp2) = split(" ", $output[7]);
($s,$d,$r,$p,$t) = split(/\//,$tmp2);
$tex="1 1 9 0 4 CM $names[$k]\n";
$tex.="0.2 0.8 9 0 4 LM Mw/dep/eps/= $tmp1\n0.2 0.6 9 0 4 LM S/D/R= $s/$d/$r";
@@ -113,7 +123,7 @@
($n3) = split(" ",`grep 'T003.*.sac.d' $meas[2] | wc | awk '{print \$1}'`);
($n6) = split(" ",`grep 'T006.*.sac.d' $meas[2] | wc | awk '{print \$1}'`);
- plot_pstext(\*BASH,$psfile,"-JX -R -N","1.5 4 12 0 4 CM $ename ($nz/$nr/$nt; $n2/$n3/$n6)");}
+ plot_pstext(\*BASH,$psfile,"-JX -R -N","0.5 4 12 0 4 CM $ename ($nz/$nr/$nt; $n2/$n3/$n6; $odep/$ndep/$ddep)");}
plot_psxy(\*BASH,$psfile,"-JX -X-$x -Y-$y","");
$k++;
}
Added: seismo/3D/GRD_CMT3D/scripts/mech_table3.pl
===================================================================
--- seismo/3D/GRD_CMT3D/scripts/mech_table3.pl (rev 0)
+++ seismo/3D/GRD_CMT3D/scripts/mech_table3.pl 2013-06-14 21:46:13 UTC (rev 22274)
@@ -0,0 +1,131 @@
+#!/usr/bin/perl
+use lib '/opt/seismo/lib/perl';
+use GMT_PLOT;
+use GMT_PLOT_SC;
+use CMT_TOOLS;
+use POSIX;
+use List::Util qw[max];
+
+## ok, this version includes both 2, 3 and 6 seconds measurements
+# parameter files
+$inv_file="INVERSION.PAR.SAVE";
+$outdir="mech"; $psfile="mech.ps";
+
+$weigh_data="2 2 1 0.5 1.15 0.55 0.78";
+$weigh_data_nocomp="1 1 1 0.5 1.15 0.55 0.78";
+$weigh_data_noaz="2 2 1 0 1.15 0.55 0.78";
+
+if (not -f "CMTSOLUTION") {die("Check cmt solution file\n");}
+($evid)=split(" ",`cat CMTSOLUTION | awk 'NR == 2 {print \$3}'`);
+$mwin="MEASUREMENT_WINDOWS";
+
+ at names=("Initial","Grid Search","7 Par+ZT",
+ "7 Par+no-SC","6 par+DC","7 Par+DC",
+ "7 Par(2-30s)","7 Par(3-30s)","7 Par(6-30s)");
+ at npars=(6,6,7, 7,6,7, 7,7,7);
+ at wdatas=($weigh_data, $weigh_data, $weigh_data, $weigh_data,$weigh_data,$weigh_data, $weigh_data,$weigh_data,$weigh_data);
+ at scorrs=("true","true",".true.", ".false.", ".true.", ".true.", ".true.",".true.",".true.");
+ at cons=(".true. true. 0.0", ".true. .true. 0.0", ".true. .false. 0.0",
+ ".true. .false. 0.0", "true. .false. 0.0", "true. .true. 0.0",
+ ".true. .false. 0.0",".true. .false. 0.0",".true. .false. 0.0");
+ at meas=("${mwin}_${evid}_ALL","${mwin}_${evid}_ALL","${mwin}_${evid}_ALL",
+ "${mwin}_${evid}_ALL","${mwin}_${evid}_ALL","${mwin}_${evid}_ALL",
+ "${mwin}_${evid}_T002_T030_m12","${mwin}_${evid}_T003_T030_m12","${mwin}_${evid}_T006_T030_m12");
+
+$nbeach=@npars;
+$ncols = 3;
+$nrows = 3;
+# plot size jx/jy, plot origin in xy
+($jx,$jy) = shift_xy("1 1","$ncols $nrows","7 7", \@xy,"0.8 0.8");
+$JX = "-JX$jx/$jy"; $R="-R0/3/0/3";
+$GMT_PLOT::paper_orient="-P";
+
+if (not -f "cmt3d_flexwin") {die("Check if cmt3d_flexwin exists or not\n");}
+
+open(INV,"$inv_file")|| die("check $inv_file\n");
+ at all=<INV>;close(INV);
+
+$cmt=$all[0]; $cmt_new=$all[1]; $npar=$all[2]; $delta=$all[3];
+$flex=$all[4];$lwdata=$all[5]; $wdata=$all[6]; $scorr=$all[7]; $con=$all[8];
+$wns=$all[9];
+chomp($cmt_new);chomp($cmt);
+(undef,undef,$ename)=split(" ",`grep 'event name' $cmt`);
+
+ at cmt=($cmt);@cmt=(@cmt,"${cmt}_GRD");
+if (-f "${cmt}_GRD") {print BASH "cp -f ${cmt}_GRD $outdir\n";}
+for ($i=2;$i<$nbeach;$i++) {@cmt=(@cmt,"CMTSOLUTION.$i");}
+
+open(BASH,">mech_table.bash");
+print BASH "mkdir -p $outdir\n";
+print BASH "cp -f $cmt $outdir\n";
+if (-f "${cmt}_GRD") {print BASH "cp -f ${cmt}_GRD $outdir\n";}
+
+for ($k=2;$k<$nbeach;$k++) { # cmt3d runs
+ print BASH "echo Running cmt3d_flexwin for INVERSION.$k ...\n";
+ open(INV,">INVERSION.$k");
+ print INV "${cmt}\n${cmt_new}\n$npars[$k]\n$delta$meas[$k]\n.true.\n$wdatas[$k]\n$scorrs[$k]\n$cons[$k]\n";
+ if ($k==2) {print INV ".true.\n";} else { print INV "${wns}\n";}
+ close(INV);
+ print BASH "# --- mech $k ----\n";
+ print BASH "cp -f INVERSION.$k INVERSION.PAR\n";
+ print BASH "cmt3d_flexwin > cmt3d.stdout\n";
+ print BASH "if [ \$? != 0 ]; then\n echo \"exiting case $k\";exit\nfi\n";
+ print BASH "perl -pi -e \"s/INV/$names[$k]/\" $cmt_new\n";
+ print BASH "mv -f $cmt_new $outdir/$cmt[$k]\n";
+ print BASH "mv -f cmt3d_flexwin.out $outdir/cmt3d_out.$k\n";
+ print BASH "mv -f cmt3d.stdout $outdir/cmt3d_stdout.$k\n";
+ print BASH "mv -f INVERSION.$k $outdir\n";
+}
+close(BASH);
+system("chmod a+x mech_table.bash; mech_table.bash");
+for ($k=2;$k<$nbeach;$k++) {
+ (@tmp)=split(" ",`grep Variance $outdir/cmt3d_stdout.$k`);
+ $var[$k]=$tmp[8]; print "$k -- VR = $var[$k]\n";}
+
+#die("here\n");
+# get depth information to display
+($odep)=split(" ",`cmtsol2faultpar.pl $cmt | grep mw | cut -d = -f 2 | cut -d / -f 2`);
+($ndep)=split(" ",`cmtsol2faultpar.pl $outdir/$cmt[2] | grep mw | cut -d = -f 2 | cut -d / -f 2`);
+($dep1)=split(" ",`cmtsol2faultpar.pl $outdir/$cmt[6] | grep mw | cut -d = -f 2 | cut -d / -f 2`);
+($dep2)=split(" ",`cmtsol2faultpar.pl $outdir/$cmt[7] | grep mw | cut -d = -f 2 | cut -d / -f 2`);
+($dep3)=split(" ",`cmtsol2faultpar.pl $outdir/$cmt[8] | grep mw | cut -d = -f 2 | cut -d / -f 2`);
+$ddep=sprintf("%5.2f",max(abs($dep1-$ndep),abs($dep2-$ndep),abs($dep3-$ndep)));
+
+open(BASH,">mech_plot.bash");
+print BASH "gmtset BASEMAP_TYPE plain ANOT_FONT_SIZE 9 HEADER_FONT_SIZE 10 MEASURE_UNIT inch PAPER_MEDIA letter TICK_LENGTH 0.1c\n";
+
+plot_psxy(\*BASH,$psfile,"$JX -K -X0 -Y0 $R","");
+$k=0;
+for ($i=0;$i<$nrows;$i++) {
+ for ($j=0;$j<$ncols;$j++) {
+ print BASH "# ---- Mech $k, $cmt[$k]------\n";
+ print "---- Mech $k, $cmt[$k]------\n";
+ $xy=$xy[$i][$j];($x,$y) = split(/\//,$xy);
+ plot_psxy(\*BASH,$psfile,"$JX -X$x -Y$y $R","");
+ if (-f "$outdir/$cmt[$k]") {
+ plot_psmeca_raw(\*BASH,$psfile,"-JX -R -Sm1.0 ","1.5 2","$outdir/$cmt[$k]");
+ @output = `cmtsol2faultpar.pl $outdir/$cmt[$k]`;
+ (undef,$tmp1) = split(" ", $output[5]);
+ (undef,$tmp2) = split(" ", $output[7]);
+ ($s,$d,$r,$p,$t) = split(/\//,$tmp2);
+ $tex="1 1 9 0 4 CM $names[$k]\n";
+ $tex.="0.2 0.8 9 0 4 LM Mw/dep/eps/= $tmp1\n0.2 0.6 9 0 4 LM S/D/R= $s/$d/$r";
+ if ($k>=2) {$tex.=sprintf("\n0.2 0.4 9 0 4 LM Var=%6.2f%",$var[$k]);}
+ plot_pstext(\*BASH,$psfile,"-JX -R -G0/0/255","$tex");
+ }
+ if ($k==1) {
+ ($nt) = split(" ",`grep 'T.sac.d' $meas[2] | wc | awk '{print \$1}'`);
+ ($nr) = split(" ",`grep 'R.sac.d' $meas[2] | wc | awk '{print \$1}'`);
+ ($nz) = split(" ",`grep 'Z.sac.d' $meas[2] | wc | awk '{print \$1}'`);
+ ($n2) = split(" ",`grep 'T002.*.sac.d' $meas[2] | wc | awk '{print \$1}'`);
+ ($n3) = split(" ",`grep 'T003.*.sac.d' $meas[2] | wc | awk '{print \$1}'`);
+ ($n6) = split(" ",`grep 'T006.*.sac.d' $meas[2] | wc | awk '{print \$1}'`);
+
+ plot_pstext(\*BASH,$psfile,"-JX -R -N","0.5 4 12 0 4 CM $ename ($nz/$nr/$nt; $n2/$n3/$n6; $odep/$ndep/$ddep)");}
+ plot_psxy(\*BASH,$psfile,"-JX -X-$x -Y-$y","");
+ $k++;
+ }
+}
+plot_psxy(\*BASH,$psfile,"-JX -O","");
+close(BASH);
+system("chmod a+x mech_plot.bash; mech_plot.bash");
Property changes on: seismo/3D/GRD_CMT3D/scripts/mech_table3.pl
___________________________________________________________________
Added: svn:executable
+ *
More information about the CIG-COMMITS
mailing list