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://researchmap.jp/jo2q8ykn1-1772020/)
+ [https://www.mk-mode.com/octopress/2015/01/06/gmt-paint-map-with-elevation-data/](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，文字は白抜きとした．国境線は，黒い点線で示してある．いずれも趣味的に決めている．

![fig_s_saudi0.png](https://qiita-image-store.s3.amazonaws.com/0/129300/3bf5163d-b684-94f6-eabd-19e71ed5d04a.png)



# ETOPO1のデータを取ってくる

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

+ [https://ngdc.noaa.gov/mgg/global/global.html](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/](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
```