[cig-commits] r14600 - in seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt: . scripts
carltape at geodynamics.org
carltape at geodynamics.org
Sun Apr 5 16:12:04 PDT 2009
Author: carltape
Date: 2009-04-05 16:12:03 -0700 (Sun, 05 Apr 2009)
New Revision: 14600
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models_three.pl
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_coverage.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models_basemap.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/ps2eps.pl
Log:
Updated GMT plotting scripts for tomography model cross sections. Added a version to allow subplots of model cross sections at different depths.
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_coverage.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_coverage.pl 2009-04-05 23:09:44 UTC (rev 14599)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_coverage.pl 2009-04-05 23:12:03 UTC (rev 14600)
@@ -120,7 +120,7 @@
$Wshelf = "-W0.75p,0/0/0,--";
# CMT sources
-$dir_sources = "/net/sierra/raid1/carltape/results/SOURCES/socal_16";
+$dir_sources = "/home/carltape/results/SOURCES/socal_16";
$outer_boundary = "${dir_sources}/SPECFEM_outer_points.dat";
$inner_boundary = "${dir_sources}/SPECFEM_inner_points.dat";
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl 2009-04-05 23:09:44 UTC (rev 14599)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl 2009-04-05 23:12:03 UTC (rev 14600)
@@ -15,7 +15,7 @@
# plot_horz_models.pl 1/50 m00/m16 0 1/1/1/0.3 # Vb
# plot_horz_models.pl 1/6 m00/m01 1 1/1/1/0.2 # Vs
#
-# SCIENCE PAPER:
+# FIGURES:
# plot_horz_models.pl 3/3 m00/m16 1 1/0/1/0.2 # Vs
# plot_horz_models.pl 11/11 m00/m16 1 0/0/1/0.2 # Vs
# plot_horz_models.pl 21/21 m00/m16 1 0/0/1/0.2 # Vs
@@ -23,6 +23,8 @@
# plot_horz_models.pl 11/11 m00/m16 0 0/0/1/0.3 # Vb
# plot_horz_models.pl 21/21 m00/m16 0 0/0/1/0.3 # Vb
#
+# plot_horz_models.pl 11/11 m00/m16 1 1/0/0/0.2 # Vs
+#
#==========================================================
if (@ARGV < 4) {die("Usage: plot_horz_models.pl xxx\n");}
@@ -40,7 +42,8 @@
#$itoplab = 0; # label the model above the plot
$iinsetlab = 1; # label the model inan inset
$iletter = 1; # label for publication figures
-$letter = "C";
+$letter = "c";
+$iextra = 0; # extra information on figures
$cgray = "-G220/220/220";
@@ -115,7 +118,7 @@
$Wshelf = "-W0.75p,0/0/0,--";
# CMT sources
-$dir_sources = "/net/sierra/raid1/carltape/results/SOURCES/socal_16";
+$dir_sources = "/home/carltape/results/SOURCES/socal_16";
$outer_boundary = "${dir_sources}/SPECFEM_outer_points.dat";
$inner_boundary = "${dir_sources}/SPECFEM_inner_points.dat";
@@ -223,7 +226,7 @@
for ($p = $pmin; $p <= $pmax; $p ++ ) {
$stip = sprintf("%3.3i",$p);
- $fname = "horz_xc_${modlab}_${smodeltag}_${stip}";
+ $fname = "horz_xc_${modlab}_${smodeltag}_${stip}_extra${iextra}";
$psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
$vsnorm = $vsnorms[$p-1]; # reference velocity for m00 and m16
@@ -353,26 +356,30 @@
#--------------------------------
- $finfo1 = "-W3p,0/0/0";
- $finfo2 = "-W1.5p,0/255/255";
+ if ($iextra==1) {
+ $finfo1 = "-W3p,0/0/0";
+ $finfo2 = "-W1.5p,0/255/255";
- # plot a line at -119 longitude
- if($p==11 && $ivs==1) {
- $xmark = -119;
- print CSH "psxy $J $R $finfo1 -K -O -V >>$psfile<<EOF\n$xmark $ymin\n$xmark $ymax\nEOF\n";
- print CSH "psxy $J $R $finfo2 -K -O -V >>$psfile<<EOF\n$xmark $ymin\n$xmark $ymax\nEOF\n";
- }
+ # plot a line at -119 longitude
+ if ($p==11 && $ivs==1) {
+ $xmark = -119;
+ print CSH "psxy $J $R $finfo1 -K -O -V >>$psfile<<EOF\n$xmark $ymin\n$xmark $ymax\nEOF\n";
+ print CSH "psxy $J $R $finfo2 -K -O -V >>$psfile<<EOF\n$xmark $ymin\n$xmark $ymax\nEOF\n";
+ }
- # plot two cross section ray paths
- if ($p==21 && $ivs==1) {
- @irs = (2,5); $nray = @irs ;
- for ($ik = 1; $ik <= $nray; $ik ++ ) {
- $ir = $irs[$ik-1];
- $stir = sprintf("%3.3i",$ir);
- $pfile = "/net/sierra/raid1/carltape/results/PLOTTING/SAVE/vert_01/vert_xc_${stir}_ray_path";
- if (not -f $pfile) {die("Check if pfile $pfile exist or not\n");}
- print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $J $R $finfo1 -K -O -V >> $psfile\n";
- print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $J $R $finfo2 -K -O -V >> $psfile\n";
+ # plot two cross section ray paths
+ if ($p==21 && $ivs==1) {
+ @irs = (2,5); $nray = @irs ;
+ for ($ik = 1; $ik <= $nray; $ik ++ ) {
+ $ir = $irs[$ik-1];
+ $stir = sprintf("%3.3i",$ir);
+ $pfile = "./INPUT/vert_01/vert_xc_${stir}_ray_path";
+ if (not -f $pfile) {
+ die("Check if pfile $pfile exist or not\n");
+ }
+ print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $J $R $finfo1 -K -O -V >> $psfile\n";
+ print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $J $R $finfo2 -K -O -V >> $psfile\n";
+ }
}
}
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models_three.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models_three.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models_three.pl 2009-04-05 23:12:03 UTC (rev 14600)
@@ -0,0 +1,416 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# plot_horz_models_three.pl
+# Carl Tape
+# 07-Dec-2008
+#
+# This script reads in a horizontal cross-section data file with five columns:
+# lon, lat, model-1, model-2, ln(model-2 / model-1)
+# and plots the three cross-sections.
+#
+# EXAMPLE:
+# plot_horz_models_three.pl 1/50 m00/m16 1 1/1/1/0.2 # Vs
+# plot_horz_models_three.pl 1/50 m00/m16 0 1/1/1/0.3 # Vb
+# plot_horz_models_three.pl 1/6 m00/m01 1 1/1/1/0.2 # Vs
+#
+# plot_horz_models_three.pl 3/11/21 m00/m16 1 1/0/0/0.2 # Vs
+#
+#==========================================================
+
+if (@ARGV < 4) {die("Usage: plot_horz_models_three.pl xxx\n");}
+($pinds,$smodels,$ivs,$parms) = @ARGV;
+
+$icolor = 1;
+$icuts = 1;
+
+($pmin,$pmid,$pmax) = split(/\//,$pinds);
+ at pinds = ($pmin,$pmid,$pmax);
+$nump = @pinds;
+$stips = sprintf("p%3.3i_p%3.3i_p%3.3i",$pmin,$pmid,$pmax);
+
+($smodel1,$smodel2) = split(/\//,$smodels);
+$smodeltag = "${smodel2}_${smodel1}";
+if($ivs==1) {$modlab = "vs"; $mtit = "Vs";} else {$modlab = "vb"; $mtit = "Vb";}
+($ilonlab,$itoplab,$isidelab,$ctick1b_km) = split(/\//,$parms);
+
+# detailed flags for publication figures
+#$ilonlab = 0; # display longitude tick labels at the top
+#$itoplab = 0; # label the model above the plot
+$iinsetlab = 1; # label the model inan inset
+$iletter = 0; # label for publication figures
+$letter = "C";
+
+$cgray = "-G220/220/220";
+
+#---------------------------------------------------------
+# USER INPUT
+
+# directory containing the data files
+$dirdat = "INPUT/horz_01";
+
+# KEY: file showing the cuts and the plotting range for the values
+$fcuts = "${dirdat}/horz_xc_elevation_cuts_mod";
+if (not -f $fcuts) {die("Check if fcuts $fcuts exist or not\n")}
+
+# load values and check that the files exist
+open(IN,$fcuts); @lines = <IN>; $nump0 = @lines;
+for ($p = 1; $p <= $nump0; $p ++ ) {
+ $stip = sprintf("%3.3i",$p);
+ (undef,$zcut,$vs_norm,$vb_norm,$vpert) = split(" ",$lines[$p-1]);
+ $vperts[$p-1] = $vpert;
+ $vsnorms[$p-1] = $vs_norm;
+ $vbnorms[$p-1] = $vb_norm;
+ $zcuts[$p-1] = -$zcut;
+ $dfiles[$p-1] = "${dirdat}/horz_xc_${modlab}_${smodeltag}_${stip}.dat";
+ #if (not -f $dfile) {die("Check if dfile $dfile exist or not\n")}
+}
+print "\nNumber of cuts is $nump0\n";
+if($pmax > $nump0) {die("pmax > nump0")}
+if($pmin < 1) {die("pmin < 1")}
+
+# subtitles
+ at titles = ("$mtit model $smodel1","$mtit model $smodel2","ln ( ${smodel2} / ${smodel1} )");
+#@labs = ("$smodel1","$smodel2","ln(${smodel2}/${smodel1})");
+
+# parameters for color scales
+$cpert1 = 0.2;
+#$cpert2 = 0.2;
+$scale_color = 65.0;
+$colorbar = "seis";
+
+# resolution of color plots
+#$interp = "-I2m/2m -S4m"; # key information
+$interp = "-I2m";
+$grdfile = "temp.grd";
+
+# position of titles and labels
+$xtx1 = 0.5; $ytx1 = 1.00;
+$xtx2 = -0.15; $ytx2 = 0.4;
+$xtx3 = 0.7; $ytx3 = -0.1;
+$xtx4 = $xtx3; $ytx4 = -0.3;
+$xtx5 = -0.15; $ytx5 = 0.85; # A, B, C, etc
+$xtx6 = 0.0; $ytx6 = 0.05; # m00, m16, ln(m16/m00)
+
+#print "$dtitle\n"; die("TESTING");
+
+#---------------------------------------------------------
+# DATA FILES
+
+$dir0 = "/home/carltape/gmt";
+$dir = "socal_2005";
+
+$fault_file = "$dir0/faults/jennings_more.xy";
+$fault_file2 = "$dir0/faults/scitex.lin";
+#$kcf_file = "$dir0/faults/kcf.xy";
+#$breck_file = "$dir0/faults/breck.xy";
+
+$plate_file = "$dir0/plates/plate_boundaries/bird_boundaries";
+$plate_labels = "$dir0/socal_2005/socal_plate_labs.xyz";
+
+# continental shelf
+$ishelf = 0;
+$shelf_file = "/home/carltape/gmt/coastlines/socal_shelf";
+$Wshelf = "-W0.75p,0/0/0,--";
+
+# CMT sources
+$dir_sources = "/home/carltape/results/SOURCES/socal_16";
+$outer_boundary = "${dir_sources}/SPECFEM_outer_points.dat";
+$inner_boundary = "${dir_sources}/SPECFEM_inner_points.dat";
+
+# cities
+$cities = "$dir0/cities/socal_tomo_cities";
+
+# check that files exist
+#if (not -f $xxx) { die("Check if $xxx exist or not\n") }
+
+# bounds and dimensions
+# 0 = LA
+# 1 = southern California
+$iregion = 1;
+
+if ($iregion==1) {
+ # southern California
+ $wid = 2.8; # width of figure (inches)
+ $xmin = -122; $xmax = -114; $ymin = 31.5; $ymax = 37;
+ $xmin = -121.2; $xmax = -114.8; $ymin = 32.3; $ymax = 36.7; # SPECFEM
+ $tick1 = 1; $tick2 = 0.5;
+ $cmax = 5000; $cmin = -$cmax;
+ $y_title = 0.98; $y_title2 = 0.92; $x_title = 0.5; $x_title2 = 0.5;
+ $iportrait = 0;
+ $stag = "Southern California";
+ $topo_labels = "$dir0/socal_2005/socal_topo_labs.xyz";
+ $fault_labels = "$dir0/socal_2005/socal_fault_labs.xyz";
+ $name = "plot_horz_models_three";
+ $origin = "-X1i -Y5i";
+
+ # inset map
+ $iinset = 1;
+ $origin_inset = "-Xa7.5 -Ya4.5";
+ $Jinset = "-JM1.5";
+ $Rinset = "-R-132/-110/25/50";
+ $Binset = "-B100wesn";
+ $coast_res = "-Df -A0/0/4";
+}
+
+$dX = $wid + 0.5;
+$hwid = $wid/2;
+#$R = "-Rg";
+$R = "-R$xmin/$xmax/$ymin/$ymax";
+#$J = "-JK180/${wid}i";
+$J = "-JM${wid}";
+#$J = "-Jm0.70";
+
+# scale for plots
+ at Bopts = ("WESN","Wesn","wesN","wEsn","weSn","WesN","wEsN","wESn","WeSn","WEsn","weSN","wESN","WESn","WeSN","WEsN","wesn");
+$B0 = "-Ba${tick1}f${tick2}d:.\" \":";
+#$B0 = "-Ba${tick1}f${tick1}g${tick1}:.\" \":"; # with gridlines
+
+$B = "$B0".$Bopts[0];
+
+if($iportrait == 1){$orient = "-P"; $rotangle = 0}
+else {$orient = " "; $rotangle = 90}
+
+# color bar
+$Dlen = $wid*0.6;
+$Dx = $Dlen*0.5;
+$Dy = -0.5;
+$Dscale = "-D$Dx/$Dy/$Dlen/0.15h";
+
+# plotting specifications
+$fsize0 = "14";
+$fsize1 = "12";
+$fsize2 = "10";
+$fsize3 = "6";
+$fontno = "1"; # 1 or 4
+$tick = "6p";
+$fpen = "1p";
+$tpen = "1p";
+
+# plotting specificiations
+# NOTE: pen attribute (-W) does not work for psvelo
+#$coast_info = "-A5000 -Dl -W1.0p -S220/255/255 -C220/255/255 -G200"; # A : smallest feature plotted, in km^2; D : resolution
+
+$coast_info = "${coast_res} -W0.75p -Na/0.75p";
+$fault_info = "-M -W0.75p,255/0/0";
+$fault_infoK = "-M -W0.75p,0/0/0";
+
+$city_info = "-Sc10p -W1.0p/0/0/0 -G255/255/0";
+
+$textinfo = "-G255 -S1.5p";
+$textinfo2 = "-G0/0/0 -S2p,255/255/255";
+$textinfo3 = "-G0/0/0 -S2p,255/255/0";
+
+# BOOLEAN: WHICH FIGURES TO PLOT
+$ixv = 0; $ipdf = 1;
+#$ixv = 1; $ipdf = 0;
+
+# plot title
+$J_title = "-JM${wid}";
+$R_title = "-R0/1/0/1";
+$fsize_title = $fsize0;
+
+#==================================================
+
+$cshfile = "plot_horz_models_three.csh";
+open(CSH,">$cshfile");
+print CSH "gmtset PAPER_MEDIA letter BASEMAP_TYPE plain PLOT_DEGREE_FORMAT D TICK_LENGTH $tick LABEL_FONT_SIZE $fsize2 ANOT_FONT_SIZE $fsize2 HEADER_FONT $fontno ANOT_FONT $fontno LABEL_FONT $fontno HEADER_FONT_SIZE $fsize1 FRAME_PEN $fpen TICK_PEN $tpen \n";
+
+ at shifts = ($origin,"-X$dX","-X$dX");
+if($ilonlab==1) {@iBs = (2,14,2);} else {@iBs = (15,9,15);}
+
+for ($w = 1; $w <= $nump; $w ++ ) {
+
+ $p = $pinds[$w-1];
+
+ $stip = sprintf("%3.3i",$p);
+ #$fname = "horz_xc_${modlab}_${smodeltag}_${stip}";
+ $fname = "horz_three_${modlab}_${smodel2}_${stips}_cuts${icuts}";
+ $psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
+
+ $vsnorm = $vsnorms[$p-1]; # reference velocity for m00 and m16
+ $vbnorm = $vbnorms[$p-1]; # reference velocity for m00 and m16
+ $vpert = $vperts[$p-1]; # perturbation from reference velocity to plot
+ $zdep = $zcuts[$p-1]; # depth of the cross-section (m)
+
+ # make less of a range for Vb compared with Vs -- factor
+ if($ivs == 0) {$vpert = 1.0*$vpert;}
+
+ $dtitle = sprintf("Depth = %.1f km",$zdep/1000);
+
+ # datafile -- COLUMNS: lon, lat, m00, m16
+ $dfile = $dfiles[$p-1];
+ if (not -f $dfile) {die("Check if dfile $dfile exist or not\n");}
+ print "Data file is $dfile\n";
+
+ if($ivs == 1) {$cnorm = $vsnorm;} else {$cnorm = $vbnorm;}
+ $cpert1 = $vpert;
+ $cpert2 = $vpert; # NOTE: over-rule the input value
+
+ # colorpoint file for m00 and m16 -- VARIES with depth -- in ln(m00/cnorm)
+ $cptfile1a = "color1a.cpt";
+ $cmin = -$cpert1*1.01; $cmax = $cpert1*1.01;
+ $dc = ($cmax-$cmin)/${scale_color};
+ $T = sprintf("-T%3.3e/%3.3e/%3.3e",$cmin,$cmax,$dc);
+ #print CSH "makecpt -C$colorbar $T -D > $cptfile1a\n";
+ print CSH "makecpt -C$colorbar $T -D > $cptfile1a\n";
+
+ # colorpoint file for m00 and m16 -- VARIES with depth -- in KILOMETERS
+ $cptfile1b = "color1b.cpt";
+ $cmin = $cnorm*exp(-$cpert1)*1e-3;
+ $cmax = $cnorm*exp($cpert1)*1e-3;
+ $dc = ($cmax-$cmin)/${scale_color};
+ $T = sprintf("-T%3.3e/%3.3e/%3.3e",$cmin,$cmax,$dc);
+ print CSH "makecpt -C$colorbar $T -D > $cptfile1b\n";
+
+ # colorpoint file for perturbation ln(m16/m00) -- FIXED with depth
+ $cptfile2 = "color2.cpt";
+ $cmin = -$cpert2*1.01; $cmax = $cpert2*1.01;
+ $dc = ($cmax-$cmin)/${scale_color};
+ $T = sprintf("-T%3.3e/%3.3e/%3.3e",$cmin,$cmax,$dc);
+ print CSH "makecpt -C$colorbar $T -D > $cptfile2\n";
+
+ #
+ $Bscale1a = sprintf("-Ba%2.2ef0.05:\" \": -E10p -A",$cpert1);
+ $Bscale1b = sprintf("-B%2.2ef0.1:\" \": -E10p",$ctick1b_km);
+ $Bscale2 = sprintf("-Ba%2.2ef0.05:\" \": -E10p -A",$cpert2);
+
+# attempting to get shaded relief superimposed
+#$gradfile0 = "/home/datalib/Topography/GLOBAL/ETOPO1/ETOPO1_Ice_g.grad";
+#if (not -f $gradfile0) {die("Check if gradfile $gradfile0 exist or not\n");}
+
+
+ $shift = $shifts[$w-1];
+ $iB = $iBs[$w-1];
+ $B = "$B0".$Bopts[$iB];
+ $title = $titles[$w-1];
+
+ if ($w==1) {
+ print CSH "psbasemap $J $R $B -K -V $orient $origin > $psfile\n"; # START
+ } else {
+ print CSH "psbasemap $J $R $B -K -O -V $shift >> $psfile\n";
+ }
+ print CSH "psbasemap $J $R $B $cgray -K -O -V >> $psfile\n";
+ if ($icolor==1) {
+ #print CSH "awk '{print \$1,\$2,log(\$4/$cnorm)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
+ print CSH "awk '{print \$1,\$2,log(\$4/$cnorm)}' $dfile | xyz2grd -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile1a $J -Q -K -O -V >> $psfile\n";
+
+ print CSH "psscale -C${cptfile1a} $Dscale $Bscale1a -K -O -V >> $psfile\n";
+ print CSH "psscale -C${cptfile1b} $Dscale $Bscale1b -K -O -V >> $psfile\n";
+ $slab1 = sprintf("ln(%s / %.2f)",$mtit,$cnorm/1000); $slab2 = "$mtit km/s";
+ print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx3 $ytx3 $fsize2 0 $fontno LM $slab1\nEOF\n";
+ print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx4 $ytx4 $fsize2 0 $fontno LM $slab2\nEOF\n";
+ }
+ print CSH "pscoast $J $R $coast_info -K -O -V >> $psfile \n";
+ if($ishelf==1) {print CSH "awk '{print \$2,\$1}' ${shelf_file} |psxy $J $R $Wshelf -K -O -P -V >> $psfile\n";}
+
+ #--------------------------------
+
+ #print CSH "psxy ${plate_file} $J $R $plate_info -K -V -O >> $psfile \n";
+ print CSH "psxy ${fault_file} $J $R $fault_infoK -K -V -O >> $psfile \n";
+ #print CSH "psxy ${kcf_file} $J $R $fault_infoK -K -V -O >> $psfile \n";
+ #print CSH "psxy ${breck_file} $J $R $fault_infoK -K -V -O >> $psfile \n";
+ #print CSH "psxy ${fault_file2} $J $R $fault_infoK -K -V -O >> $psfile \n";
+
+ if (0==1) {
+ # boundaries of simulation
+ print CSH "psxy ${outer_boundary} $J $R -W2p,0/0/0 -K -O -V >>$psfile\n";
+ print CSH "psxy ${inner_boundary} $J $R -W1.5p,0/0/0,-- -K -O -V >>$psfile\n";
+
+ $ibox = 0;
+ if ($ibox==1) {
+ $boxinfo = "-W1.5p,0/0/255"; # -A : suppress drawing line segments as great circle arcs
+
+ # Lin model (2007) box
+ $lin_boundary = "/net/denali/home2/carltape/gmt/tomography/lin_2007/lin_boundary_points.dat";
+ print CSH "psxy ${lin_boundary} $J $R $boxinfo -K -O -V >>$psfile\n";
+ }
+ }
+
+ #--------------------------------
+
+ if ($icuts==1) {
+ $finfo1 = "-W3p,0/0/0";
+ $finfo2 = "-W1.5p,0/255/255";
+
+ # plot a line at -119 longitude
+ if ($p==11 && $ivs==1) {
+ $xmark = -119;
+ print CSH "psxy $J $R $finfo1 -K -O -V >>$psfile<<EOF\n$xmark $ymin\n$xmark $ymax\nEOF\n";
+ print CSH "psxy $J $R $finfo2 -K -O -V >>$psfile<<EOF\n$xmark $ymin\n$xmark $ymax\nEOF\n";
+ }
+
+ # plot two cross section ray paths
+ if ($p==3 && $ivs==1) {
+ #if ($p==21 && $ivs==1) {
+ @irs = (4,10,25); $nray = @irs;
+ #@irs = (2,5); $nray = @irs;
+ for ($ik = 1; $ik <= $nray; $ik ++ ) {
+ $ir = $irs[$ik-1];
+ $stir = sprintf("%3.3i",$ir);
+ $pfile = "./INPUT/vert_01/vert_xc_${stir}_ray_path";
+ if (not -f $pfile) {die("Check if pfile $pfile exist or not\n");}
+ print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $J $R $finfo1 -K -O -V >> $psfile\n";
+ print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $J $R $finfo2 -K -O -V >> $psfile\n";
+ }
+ }
+ }
+
+ #--------------------------------
+
+ print CSH "psbasemap $J $R $B -K -V -O >> $psfile\n";
+
+ # plot cities
+ ##print CSH "psxy $J $R $city_info -K -O -V >>$psfile<<EOF\n -118.1717 34.1483 \nEOF\n"; # Pasadena
+ #print CSH "awk '{print \$2,\$1}' $cities | psxy ${city_info} $J $R -K -V -O >> $psfile \n";
+ #print CSH "awk '{print \$2,\$1,12,0,1,\"LB\",\$3}' $cities | pstext $textinfo3 -D0.05/0.05 $J $R -K -V -O >> $psfile \n";
+
+ #print CSH "pstext ${topo_labels} $J $R $textinfo -N -K -V -O >> $psfile\n";
+ #print CSH "pstext ${fault_labels} $textinfo2 $J $R -K -V -O >> $psfile \n";
+
+ #----------------------------
+
+ # plot horizontal title above each column
+ if ($itoplab==1) {
+ print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx1 $ytx1 $fsize0 0 $fontno CM $title\nEOF\n";
+ }
+
+ # inset label for each plot
+ if ($iinsetlab==1) {
+ $lab = $dtitle;
+ $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V >>$psfile<<EOF\n$xtx6 $ytx6 11 0 $fontno LB $lab\nEOF\n";
+ }
+
+ # plot vertical title left of the row
+ #if ($isidelab==1 && $k==1) {
+ # print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx2 $ytx2 $fsize0 90 $fontno CM $dtitle\nEOF\n";
+ #}
+
+ # plot overall label (for publication)
+ #if ($iletter==1 && $k==1) {
+ # $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ # print CSH "pstext $R_title $J_title $textinfo -K -O -V >>$psfile<<EOF\n$xtx5 $ytx5 18 0 $fontno TL $letter\nEOF\n";
+ #}
+
+} # loop over w
+
+ print CSH "pstext -N $R_title $J_title -O -V >>$psfile<<EOF\n $x_title $y_title 16 0 $fontno CM \nEOF\n"; # FINISH
+ #if($ixv==1) {print CSH "convert $psfile -rotate $rotangle $jpgfile\n";}
+ if($ixv==1) {print CSH "ghostview $psfile &\n";}
+ if($ipdf==1) {print CSH "ps2pdf $psfile\n";}
+
+#------------------------------------
+# print CSH "pstext -N $R_title $J_title -O -V >>$psfile<<EOF\n $x_title $y_title 16 0 $fontno CM \nEOF\n"; # FINISH
+
+close (CSH);
+system("csh -f $cshfile");
+
+# print "convert $psfile $jpgfile \n";
+# system("convert $psfile -rotate $rotangle $jpgfile");
+# if ($ipdf==1) {system("ps2pdf $psfile")}
+# if ($ixv==1) {system("ghostview $psfile &");}
+
+#==================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models_three.pl
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models.pl 2009-04-05 23:09:44 UTC (rev 14599)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models.pl 2009-04-05 23:12:03 UTC (rev 14600)
@@ -29,9 +29,9 @@
$icolor = 1;
$ivertlab = 0; # vertical exaggeration label
$isimbox = 1; # box around simulation region
-$iletter = 0; # label for publication figures
-$letter = "A";
-$imfinal = 0; # plot a figure with only m16 and the map
+$iletter = 1; # label for publication figures
+$letter = "b";
+$imfinal = 1; # plot a figure with only m16 and the map
($pmin,$pmax) = split(/\//,$pinds);
($smodel1,$smodel2) = split(/\//,$smodels);
@@ -81,7 +81,7 @@
$Wshelf = "-W0.75p,0/0/0,--";
# CMT sources
-$dir_sources = "/net/sierra/raid1/carltape/results/SOURCES/socal_16";
+$dir_sources = "/home/carltape/results/SOURCES/socal_16";
$outer_boundary = "${dir_sources}/SPECFEM_outer_points.dat";
$inner_boundary = "${dir_sources}/SPECFEM_inner_points.dat";
@@ -145,6 +145,7 @@
# BOOLEAN: WHICH FIGURES TO PLOT
#$ixv = 1; $ipdf = 0;
$ixv = 0; $ipdf = 1;
+$ijpg = 0;
# plot title
$J_title = "-JM1";
@@ -493,12 +494,17 @@
if($irun==3) {print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$gf $zmax0\n$gf $zmin0\nEOF\n";}
}
- #if ($p==6) {
- # $xlarse = $saf+60;
- # print CSH "psxy $J $R -W2p,0/0/0 -K -O -V $oxc >>$psfile<<EOF\n$saf -19\n$xlarse -21\nEOF\n"; # LARSE reflector
- #}
- #}
+ if ($p==6 && 1==1) {
+ $xlarse = $saf+60;
+ $zminf0 = -5;
+ $finfo = "-W1.5p,0/0/0,--";
+ print CSH "psxy $J $R -W2p,0/0/0 -K -O -V $oxc >>$psfile<<EOF\n$saf -19\n$xlarse -21\nEOF\n"; # LARSE reflector
+ print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$saf $zmax0\n$saf $zmin0\nEOF\n";
+ print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$gf $zmax0\n$gf $zmin0\nEOF\n";
+ print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$sgf $zmax0\n$sgf $zmin0\nEOF\n";
+ }
+
# plot source and receiver
if ($irun==1) {
print CSH "awk '{print \$3,-\$4}' $rfile | psxy $J $R $srcinfo_xc -K -O -V $oxc >> $psfile\n";
@@ -534,16 +540,16 @@
} # loop over k
# plot overall label (for publication)
- if ($iletter==1) {
- $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
- print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 18 0 $fontno TL $letter\nEOF\n";
- } else {
- if($irun==1) {
+ if ($iletter==1) {
+ $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 18 0 $fontno TL $letter\nEOF\n";
+ } else {
+ if ($irun==1) {
$textinfo = "-N";
print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 10 0 $fontno TC $eid\nEOF\n";
print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab2 >>$psfile<<EOF\n0 0 10 0 $fontno TC $stanet\nEOF\n";
}
- }
+ }
# plot label on each cross section
for ($k = 1; $k <= 3; $k ++ ) {
@@ -566,8 +572,7 @@
#if($ixv==1) {print CSH "convert $psfile -rotate $rotangle $jpgfile\n";}
if ($ixv==1) {print CSH "ghostview $psfile &\n";}
if ($ipdf==1) {print CSH "ps2pdf $psfile\n";}
-
-} # loop over p
+ if ($ijpg==1) {print CSH "convert $psfile $jpgfile\n";}
#==================================================
# figure with map with m16 beneath
@@ -578,11 +583,11 @@
@iBs = (5,8,4);
- for ($p = $pmin; $p <= $pmax; $p ++ ) {
+ #for ($p = $pmin; $p <= $pmax; $p ++ ) {
$stip = sprintf("%3.3i",$p);
+ #$fname = "vert_${stirun}_xc_${modlab}_${smodeltag}_${stip}_single";
$fname = "${fname}_single";
-
$psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
# datafile -- COLUMNS: lon, lat, m00, m16
@@ -621,6 +626,9 @@
$x0_label = $x2;
$y0_label = $y0_xc + $heightxc;
$olab = "-Xa${x0_label} -Ya${y0_label}";
+ $x0_label2 = $x2;
+ $y0_label2 = $y0_xc + 0.8*$heightxc;
+ $olab2 = "-Xa${x0_label2} -Ya${y0_label2}";
# map
$xgap_lab_map = 0.4;
@@ -653,15 +661,17 @@
print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $Jmap $Rmap -W3p,0/255/255 -K -O -V $omap >> $psfile\n";
# plot source and receiver lon-lat
- $rfile = "${dirdat}/vert_xc_${stip}_A_B";
- if (not -f $rfile) {die("Check if rfile $rfile exist or not\n");}
- print CSH "awk '{print \$1,\$2}' ${rfile} | psxy $Jmap $Rmap $srcinfo_map -K -O -V $omap >> $psfile\n";
- print CSH "awk '{print \$5,\$6}' ${rfile} | psxy $Jmap $Rmap $recinfo_map -K -O -V $omap >> $psfile\n";
+ if ($irun==1) {
+ $rfile = "${dirdat}/vert_xc_${stip}_A_B";
+ if (not -f $rfile) {die("Check if rfile $rfile exist or not\n");}
+ print CSH "awk '{print \$1,\$2}' ${rfile} | psxy $Jmap $Rmap $srcinfo_map -K -O -V $omap >> $psfile\n";
+ print CSH "awk '{print \$5,\$6}' ${rfile} | psxy $Jmap $Rmap $recinfo_map -K -O -V $omap >> $psfile\n";
+ }
print CSH "psbasemap $Jmap $Rmap $Bmap -K -O -V $omap >> $psfile\n";
# plot additional ray paths
- if ($p==10) {
+ if (($irun==1) && ($p==10) ) {
$x0_mapinset = $x0_map + 2;
$y0_mapinset = $y0_map + 0;
$omapinset = "-Xa${x0_mapinset} -Ya${y0_mapinset}";
@@ -743,8 +753,10 @@
#}
# plot source and receiver
- print CSH "awk '{print \$3,-\$4}' $rfile | psxy $J $R $srcinfo_xc -K -O -V $oxc >> $psfile\n";
- print CSH "awk '{print \$7,$zmax0}' $rfile | psxy $J $R $recinfo_xc -K -O -V $oxc >> $psfile\n";
+ if ($irun == 1) {
+ print CSH "awk '{print \$3,-\$4}' $rfile | psxy $J $R $srcinfo_xc -K -O -V $oxc >> $psfile\n";
+ print CSH "awk '{print \$7,$zmax0}' $rfile | psxy $J $R $recinfo_xc -K -O -V $oxc >> $psfile\n";
+ }
# details for specific cross-sections
#if ($irun==1) {
@@ -770,8 +782,11 @@
$textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 18 0 $fontno TL $letter\nEOF\n";
} else {
- $textinfo = "-N";
- print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 12 0 $fontno TR $eid to $stanet\nEOF\n";
+ if ($irun==1) {
+ $textinfo = "-N";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 10 0 $fontno TC $eid\nEOF\n";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab2 >>$psfile<<EOF\n0 0 10 0 $fontno TC $stanet\nEOF\n";
+ }
}
# plot label on each cross section
@@ -788,17 +803,13 @@
print CSH "pstext -N $R_title $J_title -O -V >>$psfile<<EOF\n 10 10 16 0 $fontno CM \nEOF\n"; # FINISH
#if($ixv==1) {print CSH "convert $psfile -rotate $rotangle $jpgfile\n";}
- if ($ixv==1) {
- print CSH "ghostview $psfile &\n";
- }
- if ($ipdf==1) {
- print CSH "ps2pdf $psfile\n";
- }
+ if ($ixv==1) {print CSH "ghostview $psfile &\n";}
+ if ($ipdf==1) {print CSH "ps2pdf $psfile\n";}
+ if ($ijpg==1) {print CSH "convert $psfile $jpgfile\n";}
- } # loop over p
+ } # imfinal == 1
+} # loop over p
-} # imfinal == 1
-
#==================================================
#------------------------------------
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models_basemap.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models_basemap.pl 2009-04-05 23:09:44 UTC (rev 14599)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models_basemap.pl 2009-04-05 23:12:03 UTC (rev 14600)
@@ -4,15 +4,21 @@
#
# plot_vert_models_basemap.pl
# Carl Tape
-# 07-Dec-2008
+# 02-April-2009
#
# This script reads in a horizontal cross-section data file with five columns:
# lon, lat, model-1, model-2, ln(model-2 / model-1)
# and plots the three cross-sections.
#
# EXAMPLE:
-# plot_vert_models_basemap.pl 1
+# plot_vert_models_basemap.pl 1 2 3
+# plot_vert_models_basemap.pl 4 5 6 7 8 9 10 12 13 15 17 18 19 20 21 23 24 26 27 28 29 30 31 32 # SAF
+# plot_vert_models_basemap.pl 1 2 3 4 5 6 7 8 9 10 # Garlock
+# plot_vert_models_basemap.pl 1 2 3 4 5 6 7 8 9 10 11 12 13 # lon lines
+# plot_vert_models_basemap.pl 1 2 3 4 5 6 7 8 9 # lat lines
#
+# plot_vert_models_basemap.pl 4 5 10 25
+#
#==========================================================
if (@ARGV < 1) {die("Usage: plot_vert_models_basemap.pl xxx\n");}
@@ -20,7 +26,7 @@
$nump = @pinds;
$stp = sprintf("%3.3i",$nump);
-$irun = 1;
+$irun = 1; # 1, 2 (SAF), 3 (Garlock)
$icolor = 1;
$ivertlab = 0; # vertical exaggeration label
$isimbox = 1; # box around simulation region
@@ -38,7 +44,7 @@
# directory containing the data files
$stirun = sprintf("%2.2i",$irun);
-$dirdat = "/net/sierra/raid1/carltape/results/PLOTTING/SAVE/vert_${stirun}";
+$dirdat = "INPUT/vert_${stirun}";
#---------------------------------------------------------
@@ -76,7 +82,7 @@
$Wshelf = "-W0.75p,0/0/0,--";
# CMT sources
-$dir_sources = "/net/sierra/raid1/carltape/results/SOURCES/socal_16";
+$dir_sources = "/home/carltape/results/SOURCES/socal_16";
$outer_boundary = "${dir_sources}/SPECFEM_outer_points.dat";
$inner_boundary = "${dir_sources}/SPECFEM_inner_points.dat";
@@ -235,11 +241,13 @@
die("Check if rfile $rfile exist or not\n");
}
+ if ($irun==1) {
# plot source
print CSH "awk '{print \$1,\$2}' ${rfile} | psxy $Jmap $Rmap $srcinfo_map -K -O -V $omap >> $psfile\n";
# plot receiver
print CSH "awk '{print \$5,\$6}' ${rfile} | psxy $Jmap $Rmap $recinfo_map -K -O -V $omap >> $psfile\n";
+}
# plot cities
##print CSH "psxy $Jmap $Rmap $city_info -K -O -V >>$psfile<<EOF\n -118.1717 34.1483 \nEOF\n"; # Pasadena
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/ps2eps.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/ps2eps.pl 2009-04-05 23:09:44 UTC (rev 14599)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/ps2eps.pl 2009-04-05 23:12:03 UTC (rev 14600)
@@ -1,16 +1,22 @@
#!/usr/bin/perl -w
-$irun = 3;
+$ivert = 0;
+$irun = 1;
$stip = sprintf("%2.2i",$irun);
- at epsfiles = glob("vert_${stip}_xc*.ps");
-$neps = @epsfiles;
-for ($k = 1; $k <= $neps; $k = $k+1) {
- $efile = $epsfiles[$k-1];
- `ps2eps -f $efile`;
-}
+#$ftag = "vert_${stip}_xc";
+#$ftag = "horz_xc_vb_m16_m00";
+$ftag = "horz_xc_vs_m01_m00";
-if ($irun==1) {
+ at epsfiles = glob("${ftag}*.ps");
+
+ $neps = @epsfiles;
+ for ($k = 1; $k <= $neps; $k = $k+1) {
+ $efile = $epsfiles[$k-1];
+ `ps2eps -f $efile`;
+ }
+
+if ($irun==1 && $ivert==1) {
`ps2ps vert_01_xc_HEC_CI_282_9968977_vs_m16_m00_003.ps temp.ps ; ps2eps temp.ps ; mv temp.eps vert_01_xc_HEC_CI_282_9968977_vs_m16_m00_003.eps`;
`ps2ps vert_01_xc_LCG_CI_122_10215753_vs_m16_m00_021.ps temp.ps ; ps2eps temp.ps ; mv temp.eps vert_01_xc_LCG_CI_122_10215753_vs_m16_m00_021.eps`;
`ps2ps vert_01_xc_LGU_CI_353_10097009_vs_m16_m00_067.ps temp.ps ; ps2eps temp.ps ; mv temp.eps vert_01_xc_LGU_CI_353_10097009_vs_m16_m00_067.eps`;
@@ -22,4 +28,10 @@
`ps2ps vert_01_xc_DPP_CI_314_9753485_vs_m16_m00_079.ps temp.ps ; ps2eps temp.ps ; mv temp.eps vert_01_xc_DPP_CI_314_9753485_vs_m16_m00_079.eps`;
}
+if ($irun==1 && $ivert==0) {
+ `ps2ps ${ftag}_001.ps temp.ps ; ps2eps temp.ps ; mv temp.eps ${ftag}_001.eps`;
+ `ps2ps ${ftag}_009.ps temp.ps ; ps2eps temp.ps ; mv temp.eps ${ftag}_009.eps`;
+ `ps2ps ${ftag}_016.ps temp.ps ; ps2eps temp.ps ; mv temp.eps ${ftag}_016.eps`;
+}
+
#==================================================
More information about the CIG-COMMITS
mailing list