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

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

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

<img src="https://qiita-image-store.s3.amazonaws.com/0/195213/d43f6490-84cf-4e1a-84b1-d09c3b2a8675.png" width=500>


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

```shell: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](http://www.data.jma.go.jp/obd/stats/etrn/index.php)のデータを加工して作成した（実際の観測値を保証するものではありません）。

```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形式](https://ja.wikipedia.org/wiki/ISO_8601)の拡張形式と呼ばれる。
+ 使用しない時間部分は省略できる。例えば、上記のデータでは秒の部分（: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軸は次のようになる。

<img src="https://qiita-image-store.s3.amazonaws.com/0/195213/b1cfa976-0ddc-9af6-1f9c-8ff9d32895b5.png" width=500>

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

<img src="https://qiita-image-store.s3.amazonaws.com/0/195213/042840b7-ece4-9a71-5d16-5e296eb5bc78.png" width=500>

+ 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、とした場合は、次のようになる。

<img src="https://qiita-image-store.s3.amazonaws.com/0/195213/f4f9a98b-067d-db16-5814-24b1872b0d26.png" width=500>


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

<img src="https://qiita-image-store.s3.amazonaws.com/0/195213/b2a69240-f70d-5852-5f0d-f2b93045a753.png" width=500>

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



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

<img src="https://qiita-image-store.s3.amazonaws.com/0/195213/9960938b-9a28-7f89-c06b-5c470ed641e5.png" width=500>

```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日間隔としている。

<img src="https://qiita-image-store.s3.amazonaws.com/0/195213/08b25ad9-e334-607c-9e20-ea7bb165fab7.png" width=500>

## 解答

### 課題(1)
+ 日付データの変換をExcelで行う場合、セルの表示形式を「yyyy-mm-ddThh:mm」とすればGMT形式となる。
+ 以下の例では、awkで変換を行う場合のスクリプトを記述しておきます。
+ awkの-Fオプションで区切り文字をtabに変更し、gsub関数で文字列を置換している。

``` shell: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の設定値を変更する。

``` shell: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
```