GMTとETOPO1で地形図作成

はじめに

ETOPO1 は，National Geophysical Data Center (Colorado) により作成された，全球の１分刻みの地形モデルである．
ちょっとアラビア半島の地形を見たかったので，ETOPO1 のデータを GMT で可視化してみた．

作業環境は以下の通り．使っているプログラムは GMT (Generic Mapping Tools) と ImageMagick である．

• MacBook Pro (Retina, 13-inch, Mid 2014)
• macOS High Sierra
• GMT 5.4.3
• ImageMagick 7.0.7-21

参考サイト

• https://researchmap.jp/jo2q8ykn1-1772020/
• https://www.mk-mode.com/octopress/2015/01/06/gmt-paint-map-with-elevation-data/

まずは作例を

ETOPO1 のデータから，標高 500m ピッチで色分けし，紅海（Red Sea），サウジアラビア（Saudi Arabia），ペルシャ湾（Persian Gulf）の文字を入れた． マスターカラーテーブルは no_green，文字は白抜きとした．国境線は，黒い点線で示してある．いずれも趣味的に決めている．

ETOPO1のデータを取ってくる

ETOPO1 のサイトはこれ．ETOPO1 Ice Surfaceのgrid-registered: netCDF 形式を選択．

• https://ngdc.noaa.gov/mgg/global/global.html

データのダウンロードは以下より．ETOPO1_Ice_g_gmt4.grd.gz をダウンロードし作業したいディレクトリに展開しておく．

• https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/netcdf/

Name    Last modified   Size
Parent Directory        -
ETOPO1_Ice_g_gdal.grd.gz    2011-03-22 20:52    377M
ETOPO1_Ice_g_gmt4.grd.gz    2013-01-02 18:00    377M  <= これをダウンロードして展開
readme_etopo1_netcdf.txt    2010-09-23 18:29    1.1K

コマンドの説明

変数

src=ETOPO1_Ice_g_gmt4.grd   # 大元の入力グリッドファイル（933.5MB）
range=30/60/10/35           # 作図範囲（経度の最小・最大，緯度の最小・最大）
size=12                     # 図面サイズ（12cm）
file=saudi0                 # 作業ファイル名
fig=fig_$file.eps           # 出力画像ファイル名

grdcut

gmt grdcut  $src -R$range -G$file.grd

• $src：入力グリッドファイル
• $range：切り取る範囲
• $file.grd：出力グリッドファイル

makecpt

gmt makecpt -Cno_green -T-2000/2000/500 > $file.cpt

• no_green：マスターカラーテーブル
• -T-2000/2000/500：カラーテーブルの範囲と増分（標高-2000mから+2000mまで500mピッチ）
• $file.cpt：出力CPTファイル名

grdimage

gmt grdimage $file.grd -R$range -JM$size -C$file.cpt -P -K > $fig

• $file.grd：入力グリッドファイル
• -R$range：描画範囲
• -JM$size：プロジェクション（メルカトル）と画像の大きさ（cm）
• -C$file.cpt：使用するCPTファイル

grdcontour

gmt grdcontour $file.grd -R -J -C500 -W0 -L0/3000 -A500 -O -K >> $fig

• $file.grd：入力グリッドファイル
• -C500：等高線ピッチ
• -W0：等高線の線の太さはゼロ（描かない）．指定しないとデフォルトで描かれる
• -L0/3000：等高線を描く範囲（標高ゼロメートル以上に等高線を描く）
• -A500：等高線に描く数値のピッチ

psscale

gmt psscale -Ba1000g500f500 -C$file.cpt -D6c/-1c/6c/0.3ch -O -K >> $fig

• -Ba1000g500f500：スケールに描く体裁
• a：文字を1000ピッチで，g：グリッド線を500ピッチで，f：補助線を500ピッチで描く
• -C$file.cpt：用いるCPTファイル名
• -D6c/-1c/6c/0.3ch：スケールの描画位置，長さ，幅，水平スケールの指定
  ここでは，x方向6cm，y方向-1cmを中心に長さ6cm，幅0.3cmのスケールを描画．hの指定で水平スケールとなる．

pscoast

gmt pscoast -R -J -N1/2p,#000000,3_3:0 -O -K >> $fig

• pscoast：海岸線の描画
• -N1/2p,#000000,3_3:0：境界線の描画（N1：国境，線の太さ2pt，線の色は黒，線種は点線）

psxy

gmt psxy -R -J -W3,#0000ff -K -O << EOT >> $fig
33 27.5
37 27.5
37 30
33 30
33 27.5
EOT

• psxy：線を引く
• -W3,#0000ff：太さ3ptで青い線を指定

pstext

gmt pstext -R -J -F+f+a+j -N -K -O << EOT >> $fig
44.8 25.2 16p,Helvetica,-=0.6p,black   0 MC Saudi Arabia
37.5 23   16p,Helvetica,-=0.6p,black -53 MC Red Sea
51.5 27   16p,Helvetica,-=0.6p,black -45 MC Persian Gulf
EOT

• pstest：文字列の描画
• -F+f+a+j：座標データのあとにフォント情報，描画角度，指定座標に対する位置を指定

 x    y                 f              a  j    text
---------------------------------------------------------
46   26   16p,Helvetica,-=0.6p,black   0 MC Saudi Arabia
37.5 23   16p,Helvetica,-=0.6p,black -53 MC Red Sea
51.5 27   16p,Helvetica,-=0.6p,black -45 MC Persian Gulf

psbasemap

gmt psbasemap -R -J -Bxa5f1 -Bya5f1 -O >> $fig

• psbasemap：枠線の描画
• -Bxa5f1：x軸方向の指定（数値を5単位で，補助線を1単位で描画）
• -Bya5f1：y軸方向の指定（数値を5単位で，補助線を1単位で描画）

ImageMagick コマンド

使いやすいように，epsをpngに変換し，幅600pxにリサイズしている．

mogrify -trim -density 300 -bordercolor 'transparent' -border 10x10 -format png *.eps
convert -resize 600x fig_saudi0.png fig_s_saudi0.png

スクリプト全文

# ====================
# GMT command
# ====================
src=ETOPO1_Ice_g_gmt4.grd
range=30/60/10/35
size=12
file=saudi0
fig=fig_$file.eps
gmt grdcut $src -R$range -G$file.grd
gmt makecpt -Cno_green -T-2000/2000/500 > $file.cpt
gmt grdimage $file.grd -R$range -JM$size -C$file.cpt -P -K > $fig
gmt grdcontour $file.grd -R -J -C500 -W0 -L0/3000 -A500 -O -K >> $fig
gmt psscale -Ba1000g500f500 -C$file.cpt -D6c/-1c/6c/0.3ch -O -K >> $fig
gmt pscoast -R -J -N1/2p,#000000,3_3:0 -O -K >> $fig
gmt psxy -R -J -W3,#0000ff -K -O << EOT >> $fig
33 27.5
37 27.5
37 30
33 30
33 27.5
EOT
gmt pstext -R -J -F+f+a+j -N -K -O << EOT >> $fig
44.8 25.2 16p,Helvetica,-=0.6p,black   0 MC Saudi Arabia
37.5 23   16p,Helvetica,-=0.6p,black -53 MC Red Sea
51.5 27   16p,Helvetica,-=0.6p,black -45 MC Persian Gulf
EOT
gmt psbasemap -R -J -Bxa5f1 -Bya5f1 -O >> $fig


# ====================
# ImageMagick command
# ====================
mogrify -trim -density 300 -bordercolor 'transparent' -border 10x10 -format png *.eps
convert -resize 600x fig_saudi0.png fig_s_saudi0.png