plane_percvp.gmt

这是一个TOMO TO GMT的源码

#!/bin/csh
#
# to plot 
#  percentage velocity change 
# in a horizontal depth section  
# 
# take output from 'tomo2gmt', stored in the VP directory,
# and named planeXXXX.X.gmt
# 
# This script has been set up to plot 4 sections per page.
#
# This script has been adapted for GMT3.3.1 and measure unit is inch!
#  
#				(c) [e-mail]

# Set Istring to correspond to your grid spacing (lon/lat) in the file 
# planeXXXX.X.gmt in DEGREES. GMT3.3.1 is very sensitive to the grid  
# spacing values, so the best way to calculate them is to use the minimum
# and maximum longitude/latitude in the planeXXXX.X.gmt file, divide the
# difference by the corresponding distance in km and multiply it by 
# the grid spacing in km. Be sure to use at least 8 digits!

# for 2x2 km gridding
set Istring=0.02292750/0.01801750

# general settings

# !! be sure to set paths to your directory !!
#  wpath = working directory with GMT scripts
#  inpath = directory where to find output files by tomo2gmt
 
set wpath=/ek1/tomoled_soft/TOMO2GMT/PACKAGE
set inpath=$wpath/example_tomores/VP

# boundaries (lon/lat)
# set to min/max values in the planeXXXX.X.xyz file
set Rstring=20.0829/21.00/38.5/39.2207

# scaling, set to absolute value
set Jstring=M3

set psfile=$wpath/GMT/plane_percvp.ps

# which results are plotted (fatomo or simulps)?
set tomo=fatomo

# label for colorscale
set clabel="% Vp change, rel. to 1D initial model"

# % velocity change is the 7th column in the planeXXXX.X.xyz files
# to plot absolute velocities set pos=6
set pos=7

# color table for grdimage
set cpt=$wpath/GMT/perc_chnge_bw.syn.cpt


### now GMT-job

#plot colorscale
psscale -C$cpt -D4.5/0.5/5/0.25h -L -B:"$clabel": -K -P -X-0.25 -Y1.0 > $psfile

set cnt=0

# the following is done for each plane listed between the parentheses
# modify for your planes to be used
foreach plane(0010.0 0015.0 0020.0 0025.0)

set cnt=`expr $cnt + 1`

set infile=$inpath/plane$plane.xyz

# get plot position on the page
switch ($cnt)
 case 1:
   set xpos=1
   set ypos=6
   breaksw
 case 2:
   set xpos=4
   set ypos=0
   breaksw
 case 3:
   set xpos=-4
   set ypos=-5
   breaksw
 case 4:
   set xpos=4
   set ypos=0
   breaksw
endsw

echo --- working on plane $plane ---

# create grid using xyz2grd
nawk 'substr($1,1,1)!="#" { print $3,$4,$pos }' pos=$pos $infile | \
xyz2grd -I$Istring -N-99 -Gpvp.grd -R$Rstring -V

# image grid
grdimage pvp.grd -J$Jstring -R$Rstring -Ba0.5WenS -C$cpt -O -K -P -V -X$xpos -Y$ypos >> $psfile

# draw coastlines
#  not in this example since this is an synthetic example
#pscoast -J$Jstring -R -Df -W2 -N1/1top -O -K -P -V >>$psfile

# plot nodes with crosses
if ($tomo == "fatomo") then
 set nodefile=$inpath/plane_nodes.xyz
else
 set nodefile=$inpath/plane_nodes.xyz
endif
nawk '{print $4,$3}' $nodefile | \
psxy  -J$Jstring -R$Rstring -Sx0.08 -W0.5/120top -O -K -V >>$psfile 

# plot stations 
psxy $inpath/sta.xyz -R$Rstring -J$Jstring -St0.1 -L -G230 -O -P  -K  >>$psfile

# print depth of horizontal section at upper left corner
# adjust lon/lat position 
pstext -J$Jstring -R -O -P -K -V -N <<END>> $psfile
20.10 39.3 12 0 1 CB z = $plane km
END

end

# dummy closing file
pstext -J$Jstring -R -N -O -P <<END >>$psfile
END

gs $psfile

# cleanup 
/bin/rm pvp.grd

exit