[cig-commits] r15033 - in seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt: . scripts
carltape at geodynamics.org
carltape at geodynamics.org
Fri May 22 06:36:23 PDT 2009
Author: carltape
Date: 2009-05-22 06:36:23 -0700 (Fri, 22 May 2009)
New Revision: 15033
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_one.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_three.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_vert_models_basemap.pl
Removed:
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models_three.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models_basemap.pl
Modified:
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
Log:
Updated perl+GMT plotting scripts for tomographic cross sections.
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-05-22 13:18:45 UTC (rev 15032)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl 2009-05-22 13:36:23 UTC (rev 15033)
@@ -13,7 +13,6 @@
# EXAMPLE:
# plot_horz_models.pl 1/50 m00/m16 1 1/1/1/0.2 # Vs
# 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
#
# FIGURES:
# plot_horz_models.pl 3/3 m00/m16 1 1/0/1/0.2 # Vs
@@ -23,7 +22,7 @@
# 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
+# plot_horz_models.pl 3/3 m00/m01 1 1/0/1/0.2 # Vs
#
#==========================================================
@@ -41,7 +40,7 @@
#$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 = 1; # label for publication figures
+$iletter = 0; # label for publication figures
$letter = "c";
$iextra = 0; # extra information on figures
@@ -91,8 +90,8 @@
# 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;
+$xtx3 = 0.7; $ytx3 = -0.2; # cbar 1
+$xtx4 = $xtx3; $ytx4 = -0.1; # cbar 2
$xtx5 = -0.15; $ytx5 = 0.85; # A, B, C, etc
$xtx6 = 0.0; $ytx6 = 0.05; # m00, m16, ln(m16/m00)
@@ -176,16 +175,20 @@
else {$orient = " "; $rotangle = 90}
# color bar
+$Dwid = 0.15;
$Dlen = $wid*0.6;
$Dx = $Dlen*0.5;
-$Dy = -0.5;
-$Dscale = "-D$Dx/$Dy/$Dlen/0.15h";
+$Dy = -0.2;
+$Dy2 = $Dy-$Dwid;
+$Dy2 = $Dy;
+$Dscale = "-D$Dx/$Dy/$Dlen/${Dwid}h";
+$Dscale2 = "-D$Dx/$Dy2/$Dlen/${Dwid}h";
# plotting specifications
$fsize0 = "14";
$fsize1 = "12";
$fsize2 = "10";
-$fsize3 = "6";
+$fsize3 = "8";
$fontno = "1"; # 1 or 4
$tick = "6p";
$fpen = "1p";
@@ -253,16 +256,17 @@
$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
+ # THIS IS FOR THE COLORBAR ONLY
$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";
+ $Bscale1b = sprintf("-B%2.2ef0.1:\" \": -E10p",$ctick1b_km);
# colorpoint file for perturbation ln(m16/m00) -- FIXED with depth
$cptfile2 = "color2.cpt";
@@ -271,10 +275,9 @@
$T = sprintf("-T%3.3e/%3.3e/%3.3e",$cmin,$cmax,$dc);
print CSH "makecpt -C$colorbar $T -D > $cptfile2\n";
- #
+ # ticks for colorbars
$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);
+ $Bscale2 = sprintf("-Ba%2.2ef0.05:\" \": -E10p",$cpert2);
# attempting to get shaded relief superimposed
#$gradfile0 = "/home/datalib/Topography/GLOBAL/ETOPO1/ETOPO1_Ice_g.grad";
@@ -296,24 +299,32 @@
print CSH "psbasemap $J $R $B $cgray -K -O -V >> $psfile\n";
if ($icolor==1) {
if ($k == 1) {
- #print CSH "awk '{print \$1,\$2,log(\$3/$cnorm)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
- print CSH "awk '{print \$1,\$2,log(\$3/$cnorm)}' $dfile | xyz2grd -G$grdfile $R $interp\n";
- print CSH "grdimage $grdfile -C$cptfile1a $J -Q -K -O -V >> $psfile\n";
+ ##print CSH "awk '{print \$1,\$2,log(\$3/$cnorm)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
+ #print CSH "awk '{print \$1,\$2,log(\$3/$cnorm)}' $dfile | xyz2grd -G$grdfile $R $interp\n";
+ #print CSH "grdimage $grdfile -C$cptfile1a $J -Q -K -O -V >> $psfile\n";
+ print CSH "awk '{print \$1,\$2,\$3/1000}' $dfile | xyz2grd -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile1b $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 "psscale -C${cptfile1a} $Dscale $Bscale1a -K -O -V >> $psfile\n";
+ print CSH "psscale -C${cptfile1b} $Dscale2 $Bscale1b -K -O -V >> $psfile\n";
+ #$slab1 = sprintf("ln(%s / %.2f)",$mtit,$cnorm/1000);
+ $slab1 = sprintf("(%.2f \261 %.0f \045)",$cnorm/1000,$cpert1*100);
+ $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";
} elsif ($k == 2) {
- #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 "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 "awk '{print \$1,\$2,\$4/1000}' $dfile | xyz2grd -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile1b $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 "psscale -C${cptfile1a} $Dscale $Bscale1a -K -O -V >> $psfile\n";
+ print CSH "psscale -C${cptfile1b} $Dscale2 $Bscale1b -K -O -V >> $psfile\n";
+ #$slab1 = sprintf("ln(%s / %.2f)",$mtit,$cnorm/1000);
+ $slab1 = sprintf("(%.2f \261 %.0f \045)",$cnorm/1000,$cpert1*100);
+ $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";
@@ -325,7 +336,7 @@
print CSH "psscale -C${cptfile2} $Dscale $Bscale2 -K -O -V >> $psfile\n";
$slab1 = "ln($smodel2 / $smodel1)";
- 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 $slab1\nEOF\n";
}
}
print CSH "pscoast $J $R $coast_info -K -O -V >> $psfile \n";
@@ -424,7 +435,7 @@
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($ixv==1) {print CSH "gv $psfile &\n";}
if($ipdf==1) {print CSH "ps2pdf $psfile\n";}
} # loop over p
@@ -438,6 +449,8 @@
# print "convert $psfile $jpgfile \n";
# system("convert $psfile -rotate $rotangle $jpgfile");
# if ($ipdf==1) {system("ps2pdf $psfile")}
-# if ($ixv==1) {system("ghostview $psfile &");}
+# if ($ixv==1) {system("gv $psfile &");}
+system("gv $psfile &");
+
#==================================================
Deleted: 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 2009-05-22 13:18:45 UTC (rev 15032)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models_three.pl 2009-05-22 13:36:23 UTC (rev 15033)
@@ -1,416 +0,0 @@
-#!/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 &");}
-
-#==================================================
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-05-22 13:18:45 UTC (rev 15032)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models.pl 2009-05-22 13:36:23 UTC (rev 15033)
@@ -21,6 +21,10 @@
#
# plot_vert_models.pl 1/3 m00/m16 -40/5 1 0.15 3.0 2 1.5 # Vs
#
+#--------------------
+# GJI figures: 15,25,2,5,10,9,6,79,19,23,7,49,11
+# plot_vert_models.pl 11/11 m00/m16 -40/5 1 0.10 3.0 1 1.5
+#
#==========================================================
if (@ARGV < 8) {die("Usage: plot_vert_models.pl xxx\n");}
@@ -29,9 +33,9 @@
$icolor = 1;
$ivertlab = 0; # vertical exaggeration label
$isimbox = 1; # box around simulation region
-$iletter = 1; # label for publication figures
-$letter = "b";
-$imfinal = 1; # plot a figure with only m16 and the map
+$iletter = 0; # label for publication figures
+$letter = "a";
+$imfinal = 1; # plot an additional figure with only m16 and the map
($pmin,$pmax) = split(/\//,$pinds);
($smodel1,$smodel2) = split(/\//,$smodels);
@@ -169,7 +173,7 @@
$cshfile = "plot_vert_models.csh";
open(CSH,">$cshfile");
-print CSH "gmtset COLOR_NAN $sky_color PAPER_MEDIA letter BASEMAP_TYPE plain PLOT_DEGREE_FORMAT D TICK_LENGTH $tlen 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";
+print CSH "gmtset COLOR_NAN $sky_color PAPER_MEDIA letter MEASURE_UNIT inch BASEMAP_TYPE plain PLOT_DEGREE_FORMAT D TICK_LENGTH $tlen 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";
@iBs = (5,8,4); # what sides of each cross section to show tick marks
@@ -579,8 +583,6 @@
if ($imfinal == 1) {
- #print CSH "gmtset COLOR_NAN $sky_color 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";
-
@iBs = (5,8,4);
#for ($p = $pmin; $p <= $pmax; $p ++ ) {
@@ -622,14 +624,27 @@
#-------------------
# absolute origins -- these depend on the WIDTH of the cross section
+ $irow = 1;
+
+ if($irow==1) {
+
+ # map
+ $xgap_lab_map = 0.3;
+ $xgap_xc_map = 0.5;
+ $x0_map = $x2 - $widmap - $xgap_xc_map;
+ $y0_map = $y2;
+ $omap = "-Xa${x0_map} -Ya${y0_map}";
+
# overall label
- $x0_label = $x2;
- $y0_label = $y0_xc + $heightxc;
+ $x0_label = $x0_map - $xgap_lab_map;
+ $y0_label = $y0_map + $heightxc;
$olab = "-Xa${x0_label} -Ya${y0_label}";
- $x0_label2 = $x2;
- $y0_label2 = $y0_xc + 0.8*$heightxc;
+ $x0_label2 = $x0_label;
+ $y0_label2 = $y0_map + 0.8*$heightxc;
$olab2 = "-Xa${x0_label2} -Ya${y0_label2}";
+ } else {
+
# map
$xgap_lab_map = 0.4;
$ygap_xc_map = 0.2;
@@ -637,6 +652,16 @@
$y0_map = $y2 + $heightxc + $ygap_xc_map;
$omap = "-Xa${x0_map} -Ya${y0_map}";
+ # overall label
+ $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}";
+
+ }
+
#---------------------
print CSH "gmtset TICK_LENGTH $tlen_map\n"; # ticks for map
@@ -740,7 +765,7 @@
print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$saf $zmax0\n$saf $zminf0\nEOF\n";
print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$gf $zmax0\n$gf $zminf0\nEOF\n";
print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$sgf $zmax0\n$sgf $zminf0\nEOF\n";
- print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$mcf $zmax0\n$mcf $zminf0\nEOF\n";
+ #print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$mcf $zmax0\n$mcf $zminf0\nEOF\n";
print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$syf $zmax0\n$syf $zminf0\nEOF\n";
print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$crf $zmax0\n$crf $zminf0\nEOF\n";
print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$ef $zmax0\n$ef $zminf0\nEOF\n";
@@ -767,7 +792,7 @@
if($saf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$saf $zmax0 $flsize 0 $fontno $talign SA\nEOF\n";}
if($gf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$gf $zmax0 $flsize 0 $fontno $talign G\nEOF\n";}
if($sgf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$sgf $zmax0 $flsize 0 $fontno $talign SG\nEOF\n";}
- if($mcf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$mcf $zmax0 $flsize 0 $fontno $talign MC\nEOF\n";}
+ #if($mcf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$mcf $zmax0 $flsize 0 $fontno $talign MC\nEOF\n";}
if($syf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$syf $zmax0 $flsize 0 $fontno $talign SY\nEOF\n";}
if($crf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$crf $zmax0 $flsize 0 $fontno $talign CR\nEOF\n";}
if($ef != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$ef $zmax0 $flsize 0 $fontno $talign E\nEOF\n";}
@@ -820,6 +845,7 @@
# print "convert $psfile $jpgfile \n";
# system("convert $psfile -rotate $rotangle $jpgfile");
# if ($ipdf==1) {system("ps2pdf $psfile")}
-# if ($ixv==1) {system("ghostview $psfile &");}
+# if ($ixv==1) {system("gv $psfile &");}
+system("gv $psfile &");
#==================================================
Deleted: 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-05-22 13:18:45 UTC (rev 15032)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models_basemap.pl 2009-05-22 13:36:23 UTC (rev 15033)
@@ -1,291 +0,0 @@
-#!/usr/bin/perl -w
-
-#==========================================================
-#
-# plot_vert_models_basemap.pl
-# Carl Tape
-# 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 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");}
- at pinds = @ARGV;
-$nump = @pinds;
-$stp = sprintf("%3.3i",$nump);
-
-$irun = 1; # 1, 2 (SAF), 3 (Garlock)
-$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
-
-#($pmin,$pmax) = split(/\//,$pinds);
-#($smodel1,$smodel2) = split(/\//,$smodels);
-#$smodeltag = "${smodel2}_${smodel1}";
-#if($ivel==1) {$modlab = "vs"; $mtit = "Vs"; $vmin = 2.0; $vmax = 4.21; $icol = 3; $vtick = 1;}
-#elsif($ivel==2) {$modlab = "vb"; $mtit = "Vb"; $vmin = 2.0; $vmax = 5.5; $icol = 6; $vtick = 1;}
-#elsif($ivel==3) {$modlab = "vp"; $mtit = "Vp"; $vmin = 2.5; $vmax = 7.5; $icol = 9; $vtick = 1;}
-#($zmin0,$zmax0) = split(/\//,$zran);
-
-# directory containing the data files
-$stirun = sprintf("%2.2i",$irun);
-$dirdat = "INPUT/vert_${stirun}";
-
-
-#---------------------------------------------------------
-
-# subtitles
-#@titles = ("$mtit $smodel1","$mtit $smodel2","ln(${smodel2} / ${smodel1})");
-
-# parameters for color scales
-#$cpert1 = 0.2;
-#$cpert2 = 0.2;
-$scale_color = 65.0;
-$colorbar = "seis";
-
-# KEY: resolution of color plots
-$interp = "-I0.5/0.5 -S1.5"; # nearneighbor
-#$interp = "-I1.0"; # xyz2grd
-$grdfile = "temp.grd";
-
- at Bopts = ("WESN","Wesn","wesN","wEsn","weSn","WesN","wEsN","wESn","WeSn","WEsn","weSN","wESN","WESn","WeSN","WEsN","wesn");
-
-$dir0 = "/home/carltape/gmt";
-$dir = "socal_2005";
-
-$fault_file = "$dir0/faults/jennings_more.xy";
-$fault_file2 = "$dir0/faults/socal_faults_xc.lonlat";
-#$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";
-
-# socal map
-#$tick1 = 100; $tick2 = 100; # no ticks for map
-$tick1 = 1; $tick2 = 0.2;
-#$Rmap = "-R-121.2/-114.8/32.3/36.7";
-$Rmap = "-R-122/-114/32/37";
-
-$topo_labels = "$dir0/socal_2005/socal_topo_labs.xyz";
-$fault_labels = "$dir0/socal_2005/socal_fault_labs.xyz";
-$coast_res = "-Df -A0/0/4";
-
-$iportrait = 1;
-if($iportrait == 1){$orient = "-P"; $rotangle = 0}
-else {$orient = " "; $rotangle = 90}
-
-# plotting specifications
-$fsize0 = "14";
-$fsize1 = "12";
-$fsize2 = "11";
-$fsize3 = "10";
-$fontno = "1"; # 1 or 4
-$tlen = "6p";
-$fpen = "1.5p";
-$tpen = "1.5p";
-
-# 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_res = "-Df -A0/0/4";
-#$coast_info = "${coast_res} -W0.75p -Na/0.75p";
-$coast_info = "${coast_res} -W1p -Na/1.0p -S220/255/255 -C220/255/255 -G200"; # -I1/1.0p for rivers
-$coast_info2 = "${coast_res} -S220/255/255 -C220/255/255 -G200"; # -I1/1.0p for rivers
-$fault_info = "-M -W0.5p,255/0/0";
-$fault_infoK = "-M -W0.5p,0/0/0";
-$fault_info2 = "-M -W2p,255/0/0";
-
-$sky_color = "100/100/100";
-$city_info = "-Sc10p -W1.0p/0/0/0 -G255/255/0";
-
-#$colinfo = "-W1p,0/0/0 -G255/255/255";
-$colinfo_map = "-W1.5p,0/0/0 -G0/255/255";
-$srcinfo_map = "-Sa20p $colinfo_map -N";
-$recinfo_map = "-Si20p $colinfo_map -N";
-$srcinfo_mapinset = "-Sa10p -W0.5p,0/0/0 -G0/255/255 -N";
-$recinfo_mapinset = "-Si10p -W0.5p,0/0/0 -G0/255/255 -N";
-$colinfo_xc = "-W2p,0/0/0";
-$srcinfo_xc = "-Sa20p $colinfo_xc -N";
-$recinfo_xc = "-Si20p $colinfo_map -N";
-
-$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 = 1; $ipdf = 0;
-$ixv = 0; $ipdf = 1;
-
-# plot title
-$J_title = "-JM1";
-$R_title = "-R0/1/0/1";
-$fsize_title = $fsize0;
-
-#---------------------------------------------------------
-
-# load values and check that the files exist
-($nump0,undef,undef) = split(" ",`ls -1 ${dirdat}/*path | wc`);
-if($nump0 == 0) {die("No ray paths available")}
-#if($pmax > $nump) {$pmax = $nump;}
-#if($pmin < 1) {$pmin = 1;}
-print "\nNumber of total cuts is $nump0\n";
-
-# load ALL fault positions along each profile
- #$ffile = "$dirdat/ALL_rays_fault_positions";
- $ffile = "$dirdat/ALL_rays_fault_positions_mod";
- if (not -f $ffile) {die("Check if ffile $ffile exist or not\n");}
- open(IN,$ffile); @faults = <IN>; close(IN);
-
- $widmap = 6.5; # width of map
- $Jmap = "-JM${widmap}";
- $B0map = "-Ba${tick1}f${tick2}d:.\" \":";
- #$Bmap = "$B0map".$Bopts[7];
- $Bmap = "$B0map".$Bopts[0];
-
- # parameters controlling the placement of subplots and labels
- $omap = "-Xa1 -Ya3";
-
- $fname = "vert_${stirun}_basemap_${stp}";
- $psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
-
-#==============================================
-
-$cshfile = "plot_vert_models_basemap.csh";
-open(CSH,">$cshfile");
-print CSH "gmtset COLOR_NAN $sky_color PAPER_MEDIA letter BASEMAP_TYPE plain PLOT_DEGREE_FORMAT D TICK_LENGTH $tlen 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";
-
- print CSH "psbasemap $Jmap $Rmap $Bmap -K -V $orient $omap > $psfile\n"; # START
-
- print CSH "pscoast $Jmap $Rmap $coast_info -K -O -V $omap >> $psfile \n";
- if ($ishelf==1) {
- print CSH "awk '{print \$2,\$1}' ${shelf_file} | psxy $Jmap $Rmap $Wshelf -K -O -V $omap >> $psfile\n";
- }
-
- print CSH "psxy ${fault_file} $Jmap $Rmap $fault_info -K -V -O $omap >> $psfile \n";
- print CSH "psxy ${fault_file2} $Jmap $Rmap $fault_info2 -K -V -O $omap >> $psfile \n";
-
- if ($isimbox==1) {
- $outer_boundary = "${dir_sources}/SPECFEM_outer_points.dat";
- if (not -f ${outer_boundary}) {
- die("Check if outer_boundary ${outer_boundary} exist or not\n");
- }
- print CSH "psxy ${outer_boundary} $Jmap $Rmap -W2p,0/0/0 -K -O -V $omap >> $psfile\n";
- }
-
-#==============================================
-
-# modify the height of each cross-section based on the length
-#for ($p = $pmin; $p <= $pmax; $p ++ ) {
-for ($ip = 1; $ip <= $nump; $ip ++ ) {
- $p = $pinds[$ip-1];
- $stip = sprintf("%3.3i",$p);
- print "-- $ip -- $nump -- $p --\n";
-
- # source and receiver points -- NOT necessarily the endpoints of the ray path
- # COLUMNS: slons(ii),slats(ii),0,rlons(ii),rlats(ii),deg2km(dq),azq,iflip
- $rfile = "${dirdat}/vert_xc_${stip}_A_B";
- if (not -f $rfile) {die("Check if rfile $rfile exist or not\n");}
- open(IN,$rfile); @line = <IN>; close(IN);
- ($slon,$slat,$sdist,$sdep,$rlon,$rlat,$rdist,$rdep,$raz,$iflip) = split(" ",$line[0]);
-
- # positions of six faults along each profile
- # (1) SAF, (2) GF, (3) SGF, (4) MCF, (5) SYF, (6) CRF
- (undef,$eid,$stanet,$azi1,$azi2,$saf,$gf,$sgf,$mcf,$syf,$crf) = split(" ",$faults[$p-1]);
- ($sta,$net) = split("\\.",$stanet);
- $sazi1 = sprintf("%3.3i",$azi1);
- $sazi2 = sprintf("%3.3i",$azi2);
- print "-- $stanet -- $sta -- $net -- $eid -- $sazi1 -- $sazi2 --\n";
-
- #---------------------------------------------------------
-
- # plot ray path
- $pfile = "${dirdat}/vert_xc_${stip}_ray_path";
- if (not -f $pfile) {
- die("Check if pfile $pfile exist or not\n");
- }
- print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $Jmap $Rmap -W6p,0/0/0 -K -O -V $omap >> $psfile\n";
- 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");
- }
-
- 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
- #print CSH "awk '{print \$2,\$1}' $cities | psxy ${city_info} $Jmap $Rmap -K -V -O >> $psfile \n";
- #print CSH "awk '{print \$2,\$1,12,0,1,\"LB\",\$3}' $cities | pstext $textinfo3 -D0.05/0.05 $Jmap $Rmap -K -V -O >> $psfile \n";
-
- #print CSH "pstext ${topo_labels} $Jmap $Rmap $textinfo -N -K -V -O >> $psfile\n";
- #print CSH "pstext ${fault_labels} $textinfo2 $Jmap $Rmap -K -V -O >> $psfile \n";
-
-# # 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 {
-# $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";
-# }
-
-} # loop over p
-
- # fault labels
- $fault_labs = "$dir0/faults/socal_fault_labs.lonlat";
- $textinfo = "-N -C4p -W255/255/255o,1.0p,255/0/0,solid -G0/0/0";
- print CSH "awk '{print \$1,\$2,10,0,1,\"CM\",\$3}' ${fault_labs} | pstext $textinfo $Jmap $Rmap -K -V -O $omap >> $psfile\n";
-
- print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n 10 10 16 0 $fontno CM \nEOF\n";
- print CSH "psbasemap $Jmap $Rmap $Bmap -O -V $omap >> $psfile\n"; # FINISH (no K)
-
- #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";}
-
-#------------------------------------
-
-close (CSH);
-system("csh -f $cshfile");
-
-system("ghostview $psfile &")
-
-#==================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_one.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_one.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_one.pl 2009-05-22 13:36:23 UTC (rev 15033)
@@ -0,0 +1,420 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# plot_horz_models_one.pl
+# Carl Tape
+# 21-May-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 m16 cross-sections only.
+#
+# FIGURES:
+# plot_horz_models_one.pl 11/11 m00/m16 1 0/0/1/0.2 # Vs
+#
+#==========================================================
+
+if (@ARGV < 4) {die("Usage: plot_horz_models_one.pl xxx\n");}
+($pinds,$smodels,$ivs,$parms) = @ARGV;
+
+$icolor = 1;
+($pmin,$pmax) = split(/\//,$pinds);
+($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 = "b";
+$iextra = 1; # extra information on figures
+
+$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>; $nump = @lines;
+for ($p = 1; $p <= $nump; $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 $nump\n";
+if($pmax > $nump) {$pmax = $nump;}
+if($pmin < 1) {$pmin = 1;}
+
+# subtitles
+ at titles = ("$mtit model $smodel2");
+ at labs = ("Vs $smodel2");
+
+# 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";
+$interp = "-I1.5m";
+$grdfile = "temp.grd";
+
+# position of titles and labels
+$xtx1 = 0.5; $ytx1 = 1.00;
+$xtx2 = -0.2; $ytx2 = 0.4;
+#$xtx3 = 1.0; $ytx3 = 0.6;
+$xtx4 = 1.04; $ytx4 = 0.6;
+$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 = 4; # 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 = 1;
+ $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_one";
+ $origin = "-X2i -Y4i";
+
+ # 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
+$Dwid = 0.15;
+$Dlen = $wid*0.5;
+$Dx = $wid + 0.25;
+$Dy = $Dlen/2;
+#$Dy2 = $Dy-$Dwid;
+$Dy2 = $Dy;
+$Dscale = "-D$Dx/$Dy/$Dlen/${Dwid}";
+$Dscale2 = "-D$Dx/$Dy2/$Dlen/${Dwid}";
+
+# plotting specifications
+$fsize0 = "16";
+$fsize1 = "12";
+$fsize2 = "11";
+$fsize3 = "8";
+$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_one.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 ($p = $pmin; $p <= $pmax; $p ++ ) {
+
+ $stip = sprintf("%3.3i",$p);
+ $fname = "horz_xc_${modlab}_${smodeltag}_${stip}_extra${iextra}_single";
+ $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";
+
+ # colorpoint file for m00 and m16 -- VARIES with depth -- in KILOMETERS
+ # THIS IS FOR THE COLORBAR ONLY
+ $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";
+ $Bscale1b = sprintf("-B%2.2ef0.1:\" \": -E10p",$ctick1b_km);
+
+ # 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";
+
+ # ticks for colorbars
+ $Bscale1a = sprintf("-Ba%2.2ef0.05:\" \": -E10p -A",$cpert1);
+ $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");}
+
+
+ for ($k = 1; $k <= 1; $k ++ ) {
+
+ $shift = $shifts[$k-1];
+ $iB = 8;
+ $B = "$B0".$Bopts[$iB];
+ $title = $titles[$k-1];
+
+ if ($k==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 "awk '{print \$1,\$2,\$4)}' $dfile | xyz2grd -G$grdfile $R $interp\n";
+ #print CSH "grdimage $grdfile -C$cptfile1a $J -Q -K -O -V >> $psfile\n";
+
+ print CSH "awk '{print \$1,\$2,\$4/1000}' $dfile | xyz2grd -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile1b $J -Q -K -O -V >> $psfile\n";
+
+ #print CSH "psscale -C${cptfile1a} $Dscale $Bscale1a -K -O -V >> $psfile\n";
+ print CSH "psscale -C${cptfile1b} $Dscale2 $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 14 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 ($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 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";
+ }
+ }
+
+ # fault labels
+ $fault_labs = "$dir0/faults/socal_fault_labs_red.lonlat";
+ $textinfo = "-N -C2p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ print CSH "awk '{print \$1,\$2,10,0,1,\"CM\",\$3}' ${fault_labs} | pstext $textinfo $J $R -K -V -O >> $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 = $labs[$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$xtx6 $ytx6 14 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 k
+
+ 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 "gv $psfile &\n";}
+ if($ipdf==1) {print CSH "ps2pdf $psfile\n";}
+
+} # loop over p
+
+#------------------------------------
+# 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");
+system("gv $psfile &");
+
+#==================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_one.pl
___________________________________________________________________
Name: svn:executable
+ *
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_three.pl (from rev 14984, seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models_three.pl)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_three.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_three.pl 2009-05-22 13:36:23 UTC (rev 15033)
@@ -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/scripts/plot_horz_models_three.pl
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_vert_models_basemap.pl (from rev 14984, seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models_basemap.pl)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_vert_models_basemap.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_vert_models_basemap.pl 2009-05-22 13:36:23 UTC (rev 15033)
@@ -0,0 +1,291 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# plot_vert_models_basemap.pl
+# Carl Tape
+# 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 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");}
+ at pinds = @ARGV;
+$nump = @pinds;
+$stp = sprintf("%3.3i",$nump);
+
+$irun = 2; # KEY: 1, 2 (SAF), 3 (Garlock)
+$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
+
+#($pmin,$pmax) = split(/\//,$pinds);
+#($smodel1,$smodel2) = split(/\//,$smodels);
+#$smodeltag = "${smodel2}_${smodel1}";
+#if($ivel==1) {$modlab = "vs"; $mtit = "Vs"; $vmin = 2.0; $vmax = 4.21; $icol = 3; $vtick = 1;}
+#elsif($ivel==2) {$modlab = "vb"; $mtit = "Vb"; $vmin = 2.0; $vmax = 5.5; $icol = 6; $vtick = 1;}
+#elsif($ivel==3) {$modlab = "vp"; $mtit = "Vp"; $vmin = 2.5; $vmax = 7.5; $icol = 9; $vtick = 1;}
+#($zmin0,$zmax0) = split(/\//,$zran);
+
+# directory containing the data files
+$stirun = sprintf("%2.2i",$irun);
+$dirdat = "INPUT/vert_${stirun}";
+
+
+#---------------------------------------------------------
+
+# subtitles
+#@titles = ("$mtit $smodel1","$mtit $smodel2","ln(${smodel2} / ${smodel1})");
+
+# parameters for color scales
+#$cpert1 = 0.2;
+#$cpert2 = 0.2;
+$scale_color = 65.0;
+$colorbar = "seis";
+
+# KEY: resolution of color plots
+$interp = "-I0.5/0.5 -S1.5"; # nearneighbor
+#$interp = "-I1.0"; # xyz2grd
+$grdfile = "temp.grd";
+
+ at Bopts = ("WESN","Wesn","wesN","wEsn","weSn","WesN","wEsN","wESn","WeSn","WEsn","weSN","wESN","WESn","WeSN","WEsN","wesn");
+
+$dir0 = "/home/carltape/gmt";
+$dir = "socal_2005";
+
+$fault_file = "$dir0/faults/jennings_more.xy";
+$fault_file2 = "$dir0/faults/socal_faults_xc.lonlat";
+#$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";
+
+# socal map
+#$tick1 = 100; $tick2 = 100; # no ticks for map
+$tick1 = 1; $tick2 = 0.2;
+#$Rmap = "-R-121.2/-114.8/32.3/36.7";
+$Rmap = "-R-122/-114/32/37";
+
+$topo_labels = "$dir0/socal_2005/socal_topo_labs.xyz";
+$fault_labels = "$dir0/socal_2005/socal_fault_labs.xyz";
+$coast_res = "-Df -A0/0/4";
+
+$iportrait = 1;
+if($iportrait == 1){$orient = "-P"; $rotangle = 0}
+else {$orient = " "; $rotangle = 90}
+
+# plotting specifications
+$fsize0 = "14";
+$fsize1 = "12";
+$fsize2 = "11";
+$fsize3 = "10";
+$fontno = "1"; # 1 or 4
+$tlen = "6p";
+$fpen = "1.5p";
+$tpen = "1.5p";
+
+# 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_res = "-Df -A0/0/4";
+#$coast_info = "${coast_res} -W0.75p -Na/0.75p";
+$coast_info = "${coast_res} -W1p -Na/1.0p -S220/255/255 -C220/255/255 -G200"; # -I1/1.0p for rivers
+$coast_info2 = "${coast_res} -S220/255/255 -C220/255/255 -G200"; # -I1/1.0p for rivers
+$fault_info = "-M -W0.5p,255/0/0";
+$fault_infoK = "-M -W0.5p,0/0/0";
+$fault_info2 = "-M -W2p,255/0/0";
+
+$sky_color = "100/100/100";
+$city_info = "-Sc10p -W1.0p/0/0/0 -G255/255/0";
+
+#$colinfo = "-W1p,0/0/0 -G255/255/255";
+$colinfo_map = "-W1.5p,0/0/0 -G0/255/255";
+$srcinfo_map = "-Sa20p $colinfo_map -N";
+$recinfo_map = "-Si20p $colinfo_map -N";
+$srcinfo_mapinset = "-Sa10p -W0.5p,0/0/0 -G0/255/255 -N";
+$recinfo_mapinset = "-Si10p -W0.5p,0/0/0 -G0/255/255 -N";
+$colinfo_xc = "-W2p,0/0/0";
+$srcinfo_xc = "-Sa20p $colinfo_xc -N";
+$recinfo_xc = "-Si20p $colinfo_map -N";
+
+$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 = 1; $ipdf = 0;
+$ixv = 0; $ipdf = 1;
+
+# plot title
+$J_title = "-JM1";
+$R_title = "-R0/1/0/1";
+$fsize_title = $fsize0;
+
+#---------------------------------------------------------
+
+# load values and check that the files exist
+($nump0,undef,undef) = split(" ",`ls -1 ${dirdat}/*path | wc`);
+if($nump0 == 0) {die("No ray paths available")}
+#if($pmax > $nump) {$pmax = $nump;}
+#if($pmin < 1) {$pmin = 1;}
+print "\nNumber of total cuts is $nump0\n";
+
+# load ALL fault positions along each profile
+ #$ffile = "$dirdat/ALL_rays_fault_positions";
+ $ffile = "$dirdat/ALL_rays_fault_positions_mod";
+ if (not -f $ffile) {die("Check if ffile $ffile exist or not\n");}
+ open(IN,$ffile); @faults = <IN>; close(IN);
+
+ $widmap = 6.5; # width of map
+ $Jmap = "-JM${widmap}";
+ $B0map = "-Ba${tick1}f${tick2}d:.\" \":";
+ #$Bmap = "$B0map".$Bopts[7];
+ $Bmap = "$B0map".$Bopts[0];
+
+ # parameters controlling the placement of subplots and labels
+ $omap = "-Xa1 -Ya3";
+
+ $fname = "vert_${stirun}_basemap_${stp}";
+ $psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
+
+#==============================================
+
+$cshfile = "plot_vert_models_basemap.csh";
+open(CSH,">$cshfile");
+print CSH "gmtset COLOR_NAN $sky_color PAPER_MEDIA letter BASEMAP_TYPE plain PLOT_DEGREE_FORMAT D TICK_LENGTH $tlen 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";
+
+ print CSH "psbasemap $Jmap $Rmap $Bmap -K -V $orient $omap > $psfile\n"; # START
+
+ print CSH "pscoast $Jmap $Rmap $coast_info -K -O -V $omap >> $psfile \n";
+ if ($ishelf==1) {
+ print CSH "awk '{print \$2,\$1}' ${shelf_file} | psxy $Jmap $Rmap $Wshelf -K -O -V $omap >> $psfile\n";
+ }
+
+ print CSH "psxy ${fault_file} $Jmap $Rmap $fault_info -K -V -O $omap >> $psfile \n";
+ print CSH "psxy ${fault_file2} $Jmap $Rmap $fault_info2 -K -V -O $omap >> $psfile \n";
+
+ if ($isimbox==1) {
+ $outer_boundary = "${dir_sources}/SPECFEM_outer_points.dat";
+ if (not -f ${outer_boundary}) {
+ die("Check if outer_boundary ${outer_boundary} exist or not\n");
+ }
+ print CSH "psxy ${outer_boundary} $Jmap $Rmap -W2p,0/0/0 -K -O -V $omap >> $psfile\n";
+ }
+
+#==============================================
+
+# modify the height of each cross-section based on the length
+#for ($p = $pmin; $p <= $pmax; $p ++ ) {
+for ($ip = 1; $ip <= $nump; $ip ++ ) {
+ $p = $pinds[$ip-1];
+ $stip = sprintf("%3.3i",$p);
+ print "-- $ip -- $nump -- $p --\n";
+
+ # source and receiver points -- NOT necessarily the endpoints of the ray path
+ # COLUMNS: slons(ii),slats(ii),0,rlons(ii),rlats(ii),deg2km(dq),azq,iflip
+ $rfile = "${dirdat}/vert_xc_${stip}_A_B";
+ if (not -f $rfile) {die("Check if rfile $rfile exist or not\n");}
+ open(IN,$rfile); @line = <IN>; close(IN);
+ ($slon,$slat,$sdist,$sdep,$rlon,$rlat,$rdist,$rdep,$raz,$iflip) = split(" ",$line[0]);
+
+ # positions of six faults along each profile
+ # (1) SAF, (2) GF, (3) SGF, (4) MCF, (5) SYF, (6) CRF
+ (undef,$eid,$stanet,$azi1,$azi2,$saf,$gf,$sgf,$mcf,$syf,$crf) = split(" ",$faults[$p-1]);
+ ($sta,$net) = split("\\.",$stanet);
+ $sazi1 = sprintf("%3.3i",$azi1);
+ $sazi2 = sprintf("%3.3i",$azi2);
+ print "-- $stanet -- $sta -- $net -- $eid -- $sazi1 -- $sazi2 --\n";
+
+ #---------------------------------------------------------
+
+ # plot ray path
+ $pfile = "${dirdat}/vert_xc_${stip}_ray_path";
+ if (not -f $pfile) {
+ die("Check if pfile $pfile exist or not\n");
+ }
+ print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $Jmap $Rmap -W6p,0/0/0 -K -O -V $omap >> $psfile\n";
+ 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");
+ }
+
+ 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
+ #print CSH "awk '{print \$2,\$1}' $cities | psxy ${city_info} $Jmap $Rmap -K -V -O >> $psfile \n";
+ #print CSH "awk '{print \$2,\$1,12,0,1,\"LB\",\$3}' $cities | pstext $textinfo3 -D0.05/0.05 $Jmap $Rmap -K -V -O >> $psfile \n";
+
+ #print CSH "pstext ${topo_labels} $Jmap $Rmap $textinfo -N -K -V -O >> $psfile\n";
+ #print CSH "pstext ${fault_labels} $textinfo2 $Jmap $Rmap -K -V -O >> $psfile \n";
+
+# # 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 {
+# $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";
+# }
+
+} # loop over p
+
+ # fault labels
+ $fault_labs = "$dir0/faults/socal_fault_labs.lonlat";
+ $textinfo = "-N -C4p -W255/255/255o,1.0p,255/0/0,solid -G0/0/0";
+ print CSH "awk '{print \$1,\$2,10,0,1,\"CM\",\$3}' ${fault_labs} | pstext $textinfo $Jmap $Rmap -K -V -O $omap >> $psfile\n";
+
+ print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n 10 10 16 0 $fontno CM \nEOF\n";
+ print CSH "psbasemap $Jmap $Rmap $Bmap -O -V $omap >> $psfile\n"; # FINISH (no K)
+
+ #if($ixv==1) {print CSH "convert $psfile -rotate $rotangle $jpgfile\n";}
+ if ($ixv==1) {print CSH "gv $psfile &\n";}
+ if ($ipdf==1) {print CSH "ps2pdf $psfile\n";}
+
+#------------------------------------
+
+close (CSH);
+system("csh -f $cshfile");
+
+system("gv $psfile &")
+
+#==================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_vert_models_basemap.pl
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:mergeinfo
+
More information about the CIG-COMMITS
mailing list