GMT5で地図にスケールをつける

【更新履歴】2020.2.2　スクリプト内のコメントを修正

公式サイトの情報を参考に、GMT5でスケール付きの地図を作成してみたので紹介します。地形データにはETOPO1を使用しました。

#!/bin/sh

# parameter setting
# map region west/east/south/north
region=128/140/28/36
# map projection and scale 
proj=M18
# frame
frame=WeSn
# xaxis (f: minor tick spacing)
xaxis=x4f2
# yaxis (f: minor tick spacing)
yaxis=y4f2
# master color table
cpt=seafloor
# color table min/max
cptlimit=-6000/0
# input bathymetry grid file name                                   
grdfile=ETOPO1_Bed_g_gmt4.grd
# tmp bathymetry grid file name     
tmpgrdfile=ETOP1_tmp.grd
# output color table file name                 
cptfile=my_seafloor.cpt
# output postscript file name 
psfile=map.ps 

# set default parameter
gmt set FORMAT_GEO_MAP dddF

# cut grid data
gmt grdcut $grdfile -R$region -G$tmpgrdfile

# make color table
gmt makecpt -C$cpt -T$cptlimit > $cptfile

# plot map
gmt grdimage $tmpgrdfile -R$region -J$proj -C$cptfile -K -V > $psfile
gmt pscoast -R -J -Da -Ggray -Wthinnest,black -K -O -V >> $psfile

# plot frame
gmt psbasemap -R -J -B$frame -B$xaxis -B$yaxis -K -O -V >> $psfile

# plot scale
gmt psbasemap -R -J -Lg136.4/28.5+c32N+w300k+f+l+jBL --MAP_SCALE_HEIGHT=10p --MAP_TICK_PEN_PRIMARY=white --FONT_ANNOT_PRIMARY=white --FONT_LABEL=white -O -V >> $psfile

# convert to raseter file
gmt psconvert $psfile -A -Tj -V

# remove tmpfile
rm -f $tmpgrdfile

私の理解が正しければ、"-Lg136.4/28.5+c32N+w300k+f+l+jBL"が意味するところは、以下のようになっています。
・北緯32度の距離での（+c32N）
・300 kmの長さのスケールバーを（+w300k）
・左下が（+jBL）
・緯度28.5経度136.4となるように（-Lg136.4/28.5）
・fancyモードで（+f）
・距離の単位（ここでは"km"）をバーの上につけて（+l）描く

デフォルトだと黒で描かれて見にくいので、以下のパラメータを一時的に白にすることで見やすくしています。
・MAP_TICK_PEN_PRIMARY
・FONT_ANNOT_PRIMARY
・FONT_LABEL
また、デフォルトだとスケールバーの厚みが小さい印象があったので、MAP_SCALE_HEIGHTを変更しています。