profile.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 user defined profiles
#
# take output from 'tomo2gmt', stored in the VP directory,
# and named profileX.gmt
# 
#
# This script has been adapted for GMT3.3.1 and measure unit is inch!
#
#  call by: profile.gmt <1,2,3,4>  
#					(c) [e-mail]


# Set Istring to correspond to your grid spacing (in(=km along profile)/z) in the file 
# latsecX.gmt in KM/KM. 

# for 2x2 km gridding
set Istring=2.0/2.0

# !! 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 (in/z)
# set to min/max values in the latsecXX.xyz file
set Rstring=0/80/-4.0/40.0

# scaling, VE = 1.0
set Jstring=x0.06/-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

switch ($1)
  case  1:
   set cpt=$wpath/grey_hit_simulps.cpt
   set psfile=$wpath/latsec.khit.ps
   set clabel="KHIT"
   set pos=11
   breaksw
  case  2:
   set cpt=$wpath/grey_dws.cpt
   set psfile=$wpath/latec.dws.ps
   set clabel="Derivative Weigthed Sum (DWS)"
   set pos=12
   breaksw
  case  3:
   set cpt=$wpath/grey1.cpt
   set psfile=$wpath/latsec.rde.ps
   set clabel="Resolution Diag. Element"
   set pos=13
   breaksw
  case  4:
   set cpt=$wpath/sprd2.cpt
   set psfile=$wpath/latsec.sprd.ps
   set clabel="spread value"
   set pos=14
   breaksw
endsw


### now start GMT job

set prof=01

set infile1=$inpath/profile$prof.xyz

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

# image grid
grdimage prof.grd -J$Jstring -R$Rstring -Ba10f5:"km along profile":/a20f10:"depth [km]":WeNs -C$cpt -K -P -V -X1.5 -Y1.5 > $psfile


# plot title at bottom, adjust to your model size
pstext -J$Jstring -R -O -P -K -V -N <<END >>$psfile
10  38 14 0 1 CM  Profile $prof
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