対数グラフ＋関数の作図

GMTは習うより慣れろ。例題と課題でGMTの使い方を理解しよう。
今回は対数グラフを作図します。また、関数をプロットする方法も取り扱います。

例題

以下の関数をGMT(5.4.3)で図化する。

y = \exp(x) \\
y = \exp(1.1x) \\
y = \exp(1.2x) \\
y = \exp(1.3x)

GMTのスクリプト例を次に示す。

plot_exam03.sh

# !/bin/bash
# --------------------------------------------------
# plot for GMT5
# 1次元対数グラフ
# --------------------------------------------------

# ----- 前処理・設定 ------------------------------

# GMTパラメーター
gmt gmtset FONT_ANNOT_PRIMARY  9p,Helvetica # 主目盛りのフォント
gmt gmtset FONT_LABEL         10p,Helvetica # 軸タイトルのフォント

# output
ps=plot_exam03.ps               # 出力psファイル

# グラフのサイズ
mp=X                            # 地図の投影法(X:xy, M:メルカトル図法)
gx=8                            # グラフの幅[cm]
gy=5                            # グラフの高さ[cm]
# log_x=l                         # x軸を対数軸にする(*)
log_y=l                         # y軸を対数軸にする(*)
jopt=$mp${gx}c${log_x}/${gy}c${log_y} # -Jオプション

# グラフの範囲(range)（-Rオプション）
rw=0                            # x軸の最小値(west)
re=5                            # x軸の最大値(east)
rs=1                            # y軸の最小値(south)
rn=1000                         # y軸の最大値(north)
ropt=$rw/$re/$rs/$rn            # -Rオプション

# 軸・目盛りの設定
ti_x="x"                    # x軸タイトル
ti_y="y"                    # y軸タイトル
tics_x=a1f0.5               # x軸の目盛り(a=主,f=副)
grid_x=                     # x軸のグリッド(g)(*)
tics_y=a1pf3                # y軸の目盛り(a=主,f=副)
grid_y=g1                   # y軸のグリッド(g)(*)
axis=WeSn                   # 軸表示（大文字：軸＋目盛り、小文字：軸）

# 線の設定
line01=0.8p,black               # 線：太さ,色,線種
line02=0.8p,blue,-              # 線：太さ,色,線種
line03=0.8p,red,.               # 線：太さ,色,線種
line04=0.8p,darkgreen,.-        # 線：太さ,色,線種

# 凡例(pslegend)の設定
pl_D_j=TL                       # 凡例枠の配置(縦:TMB, 横:LCR)
pl_D_r=TL                       # 凡例枠の配置に合わせるポジション
pl_D_o=0.2c/0.2c                # 凡例枠のオフセット(x/y)
pl_D_s=3.4c/2.3c                # 凡例枠のサイズ(x/y)
pl_D_l=1.5                      # 凡例の行間隔（比率,default=1.1）
pl_D=j$pl_D_j+j$pl_D_r+o$pl_D_o+w$pl_D_s+l$pl_D_l # 凡例枠のまとめ
pl_F_p=+p0.8p,black             # 凡例枠の線(*)
pl_F_g=+gwhite                  # 凡例枠の塗りつぶし色(*)
pl_F_s=+s                       # 凡例枠の影(*)
pl_F_r=+r                       # 凡例枠の角を丸くする(*)
pl_F=-F${pl_F_p}${pl_F_g}${pl_F_s}${pl_F_s}${pl_F_r} # 凡例枠の装飾まとめ(*)
pl_font=9p,Helvetica            # 凡例のフォント
pl_clms=1                       # 凡例の列数(または列の相対的なサイズ)
pl_off1=0.6c                    # 凡例のオフセット(1)記号までの距離
pl_len=1.0c                     # 凡例の線の長さ
pl_off2=1.3c                    # 凡例のオフセット(2)説明文までの距離
pl_text01="y=exp(x)"            # 凡例の説明文(*)
pl_text02="y=exp(1.1x)"         # 凡例の説明文(*)
pl_text03="y=exp(1.2x)"         # 凡例の説明文(*)
pl_text04="y=exp(1.3x)"         # 凡例の説明文(*)

# 画像変換(psconvert)の設定
pc_fmt=g            # g=png, b=bmp, j=jpg, t=tiff, e=eps, f=PDF, s=SVG
pc_res=300          # 画像の解像度(DPI)
pc_mgn=0.1c         # グラフ外側のマージン
pc_aa=2             # アンチエイリアシングの設定(1,2,4)

# ----- plot ------------------------------

echo "plot $ps"

# basemap(グリッド線のみ)
gmt psbasemap -J$jopt -R$ropt -Bx$grid_x -By$grid_y -B+n -P -K > $ps

# psxy(line)
awk "BEGIN{ N=100;
  for (i=0; i<=N; i++) {
    x = $rw + ($re-($rw))/N*i;
    print x, exp(x)
  } }" \
    | gmt psxy -J -R -W$line01 -K -O >> $ps

# psxy(line)
awk "BEGIN{ N=100;
  for (i=0; i<=N; i++) {
    x = $rw + ($re-($rw))/N*i;
    print x, exp(1.1*x)
  } }" \
    | gmt psxy -J -R -W$line02 -K -O >> $ps

# psxy(line)
awk "BEGIN{ N=100;
  for (i=0; i<=N; i++) {
    x = $rw + ($re-($rw))/N*i;
    print x, exp(1.2*x)
  } }" \
    | gmt psxy -J -R -W$line03 -K -O >> $ps

# psxy(line)
awk "BEGIN{ N=100;
  for (i=0; i<=N; i++) {
    x = $rw + ($re-($rw))/N*i;
    print x, exp(1.3*x)
  } }" \
    | gmt psxy -J -R -W$line04 -K -O >> $ps

# 凡例(line)
gmt pslegend -R -J -D$pl_D $pl_F --FONT_ANNOT_PRIMARY=$pl_font -K -O << EOF >> $ps
# ----------------------------------------------------------------------
# N 列数(または列の相対的なサイズ)
  N $pl_clms
# ----------------------------------------------------------------------
# S offset1  symbol size    fill  line    offset2  text
  S $pl_off1 -      $pl_len -     $line01 $pl_off2 $pl_text01
  S $pl_off1 -      $pl_len -     $line02 $pl_off2 $pl_text02
  S $pl_off1 -      $pl_len -     $line03 $pl_off2 $pl_text03
  S $pl_off1 -      $pl_len -     $line04 $pl_off2 $pl_text04
EOF

# basemap(枠線)
gmt psbasemap -J -R -Bx$tics_x+l"$ti_x" -By$tics_y+l"$ti_y" -B$axis -O >> $ps

# ----- 後処理 ------------------------------

# png変換
gmt psconvert -A$pc_mgn -E$pc_res -T$pc_fmt -Qt$pc_aa -Qg$pc_aa $ps

作成されたグラフ。

解説

対数軸にする方法

軸を対数軸にするには、-Jオプションでサイズの次に「l」を記述する。
例では、「-JX8c/5cl」であり、y軸を対数に指定している。

対数軸の目盛りの設定

対数軸での目盛り(annotation, tics, grid)は、レベル(1〜3)を指定する。
- レベル1：1, 10, 100...という風に、10の累乗($10^n$)の位置。
- レベル2：1, 2, 5, 10, 20, 50, 100...という風に、10の累乗の倍($210^n$)と5倍（または半分）（$510^n$）の位置。概ねレベル1を3等分した位置になる。
- レベル3：1,2,3,4,5,6,7,8,9,10,20...という風に、$m*10^n~~(m=1..9)$の位置。
さらに、主目盛り(a)の後ろ添字に「p」を付けると、目盛り数値が累乗表記（$10^n$）になる。
例ではy軸を、「-Byg1」「-Bya1pf3」としている。これは、主目盛り(a)はレベル1で累乗表記、副目盛り(f)はレベル3、グリッド(g)はレベル1を意味する。

関数のプロット方法

GMTには関数をプロットするコマンド・オプションはない！
したがって、関数をプロットする方法には、以下のように行う。
1. プロットする関数のx,yの値をファイルで作成しておき、それを読み込む。
2. プロットする関数のx,yの値を、psxyの標準入力として受け取る。
例では、awkを用いて、関数のデータを出力し、それをパイプとしてpsxyの標準入力に繋いでいる（上記の方法2.になる）。

awkによる関数作成

例では、x軸の範囲を100分割して関数結果(x,y)を出力している。
awkは、ファイルを読み込みながら、各行に置換等の処理を行ったりする言語だが、"BEGIN{}"だけだと入力ファイルが不要となり、出力のみとなる。
awkスクリプトの詳細は割愛するが、awkはGMTと非常に相性が良く、今回のような関数のプロットに適している。

点線・破線のプロット方法

psxy(line)で、点線や破線、一点鎖線をプロットするには、-Wオプションで指定する。
具体的には、「-W線の太さ,色,線種」とする。線種の書き方は以下の通り（実線は線種の指定なし）。
- 破線：「-」
- 点線：「.」
- 一点鎖線：「.-」
- 二点鎖線：「..-」
破線の間隔などを細かく設定することもできるが、大抵は上記の4種類で十分であろう。

課題

課題(1)

以下のようなグラフを作成せよ。
変更点：y軸の目盛り数値・グリッド。x軸の範囲・目盛り・グリッド。

課題(2)

以下のようなグラフを作成せよ。
変更点：関数、x軸・y軸全般。
ヒント：awkでの、関数$\ln(x)$は、"log(x)"となる。

解答

課題(1)

目盛り数値は「p」を外すと、累乗表記でなくなる。

plot_assign03-01.sh（変更点）

grid_x=g1                   # x軸のグリッド(g)(*)
tics_y=a2f3                 # y軸の目盛り(a=主,f=副)
grid_y=g3                   # y軸のグリッド(g)(*)

課題(2)

plot_assign03-02.sh

# !/bin/bash
# --------------------------------------------------
# plot for GMT5
# 1次元対数グラフ
# --------------------------------------------------

# ----- 前処理・設定 ------------------------------

# GMTパラメーター
gmt gmtset FONT_ANNOT_PRIMARY  9p,Helvetica # 主目盛りのフォント
gmt gmtset FONT_LABEL         10p,Helvetica # 軸タイトルのフォント

# output
ps=plot_assign03-02.ps          # 出力psファイル

# グラフのサイズ
mp=X                            # 地図の投影法(X:xy, M:メルカトル図法)
gx=8                            # グラフの幅[cm]
gy=5                            # グラフの高さ[cm]
log_x=l                         # x軸を対数軸にする(*)
# log_y=l                         # y軸を対数軸にする(*)
jopt=$mp${gx}c${log_x}/${gy}c${log_y} # -Jオプション

# グラフの範囲(range)（-Rオプション）
rw=1                            # x軸の最小値(west)
re=1000                         # x軸の最大値(east)
rs=0                            # y軸の最小値(south)
rn=10                        　 # y軸の最大値(north)
ropt=$rw/$re/$rs/$rn            # -Rオプション

# 軸・目盛りの設定
ti_x="x"                    # x軸タイトル
ti_y="y"                    # y軸タイトル
tics_x=a1pf3                # x軸の目盛り(a=主,f=副)
grid_x=g3                   # x軸のグリッド(g)(*)
tics_y=a2f1                 # y軸の目盛り(a=主,f=副)
grid_y=g2                   # y軸のグリッド(g)(*)
axis=WeSn                   # 軸表示（大文字：軸＋目盛り、小文字：軸）

# 線の設定
line01=0.8p,black               # 線：太さ,色,線種
line02=0.8p,blue,-              # 線：太さ,色,線種
line03=0.8p,red,.               # 線：太さ,色,線種
line04=0.8p,darkgreen,.-        # 線：太さ,色,線種

# 凡例(pslegend)の設定
pl_D_j=TL                       # 凡例枠の配置(縦:TMB, 横:LCR)
pl_D_r=TL                       # 凡例枠の配置に合わせるポジション
pl_D_o=0.2c/0.2c                # 凡例枠のオフセット(x/y)
pl_D_s=3.4c/2.3c                # 凡例枠のサイズ(x/y)
pl_D_l=1.5                      # 凡例の行間隔（比率,default=1.1）
pl_D=j$pl_D_j+j$pl_D_r+o$pl_D_o+w$pl_D_s+l$pl_D_l # 凡例枠のまとめ
pl_F_p=+p0.8p,black             # 凡例枠の線(*)
pl_F_g=+gwhite                  # 凡例枠の塗りつぶし色(*)
pl_F_s=+s                       # 凡例枠の影(*)
pl_F_r=+r                       # 凡例枠の角を丸くする(*)
pl_F=-F${pl_F_p}${pl_F_g}${pl_F_s}${pl_F_s}${pl_F_r} # 凡例枠の装飾まとめ(*)
pl_font=9p,Helvetica            # 凡例のフォント
pl_clms=1                       # 凡例の列数(または列の相対的なサイズ)
pl_off1=0.6c                    # 凡例のオフセット(1)記号までの距離
pl_len=1.0c                     # 凡例の線の長さ
pl_off2=1.3c                    # 凡例のオフセット(2)説明文までの距離
pl_text01="y=ln(x)"             # 凡例の説明文(*)
pl_text02="y=1.1*ln(x)"         # 凡例の説明文(*)
pl_text03="y=1.2*ln(x)"         # 凡例の説明文(*)
pl_text04="y=1.3*ln(x)"         # 凡例の説明文(*)

# 画像変換(psconvert)の設定
pc_fmt=g            # g=png, b=bmp, j=jpg, t=tiff, e=eps, f=PDF, s=SVG
pc_res=300          # 画像の解像度(DPI)
pc_mgn=0.1c         # グラフ外側のマージン
pc_aa=2             # アンチエイリアシングの設定(1,2,4)

# ----- plot ------------------------------

echo "plot $ps"

# basemap(グリッド線のみ)
gmt psbasemap -J$jopt -R$ropt -Bx$grid_x -By$grid_y -B+n -P -K > $ps

# psxy(line)
awk "BEGIN{ N=100;
  for (i=0; i<=N; i++) {
    x = $rw + ($re-($rw))/N*i;
    print x, log(x)
  } }" \
    | gmt psxy -J -R -W$line01 -K -O >> $ps

# psxy(line)
awk "BEGIN{ N=100;
  for (i=0; i<=N; i++) {
    x = $rw + ($re-($rw))/N*i;
    print x, 1.1*log(x)
  } }" \
    | gmt psxy -J -R -W$line02 -K -O >> $ps

# psxy(line)
awk "BEGIN{ N=100;
  for (i=0; i<=N; i++) {
    x = $rw + ($re-($rw))/N*i;
    print x, 1.2*log(x)
  } }" \
    | gmt psxy -J -R -W$line03 -K -O >> $ps

# psxy(line)
awk "BEGIN{ N=100;
  for (i=0; i<=N; i++) {
    x = $rw + ($re-($rw))/N*i;
    print x, 1.3*log(x)
  } }" \
    | gmt psxy -J -R -W$line04 -K -O >> $ps

# 凡例(line)
gmt pslegend -R -J -D$pl_D $pl_F --FONT_ANNOT_PRIMARY=$pl_font -K -O << EOF >> $ps
# ----------------------------------------------------------------------
# N 列数(または列の相対的なサイズ)
  N $pl_clms
# ----------------------------------------------------------------------
# S offset1  symbol size    fill  line    offset2  text
  S $pl_off1 -      $pl_len -     $line01 $pl_off2 $pl_text01
  S $pl_off1 -      $pl_len -     $line02 $pl_off2 $pl_text02
  S $pl_off1 -      $pl_len -     $line03 $pl_off2 $pl_text03
  S $pl_off1 -      $pl_len -     $line04 $pl_off2 $pl_text04
EOF

# basemap(枠線)
gmt psbasemap -J -R -Bx$tics_x+l"$ti_x" -By$tics_y+l"$ti_y" -B$axis -O >> $ps

# ----- 後処理 ------------------------------

# png変換
gmt psconvert -A$pc_mgn -E$pc_res -T$pc_fmt -Qt$pc_aa -Qg$pc_aa $ps