GMTでグラフを描く

# 1. この記事の概要

+ GMT (Generic Mapping Tools)でグラフを描画するスクリプトを紹介します．
+ 紹介するスクリプトは Mac のターミナルで使用することを前提としています．
+ GMT コマンドは GMT5 で使用することを前提としています．GMT5 がインストールされていなければ動きません．
+ 紹介するスクリプトは，片対数グラフ作成，２軸を利用したグラフの作成です．
+ おまけとして，ImageMagickによるepsからpngにフォーマット変換するコマンドを載せておきます．



# 2. 片対数グラフ

## データの準備
プロットするデータを作成．これは適当な数値を打ち込みます．
この事例では，\_inp\_a.txt と \_inp\_b.txt の２つのファイルを準備します．

```
cat << EOT > _inp_a.txt
0.01 -6
0.1  -4
  1   0
 20   8
 50   4
100   3
EOT

cat << EOT > _inp_b.txt
0.01  3
0.1   4
  1   5
  7   2
 50  -4
100  -8
EOT
```

## GMTコマンドの作成
スクリプト構成の大きな流れは以下のとおりです．使い回すことを意識して，シェルスクリプトの変数を使っています．

+ 作図範囲，グラフの大きさ，軸ラベル，出力ファイル名などを定義
+ gmt set で共通のフォーマットを定義
+ gmt psbasemap でグラフの軸を作成
+ gmt psxy でプロット実行

スクリプト全文は以下のとおりです．

```
range=0.01/100/-10/10
scale=5l/3
xga=a1g3
yga=g2a4
xlabel='x-label'
ylabel='y-label'
fig=fig_gra_slog.eps

gmt set FONT_ANNOT_PRIMARY 12
gmt set FONT_LABEL 8
gmt set MAP_TICK_LENGTH_PRIMARY 0c

gmt psbasemap -R$range -JX$scale -Bx$xga+l"$xlabel" -By$yga+l"$ylabel" -BWSen -P -K > $fig
gmt psxy _inp_a.txt -R -J -W1.5 -K -O >> $fig
gmt psxy _inp_a.txt -R -J -SC0.2 -G#000000 -N -K -O >> $fig
gmt psxy _inp_b.txt -R -J -W1.5,3_3:0 -K -O >> $fig
gmt psxy _inp_b.txt -R -J -SC0.2 -W1,#0000ff -G#ffffff -N -K -O >> $fig
echo '0.01 -10' | gmt psxy -R -J -Sp -O >> $fig
```



### 用いている変数の意味

+ **range=0.01/100/-10/10**：横軸範囲を0.01から100に，縦軸範囲を-10から10にセットします．
+ **scale=5l/3**：描画寸法を横5cm，縦3cmにセットします．なお横軸は対数にしたいので「l」（エル）をつけています．
+ **xga=a1g3**：横軸に対し1ピッチで目盛りラベルを，3ピッチでグリッド線を表示するよう指定，という意味ですが，対数軸の場合の扱いは普通軸と異なっています．詳細はGMTの解説書を確認願います．
+ **yga=g2a4**：縦軸に対し2ピッチでグリッド線を，4ピッチで目盛りラベルを表示するよう指定します．
+ **xlabel='x-label'**：横軸ラベルをセット
+ **ylabel='y-label'**：縦軸ラベルをセット
+ **fig=fig\_gra\_slog.eps**：出力ファイル名をセット．epsもしくはpsです．

### プロット方法について少し解説
#### **-K** と **-O** について
コマンドの中に，**-K** とか **-O** という記載がありますが，これはGMTのスクリプトを描く上で大変重要な意味を持ちます．

+ **-K** は複数行の命令を描く上で，最初の行に必要なオプションで，次の行にコマンドが続くことを示します．
+ **-O** はこの行の前にもコマンドがあることを示すオプションです．よって複数行の命令の最後の行には必ず必要です．
+ 途中行については**-K -O**をつけます．この後ろにもコマンド行があり，前にもコマンド行があるよ，という意味でしょうか．
+ 以下の行は，**-O** を含む命令を行うことのみを目的としたダミーです．最後の行には**-O**をつけなければならないのですが，色々挿入していて，**-K** や **-O** を書き換えていくのが面倒なので，グラフの迷惑にならないよう，左下に点を打つようにしています．

```
echo '0.01 -10' | gmt psxy -R -J -Sp -O >> $fig
```

#### psxy によるプロットの実行
psxy というコマンドにより，\_inp\_a.txt と \_inp\_b.txtの内容をプロットしています．

##### \_inp\_a.txtに対して
プロットを黒丸で，プロット同士を実線で結ぶことを指定しています．ただし，コマンドは上から順に実行されるので，先に線を引き，丸印を描くよう命令します．

```
gmt psxy _inp_a.txt -R -J -W1.5 -K -O >> $fig
gmt psxy _inp_a.txt -R -J -SC0.2 -G#000000 -N -K -O >> $fig
```

##### \_inp\_b.txtに対して
プロットを丸で，丸の縁を青線で，中を白で塗りつぶし，プロット同士を破線で結ぶことを指定しています．ただし，コマンドは上から順に実行されるので，先に線を引き，丸印を描くよう命令します．

```
gmt psxy _inp_b.txt -R -J -W1.5,3_3:0 -K -O >> $fig
gmt psxy _inp_b.txt -R -J -SC0.2 -W1,#0000ff -G#ffffff -N -K -O >> $fig
```

### 成果品
<img src="https://qiita-image-store.s3.amazonaws.com/0/129300/2c0447fc-14fa-9902-4799-57ea2828d62a.png" width=55%>


# 3. 2軸の利用
## データの準備

```
cat << EOT > _inp_a.txt
 2.5  1
 7.5  5
12.5  8
17.5  3
EOT

cat << EOT > _inp_b.txt
-75  40
-25  20
 25 -30
 70  30
EOT
```

## GMT コマンドの作成
スクリプト構成は，基本的には片対数の場合と同じですが，2軸を利用するための定義が加わっています．

+ 作図範囲を示す変数rangeでは，rangeSWとrangeNEの２つを定義しています．これは地図でいうSW，下横軸と左縦軸に対しプロットを行うものです．対してNEは上横軸と右縦軸に対しプロットを行うものです．
+ scaleはグラフの大きさを示すものであり，SWとNEに対し共通です．ここでは横5cm，縦3cmのグラフ寸法を指定します．
+ 変数 afgS, afgW, afgN, afgE はそれぞれの軸でのラベル表示間隔・軸線表示間隔・グリッド表示間隔を指定します．
+ 軸ラベルについても東西南北で指定できます．
+ GMTは基本的に重ね書きにより複雑な描画を行うようにできているので，この考え方に従い，プロット用コマンドを書きます．
+ ここではgmt psbasemapと gmt psxyにより，SW軸に対するプロットを行います．
+ 続いてNE軸に対するプロットを行っています．

```
rangeSW=0/20/0/10
rangeNE=-100/100/-50/50
scale=5/3
afgS=a5f1g5
afgW=a2f1g2
afgN=a50f10
afgE=a20f10
labelS='Label-S'
labelW='Label-W'
labelN='Label-N'
labelE='Label-E'
fig=fig_gra_2axes.eps

gmt set FONT_ANNOT_PRIMARY 8
gmt set FONT_LABEL 12
gmt set MAP_TICK_LENGTH_PRIMARY 0.2c

gmt psbasemap -R$rangeSW -JX$scale -Bx$afgS+l"$labelS" -By$afgW+l"$labelW" -BWS -P -K > $fig
gmt psxy _inp_a.txt -R -J -W2 -K -O >> $fig
gmt pstext -R -J -F+f+a+j -N -K -O << EOT >> $fig
17.5 3 10p 0 LM SW
EOT

gmt psbasemap -R$rangeNE -JX$scale -Bx$afgN+l"$labelN" -By$afgE+l"$labelE" -BNE -K -O >> $fig
gmt psxy _inp_b.txt -R -J -W2,5_5:0 -K -O >> $fig
gmt pstext -R -J -F+f+a+j -N -K -O << EOT >> $fig
70 30 10p 0 LM NE
EOT

echo '-100 -50' | gmt psxy -R -J -Sp -O >> $fig
```

### 成果品
<img src="https://qiita-image-store.s3.amazonaws.com/0/129300/4fa9c763-bfeb-c579-7e30-0110475d485b.png" width=55%>


# 4. eps から png への画像変換(ImageMagick利用)
下のコマンドは，GMTで作成したeps画像をpngに変換するImageMagickのコマンドです．

そのディレクトリにある全てのepsファイルを，余白を調整してpng画像に変換するものです．

２行目のコメント行にしてあるものは，psファイルをpngに変換するものです．psではpngに変換する際画像の向きをあわせるため，90度回転させた処理を行うようにしています．

```
mogrify -trim -density 300 -bordercolor 'transparent' -border 10x10 -format png *.eps
#mogrify -rotate 90 -trim -density 300 -bordercolor 'transparent' -border 10x10 -format png *.ps
```

この処理を行うことにより，GMTで作った綺麗な画像を，wordでもpowerpointでも，TeXでも，好きなものに貼り付けできるようになります．

以上