LoginSignup
4
3

More than 5 years have passed since last update.

@s4s

例題演習GMT（４）時系列グラフ

Posted at

時系列データの作図

GMTは習うより慣れろ。例題と課題でGMTの使い方を理解しよう。
今回は時系列データを作図します。

例題

以下のようなグラフをGMT(5.4.3)で作成する。

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

plot_exam04.sh
#!/bin/bash
#--------------------------------------------------
# plot for GMT5
# 1次元時系列グラフ
#--------------------------------------------------

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

# GMTパラメーター
gmt gmtset FONT_ANNOT_PRIMARY=9p,Helvetica   # 主目盛りのフォント
gmt gmtset FONT_LABEL=10p,Helvetica          # 軸タイトルのフォント
# 日付関係(=前後にスペースを入れない)
gmt gmtset FORMAT_DATE_MAP=mm/dd # 日付目盛り(yy[yy]=年,mm=月,dd=日,o=月名,jjj=日数)
gmt gmtset FORMAT_TIME_PRIMARY_MAP=a # 主目盛りの月・週表示：[f|a|c]

# input
in01=chiba_wea.txt              # 入力ファイル
in01_h=1                        # ヘッダー行数
in01_c=0,1                      # 入力する列(0始まり）

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

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

# グラフの範囲(range)（-Rオプション）
rw=2017-01-01T00:00             # x軸の最小値(west)
re=2017-01-05T00:00             # x軸の最大値(east)
rs=980                          # y軸の最小値(south)
rn=1030                         # y軸の最大値(north)

# 軸・目盛りの設定
ti_x="Date (2017)"          # x軸タイトル
ti_y="Pressure [hPa]"       # y軸タイトル
tics_x=a1Df6h               # x軸の目盛り(a=主,f=副)
grid_x=                     # x軸のグリッド(g)(*)
tics_y=a10f5                # y軸の目盛り(a=主,f=副)
grid_y=g10                  # y軸のグリッド(g)(*)
axis=WeSn                   # 軸表示（大文字：軸＋目盛り、小文字：軸）

# 線の設定
line01=0.8p,black               # 線：太さ,色,線種


#----- plot ------------------------------
echo "plot $ps"

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

# psxy(line)
gmt psxy $in01 -h$in01_h -i$in01_c -J -R -W$line01 -K -O >> $ps

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


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

# 画像変換(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)

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

# 中間ファイルの削除
rm -f gmt.history gmt.conf

使用したデータファイル。
なお、データは気象庁HPのデータを加工して作成した（実際の観測値を保証するものではありません）。

chiba_wea.txt
Date    pres[hPa]
2017-01-01T00:00    1021.8
2017-01-01T01:00    1021.7
2017-01-01T02:00    1021.8
2017-01-01T03:00    1021.7
2017-01-01T04:00    1021.6
2017-01-01T05:00    1021.7
2017-01-01T06:00    1021.8
2017-01-01T07:00    1022.4
2017-01-01T08:00    1022.9
2017-01-01T09:00    1023.4
2017-01-01T10:00    1023.4
2017-01-01T11:00    1022.7
2017-01-01T12:00    1021.7
2017-01-01T13:00    1021.0
2017-01-01T14:00    1020.5
2017-01-01T15:00    1021.1
2017-01-01T16:00    1021.7
2017-01-01T17:00    1022.2
2017-01-01T18:00    1022.5
2017-01-01T19:00    1022.9
2017-01-01T20:00    1022.9
2017-01-01T21:00    1022.6
2017-01-01T22:00    1022.3
2017-01-01T23:00    1022.3
2017-01-02T00:00    1021.8
2017-01-02T01:00    1021.5
2017-01-02T02:00    1020.8
2017-01-02T03:00    1020.2
2017-01-02T04:00    1020.1
2017-01-02T05:00    1020.1
2017-01-02T06:00    1020.1
2017-01-02T07:00    1019.4
2017-01-02T08:00    1020.1
2017-01-02T09:00    1019.8
2017-01-02T10:00    1019.0
2017-01-02T11:00    1017.5
2017-01-02T12:00    1016.1
2017-01-02T13:00    1015.2
2017-01-02T14:00    1014.3
2017-01-02T15:00    1013.8
2017-01-02T16:00    1013.6
2017-01-02T17:00    1013.9
2017-01-02T18:00    1013.7
2017-01-02T19:00    1013.7
2017-01-02T20:00    1013.5
2017-01-02T21:00    1012.7
2017-01-02T22:00    1012.3
2017-01-02T23:00    1012.1
2017-01-03T00:00    1011.7
2017-01-03T01:00    1011.5
2017-01-03T02:00    1011.5
2017-01-03T03:00    1011.3
2017-01-03T04:00    1011.2
2017-01-03T05:00    1011.2
2017-01-03T06:00    1011.8
2017-01-03T07:00    1013.0
2017-01-03T08:00    1013.6
2017-01-03T09:00    1014.4
2017-01-03T10:00    1014.6
2017-01-03T11:00    1013.7
2017-01-03T12:00    1012.7
2017-01-03T13:00    1012.3
2017-01-03T14:00    1012.8
2017-01-03T15:00    1013.2
2017-01-03T16:00    1014.0
2017-01-03T17:00    1014.5
2017-01-03T18:00    1015.4
2017-01-03T19:00    1016.0
2017-01-03T20:00    1016.5
2017-01-03T21:00    1016.8
2017-01-03T22:00    1016.8
2017-01-03T23:00    1016.8
2017-01-04T00:00    1016.6
2017-01-04T01:00    1016.6
2017-01-04T02:00    1016.3
2017-01-04T03:00    1016.0
2017-01-04T04:00    1015.5
2017-01-04T05:00    1015.3
2017-01-04T06:00    1015.1
2017-01-04T07:00    1014.7
2017-01-04T08:00    1014.7
2017-01-04T09:00    1014.6
2017-01-04T10:00    1013.8
2017-01-04T11:00    1011.9
2017-01-04T12:00    1010.6
2017-01-04T13:00    1009.7
2017-01-04T14:00    1009.8
2017-01-04T15:00    1010.2
2017-01-04T16:00    1010.6
2017-01-04T17:00    1011.6
2017-01-04T18:00    1012.5
2017-01-04T19:00    1013.1
2017-01-04T20:00    1013.5
2017-01-04T21:00    1013.4
2017-01-04T22:00    1013.8
2017-01-04T23:00    1014.1
2017-01-05T00:00    1013.9

解説

時系列データのフォーマット

  • 時系列データは、上記のデータのように、「YYYY-MM-DDThh: mm: ss」とする（日付と時刻の間に「T」を挟む）（絵文字に変換されるので、「:」のあとにスペースを挿入しているが、実際には不要）。
  • この書式は、ISO8601形式の拡張形式と呼ばれる。
  • 使用しない時間部分は省略できる。例えば、上記のデータでは秒の部分（:ss）を省略している（0秒とみなされる）。

目盛りの設定

  • 時系列データの目盛りは、(1)-Bオプション（psbasemap）と、(2)GMTパラメータ（FORMAT_DATE_MAP、FORMAT_TIME_PRIMARY_MAP）で設定する。

-Bオプションの日付設定

  • 軸の目盛り(annotation, tics, grid)の設定方法は、-Ba[目盛り数値の間隔]t[目盛り線の間隔]g[グリッド線の間隔]、という風に設定する（普通のデータと同様）。
  • 時系列データの場合、[間隔]の設定には時間の単位も付ける。例えば、「-Bf1d」では、「目盛り線(tics)が1日(day)間隔」という意味になる。
  • 間隔の単位は、y=年、o=月、d=日、h=時、m=分、s=秒、u=週（月曜はじまり）となる。月はmではなくoなので注意。
  • 間隔の単位(y,m,d等)は大文字と小文字で、annotation（目盛り数値）の表記が異なることに注意！
  • 間隔の単位を小文字にすると、annotationは数値で表記される。

  • 間隔の単位を大文字にすると、annotationの表記は、パラメータのFORMAT_DATE_MAPの設定値になる。

  • 例えば、-Bx1D6h（annotation部分を大文字）とした場合、x軸は次のようになる。

  • これを、-Bx1d6h（annotation部分を小文字）とした場合は次のようになる。

  • annotation部分を小文字にした場合は、FORMAT_DATE_MAPの設定に関わらず、目盛り数値が「01, 02」と日数表記になっている。

日付表示の設定（FORMAT_DATE_MAP）

  • annotationなどの間隔の設定は、-Bオプションで設定するが、annotationの表記は、FORMAT_DATE_MAPで設定する。
  • 表記のフォーマットは、yyyy=年（4桁）、yy=年（下2桁）、mm=月（数値）、o=月（月名）、dd=日、jjj=年初からの日数、を使って指定する。（※）週の表記はww=週数なのだが、自分の環境では表示できなかった（GMTのバグと思われる）。
  • 例えば、「FORMAT_DATE_MAP=mm/dd」とすれば、1月10日は"01/10"と表記される（-Bオプションの間隔単位が大文字であること）。
  • フォーマットの最初に、-（ハイフン）を付けると、2桁目の0がなくなる。例えば、1月10日を、「-mm/dd」とすれば"1/10"となる。
  • 月を数値ではなく、月名（January, Februaryなど）で表示する場合は、パラメータのFORMAT_TIME_PRIMARY_MAPを設定する必要がある。

月名表示の設定（FORMAT_TIME_PRIMARY_MAP）

  • 月名表示を設定するパラメータFORMAT_TIME_PRIMARY_MAPの設定値は、3種類ある。（※）週名も同様に設定できるらしいが、うまくいかなかったので、省略。
  • f(=full)：「January」のように、省略なしで表示する。
  • a(=abbreviated)：「Jan」のように、3文字の短縮形で表示する。
  • c(=character)：「J」のように、1文字の頭文字で表示する。

  • f, aを大文字で指定した場合は、単語が全部大文字になる。Fの場合は「JANUARY」のようになる。

  • 例えば、FORMAT_DATE_MAP=-o.dd、FORMAT_TIME_PRIMARY_MAP=a、とした場合は、次のようになる。

時刻表示の設定（FORMAT_CLOCK_MAP）

  • 日付の表示はFORMAT_DATE_MAPだが、時刻の表示はFORMAT_CLOCK_MAPで設定する。
  • 表記のフォーマットは、hh=時、mm=分、ss=秒、を使って指定する。
  • 例えば、「FORMAT_CLOCK_MAP=hh:mm」とすれば、2時15分は"02:15"となる。

  • なお、日付と時刻を同じ軸で同時に表示することはできない（"MM/DD hh:mm"のようなフォーマットは無理）。どちらかを副軸を利用するなど、工夫が必要となる。

課題

課題(1)

  • 以下に示すtemp.txtのデータを使って、下記のような気温のグラフを作成せよ。
  • データはtab区切りで、日付(YYYY/MM/DD hh:mm)と気温[℃]のデータとなっている。
  • 日付はGMTの形式（ISO8601形式）と異なるので、変換が必要。
  • 温度記号「℃」は、"@~\260@~C"で出力できる。
  • なお、データは気象庁HPのデータを加工して作成した（実際の観測値を保証するものではありません）。

temp.txt
Date-time   Temp[℃]
2017/08/01 00:00    27.3
2017/08/01 01:00    27.2
2017/08/01 02:00    27.2
2017/08/01 03:00    26.8
2017/08/01 04:00    26.9
2017/08/01 05:00    26.8
2017/08/01 06:00    27.0
2017/08/01 07:00    27.6
2017/08/01 08:00    27.6
2017/08/01 09:00    28.4
2017/08/01 10:00    29.1
2017/08/01 11:00    28.9
2017/08/01 12:00    29.8
2017/08/01 13:00    30.6
2017/08/01 14:00    29.9
2017/08/01 15:00    27.7
2017/08/01 16:00    25.6
2017/08/01 17:00    25.1
2017/08/01 18:00    24.6
2017/08/01 19:00    24.2
2017/08/01 20:00    24.2
2017/08/01 21:00    24.2
2017/08/01 22:00    24.1
2017/08/01 23:00    23.6
2017/08/02 00:00    23.0
2017/08/02 01:00    22.7
2017/08/02 02:00    22.5
2017/08/02 03:00    21.9
2017/08/02 04:00    21.7
2017/08/02 05:00    21.6
2017/08/02 06:00    22.3
2017/08/02 07:00    22.7
2017/08/02 08:00    23.3
2017/08/02 09:00    23.9
2017/08/02 10:00    24.7
2017/08/02 11:00    24.5
2017/08/02 12:00    24.3
2017/08/02 13:00    24.7
2017/08/02 14:00    24.3
2017/08/02 15:00    24.8
2017/08/02 16:00    24.8
2017/08/02 17:00    24.4
2017/08/02 18:00    23.9
2017/08/02 19:00    23.4
2017/08/02 20:00    23.1
2017/08/02 21:00    22.7
2017/08/02 22:00    22.3
2017/08/02 23:00    21.7
2017/08/03 00:00    21.3

課題(2)

  • 以下のようなグラフを作成せよ（テスト用として枠線のみとし、データのプロットは省略）。
  • x軸は、annotationが月名（頭文字のみ）で、ticsが10日間隔としている。

解答

課題(1)

  • 日付データの変換をExcelで行う場合、セルの表示形式を「yyyy-mm-ddThh:mm」とすればGMT形式となる。
  • 以下の例では、awkで変換を行う場合のスクリプトを記述しておきます。
  • awkの-Fオプションで区切り文字をtabに変更し、gsub関数で文字列を置換している。
plot_assign04-01.sh
#!/bin/bash
#--------------------------------------------------
# plot for GMT5
# 1次元時系列グラフ
#--------------------------------------------------

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

# GMTパラメーター
gmt gmtset FONT_ANNOT_PRIMARY=9p,Helvetica   # 主目盛りのフォント
gmt gmtset FONT_LABEL=10p,Helvetica          # 軸タイトルのフォント
# 日付関係
gmt gmtset FORMAT_DATE_MAP=mm/dd # 日付目盛り(yy[yy]=年,mm=月,dd=日,o=月名,jjj=日数)
gmt gmtset FORMAT_TIME_PRIMARY_MAP=a # 主目盛りの月・週表示：[f|a|c]

# input
in01=temp.txt                   # 入力ファイル
#in01_h=1                        # ヘッダー行数
#in01_c=0,1                      # 入力する列(0始まり）

# output
ps=plot_assign04-01.ps               # 出力psファイル

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

# グラフの範囲(range)（-Rオプション）
rw=2017-08-01T00:00             # x軸の最小値(west)
re=2017-08-03T00:00             # x軸の最大値(east)
rs=0                            # y軸の最小値(south)
rn=40                           # y軸の最大値(north)

# 軸・目盛りの設定
ti_x="Date (2017)"          # x軸タイトル
ti_y="Temperature [@~\260@~C]"  # y軸タイトル
tics_x=a1Df6h               # x軸の目盛り(a=主,f=副)
grid_x=g0                   # x軸のグリッド(g)(*)
tics_y=a10f2                # y軸の目盛り(a=主,f=副)
grid_y=g10                  # y軸のグリッド(g)(*)
axis=WeSn                   # 軸表示（大文字：軸＋目盛り、小文字：軸）

# 線の設定
line01=0.8p,black               # 線：太さ,色,線種

#----- plot ------------------------------
echo "plot $ps"

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

# psxy(line)
awk -F"\t" 'NR>1 {gsub("/", "-", $1); gsub(" ","T",$1); print $1, $2}' $in01 | \
    gmt psxy -J -R -W$line01 -K -O >> $ps

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

#----- 後処理 ------------------------------
# 画像変換(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)

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

# 中間ファイルの削除
rm -f gmt.history gmt.conf

課題(2)

  • psabasemapの-Bオプションに加え、FORMAT_DATE_MAP, FORMAT_TIME_PRIMARY_MAPの設定値を変更する。
plot_assign04-02.sh
#!/bin/bash
#--------------------------------------------------
# plot for GMT5
# 1次元時系列グラフ
#--------------------------------------------------

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

# GMTパラメーター
gmt gmtset FONT_ANNOT_PRIMARY=9p,Helvetica   # 主目盛りのフォント
gmt gmtset FONT_ANNOT_SECONDARY=9p,Helvetica # 副目盛りのフォント
gmt gmtset FONT_LABEL=10p,Helvetica          # 軸タイトルのフォント
# 日付関係
gmt gmtset FORMAT_DATE_MAP=o # 日付目盛り(yy[yy]=年,mm=月,dd=日,o=月名,jjj=日数)
gmt gmtset FORMAT_TIME_PRIMARY_MAP=c # 主目盛りの月・週表示：[f|a|c]

# input
#in01=temp.txt                   # 入力ファイル
#in01_h=1                        # ヘッダー行数
#in01_c=0,1                      # 入力する列(0始まり）

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

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

# グラフの範囲(range)（-Rオプション）
rw=2017-01-01T00:00             # x軸の最小値(west)
re=2018-01-01T00:00             # x軸の最大値(east)
rs=0                            # y軸の最小値(south)
rn=40                           # y軸の最大値(north)

# 軸・目盛りの設定
ti_x="Date (2017)"          # x軸タイトル
ti_y="Temperature [@~\260@~C]"  # y軸タイトル
tics_x=a1Of10d               # x軸の目盛り(a=主,f=副)
grid_x=g0                   # x軸のグリッド(g)(*)
tics_y=a10f2                # y軸の目盛り(a=主,f=副)
grid_y=g10                  # y軸のグリッド(g)(*)
axis=WeSn                   # 軸表示（大文字：軸＋目盛り、小文字：軸）

# 線の設定
line01=0.8p,black               # 線：太さ,色,線種

#----- plot ------------------------------
echo "plot $ps"

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

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

#----- 後処理 ------------------------------
# 画像変換(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)

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

# 中間ファイルの削除
rm -f gmt.history gmt.conf
4
3
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
4
3