plane.gmt

这是一个TOMO TO GMT的源码

#!/bin/csh
#
# to plot values of 
#   1:  KHIT matrix 
#   2:  DWS
#   3:  diag. reso element (RDE)
#   4:  spread values
# in horizontal depth sections 
#
# 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!
#
#  call by: plane.gmt <1,2,3,4>  
#					(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

# !! 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/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=X3

# 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/plane.khit.ps
   set clabel="KHIT"
   set pos=10
   breaksw
  case  2:
   set cpt=$wpath/grey_dws.cpt
   set psfile=$wpath/plane.dws.ps
   set clabel="Derivative Weigthed Sum (DWS)"
   set pos=11
   breaksw
  case  3:
   set cpt=$wpath/grey1.cpt
   set psfile=$wpath/plane.rde.ps
   set clabel="Resolution Diag. Element"
   set pos=12
   breaksw
  case  4:
   set cpt=$wpath/sprd2.cpt
   set psfile=$wpath/plane.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/plane.khit.sml.ps
   set clabel="no. of rays (KHIT)"
   set pos=10
   breaksw
  case  2:
   set cpt=$wpath/grey_dws.cpt
   set psfile=$wpath/plane.dws.sml.ps
   set clabel="Derivative Weigthed Sum (DWS)"
   set pos=11
   breaksw
  case  3:
   set cpt=$wpath/grey1.cpt
   set psfile=$wpath/plane.rde.sml.ps
   set clabel="Reso. Diag. Element"
   set pos=12
   breaksw
  case  4:
   set cpt=$wpath/sprd2.cpt
   set psfile=$wpath/plane.sprd.sml.ps
   set clabel="spread value"
   set pos=13
 endsw
endif

### now start GMT job

#plot colorscale
psscale -C$cpt -D4.5/0.5/5/0.25h -L -B:"$clabel": -K -P -X-0.25 -Y0.5 > $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 infile1=$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 'NR>2 {print $3,$4,$pos}' pos=$pos $infile1 |\
xyz2grd -I$Istring -N0 -Gsprd.grd -R$Rstring -V

#image grid
grdimage sprd.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_invblo.xyz
 nawk 'BEGIN {print ">"} { if (NF==4) print $2,$1; if (NF==1) print ">" }' $nodefile |\
 psxy -J$Jstring -R -W0.5/200 -M -O -K -P >>$psfile
else
 set nodefile=$inpath/plane_nodes.xyz
 nawk '{print $4,$3}' $nodefile | \
 psxy  -J$Jstring -R$Rstring -Sx0.05 -W0.5/200 -O -K -V -P >>$psfile 
endif


# 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 sprd.grd

exit