GMT
ETOPO1

GMTとETOPO1で地形図作成

はじめに

ETOPO1 は,National Geophysical Data Center (Colorado) により作成された,全球の1分刻みの地形モデルである.
ちょっとアラビア半島の地形を見たかったので,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

参考サイト

まずは作例を

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

fig_s_saudi0.png

ETOPO1のデータを取ってくる

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

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

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