LoginSignup
1
0

More than 5 years have passed since last update.

@s4s

例題演習GMT（３）対数グラフ＋関数の作図

Posted at

対数グラフ＋関数の作図

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の累乗の倍($2*10^n$)と5倍（または半分）（$5*10^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
1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
0