時系列データの作図

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