latsec.gmt

这是一个TOMO TO GMT的源码

#!/bin/csh
#
# to plot values of 
#   1:  KHIT matrix 
#   2:  DWS
#   3:  diag. reso element (RDE)
#   4:  spread values
# in vertical depth sections along constant latitude
#
# take output from 'tomo2gmt', stored in the VP directory,
# and named latsecX.gmt
# 
#
# This script has been adapted for GMT3.3.1 and measure unit is inch!
#
#  call by: latsec.gmt <1,2,3,4>  
#					(c) [e-mail]

# Set Istring to correspond to your grid spacing (lon/z) in the file 
# latsecX.gmt in DEGREES/KM. 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 in the latsecX.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/2

# !! 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/GMT
set inpath=/ek1/tomoled_soft/TOMO2GMT/PACKAGE/example_tomores/VP

# boundaries (lon/z)
# set to min/max values in the latsecXX.xyz file
set Rstring=20.0829/21.00/-4.0/40.0

# scaling, VE ~1.0
set Jstring=x6/0.06

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

# choose color table, name of psfile, label of color scale,
# and number of column according to what is desired to be plotted

if ($tomo == "fatomo") then
switch ($1)
  case  1:
   set cpt=$wpath/grey_hit_simulps.cpt
   set psfile=$wpath/latsec.khit.ps
   set clabel="KHIT"
   set pos=10
   breaksw
  case  2:
   set cpt=$wpath/grey_dws.cpt
   set psfile=$wpath/latec.dws.ps
   set clabel="Derivative Weigthed Sum (DWS)"
   set pos=11
   breaksw
  case  3:
   set cpt=$wpath/grey1.cpt
   set psfile=$wpath/latsec.rde.ps
   set clabel="Resolution Diag. Element"
   set pos=12
   breaksw
  case  4:
   set cpt=$wpath/sprd2.cpt
   set psfile=$wpath/latsec.sprd.ps
   set clabel="spread value"
   set pos=13
   breaksw
endsw
else
switch ($1)
  case  1:
   set cpt=$wpath/grey_hit_simulps.cpt
   set psfile=$wpath/latec.khit.sml.ps
   set clabel="no. of rays (KHIT)"
   set pos=10
   breaksw
  case  2:
   set cpt=$wpath/grey_dws.cpt
   set psfile=$wpath/latsec.dws.sml.ps
   set clabel="Derivative Weigthed Sum (DWS)"
   set pos=11
   breaksw
  case  3:
   set cpt=$wpath/grey1.cpt
   set psfile=$wpath/latsec.rde.sml.ps
   set clabel="Reso. Diag. Element"
   set pos=12
   breaksw
  case  4:
   set cpt=$wpath/sprd2.cpt
   set psfile=$wpath/latsec.sprd.sml.ps
   set clabel="spread value"
   set pos=13
 endsw
endif

### now start GMT job

set latsec=38.25

set infile1=$inpath/latsec$latsec.xyz

# create grid 
nawk 'substr($1,1,1)!="#" {print $3,$5,$pos}' pos=$pos $infile1 | \
xyz2grd -I$Istring -N-99 -Glatsec.grd -R$Rgrid -V

# image grid
grdimage latsec.grd -J$Jstring -R$Rstring -Baf0.5:"latitude [degrees]":/a20f10:"depth [km]":WeNs -C$cpt -K -P -V -X1.5 -Y1.5 > $psfile


# plot nodes with crosses
psxy $nodefile -J$Jstring -R -Sx0.05 -N -W0.5top -O -P -K -V >>$psfile 

# plot title at bottom, adjust to your model size
pstext -J$Jstring -R -O -P -K -V -N <<END >>$psfile
20.5  38 12 0 1 CM Latitude Depth Section along $lat S
END

#plot colorscale
psscale -C$cpt -D4.5/1.5/5/0.25h -L -B:"$clabel": -O -P -X-2 -Y-2 >> $psfile


gs $psfile

# cleanup 
/bin/rm latsec.grd

exit