plot profile - グリッドデータから任意の断面図を作成する 

  使うコマンド  makecpt, grdimage, grdcontour, pscoast,  psbasemap




bash スクリプト例

#  plot bathymetry profile and color index map

# parameter setting
region=130/145/25/35                    # map region east/west/south/north
proj=M10                                # map projection and scale
ticks=a5f1g5                            # boundary tick info
frame=WSne+tcolor_fill                     # boundary frame info
climit=-8000/0/1000                  # color table min/max/interval
cint=500                                # contour interval
limit=-10000/-500                       # contour min/max
aint=2000                               # annotation contour interval
grdfile=JTOPO1_30.grd                   # input bathymetry grid file
start_loc=137:00/28:30                  # profile start point
end_loc=144/30                          # profile end point
dist=0.5                                # sample interval along profile
prodatfile=profile1.xyz                 # profile data file
proregion=0/400/-5000/-1000             # profile region
proproj=x0.05/0.002                      # profile projection (cartesian) and scale
pticks=f10a50g50/f100a500g500:Depth:         # profile boundarytick info
cptfile=haxby_grad.cpt                  # color table
psfile=profile+index.ps                  # output postscript file name
#
# making color table
gmt makecpt -Chaxby -T$climit -Z > $cptfile
#
# extract values along track from grid file
gmt project -C$start_loc -E$end_loc -G$dist -Q -V > tmpfile
gmt grdtrack tmpfile -G$grdfile -L -V > $prodatfile
#
# plot profile
awk '{print $3,$4}' $prodatfile | gmt psxy -R$proregion -J$proproj -B$pticks -Wthick,red -V -K > $psfile
#
# plot index map
gmt grdimage $grdfile -R$region -J$proj -C$cptfile -K -O -Y10 -V >> $psfile
gmt grdcontour $grdfile -R$region -J$proj -C$cint -A$aint -L$limit -Wthinner -K -V -O >> $psfile
gmt pscoast -R$region -J$proj -Di -Ggray -Wthin,black -K -V -O >> $psfile
gmt psxy $prodatfile -R$region -J$proj -Wthick,red -K -O -V >> $psfile
gmt psbasemap -R$region -J$proj -B$ticks -O  -V >> $psfile
#


Tips