時系列データをGMTで描画する
みなさんこんにちは。おこめまんです。
今回は「時系列グラフ」をGMTで書いてみます。
はじめに
みなさんがこのブログにたどり着いたということは、「時系列ってできるの？」とか「時刻上手くGMTで書けないんだが…」と思っていることだと思います。
気象データを時系列グラフにすることで、時間変化を追うことが出来たり、新しい発見があったりするものです。ぜひ、書けるようにしましょう。
さて、みなさんが書きたい図はこんな感じでしょうか。
この図は、ある地点のアメダスデータの風向・風速の時系列変化を表しています。
風向(右軸)はシンボル、風速(左軸)は実線です。
肝心の「時刻」(横軸)に注目すると、3/27 00:00 ~ 3/30 00:00となっています。
おそらく、皆さんがお困りなのは、この時刻をどう書かせたのかではないでしょうか。
シェルスクリプト
※fig_timeseries.shというシェルスクリプトを作成して、その中身を以下の通りとしました。赤色の部分は各々変更しましょう。
#
#!/bin/bash
#
######################################################################
### 作図日時指定 ###
SYEAR=2000
SMONTH=03
SDAY=27
SHOUR=00
SMINUTE=00
SSECOND=00
EYEAR=2000
EMONTH=03
EDAY=30
EHOUR=00
EMINUTE=00
ESECOND=00
### AMeDASデータ ###
ifile='amedas_200003.txt'
ifn=./${ifile}
### 風速の下限/上限指定 ###
minws=0
maxws=20
### 風向上限/下限指定 ###
minwd=0
maxwd=360
######################################################################
### 日時の結合 ###
starttime=${SYEAR}-${SMONTH}-${SDAY}T${SHOUR}:${SMINUTE}:${SSECOND}
endtime=${SYEAR}-${EMONTH}-${EDAY}T${EHOUR}:${EMINUTE}:${ESECOND}
### 作成する図＜風＞ ###
ofile='wind_'$SYEAR$SMONTH$SDAY'.ps'
mkdir -p ./figures/
ofn=./figures/${ofile}
### GMTのデフォルト変更 ###
#gmtdefaults -D
gmtset ANNOT_FONT_SIZE_PRIMARY 28p
gmtset HEADER_FONT_SIZE 35p
gmtset LABEL_FONT_SIZE 30p
gmtset INPUT_DATE_FORMAT yyyy-mm-dd
gmtset INPUT_CLOCK_FORMAT hh:mm:ss
gmtset PLOT_DATE_FORMAT mm-dd
gmtset PLOT_CLOCK_FORMAT hh:mm
############ 作図 ############
#####1軸＜風速＞ #####
psbasemap -R$starttime/$endtime/$minws/$maxws -Jx2.5i/0.5 -Bs1D::/:: -Ba12h:"JST":/2.5:"Wind Speed [m/s]":WS -X5 -Y5 -K > $ofn
### 風速時系列描画 ###
awk '$9 > 0.0 {print $7" "$9}' $ifn | \
psxy -Jx -R -W8 -K -O >> $ofn
### 2軸＜風向＞ ###
psbasemap -R$starttime/$endtime/$minwd/$maxwd -Jx2.5i/0.02777 -Bs1D::/:: -Ba12h:"JST":/60:"Wind Degree [\260]":ES -K -O >> $ofn
### 風向時系列描画 ###
awk '$10 > 0.0 {print $7" "$10}' $ifn | \
psxy -Jx -R -Sc0.5 -G100 -N -O >> $ofn
ps2raster -P -A -Tg $ofn
シェルスクリプトの説明
※大事なところを説明します！
### 作図日時指定 ###
SYEAR=2000
SMONTH=03
SDAY=27
SHOUR=00
SMINUTE=00
SSECOND=00
EYEAR=2000
EMONTH=03
EDAY=30
EHOUR=00
EMINUTE=00
ESECOND=00
→いきなりですが、この部分で描画したい範囲での時刻を指定します。今回は2000/3/27 00:00:00～2000/3/30 00:00:00の時系列グラフを作成します。いろんな日付の作図を行う際にはいちいち指定するのは大変なので、以前書いたこちらを参照して適宜書き換えましょう。
*****************************************************************************************************
### AMeDASデータ ###
ifile='amedas_200003.txt'
ifn=./${ifile}
→ここで作図するinputfileを指定しました。中身ですが、2000/3/27 00:00:00～2000/3/30 00:00:00のデータが1行ずつ記述されているものとします。シェルスクリプトと同じディレクトリに置きましょう。
アメダスデータ10分間値の取得方法はこちらを参照してください。
今回は、左から「年、月、日、時、分、秒、日時を結合したやつ、気温、風速、風向、降水、日照」となっているテキストデータ(コンマ区切り)を用意しました。
先ほどの図はこのG列、I列、J列を使用して作図しました。
はい。でました。「日時を結合したやつ」！！！！
ここがとっても大事です。詳しくみていきましょう。
1行目に注目。
2000-03-27T00:00:00
これは2000/3/27 00:00:00を表しています。注意したいのは間をつないでいる「T」です。この形にそろえましょう。
たとえば、Excelを使用して整形する場合…
整形前のifileには「日時を結合したやつ」なんて入っていないので付け足します。
①新しく列を作る(今回はG列を追加しました。)
②G1に「=DATE(A1,B1,C1)+TIME(D1,E1,F1)」を書き込む
→ここで、日時を結合しています。
③G列についてホーム→数値→その他の表示形式→ユーザー定義から「yyyy-mm-ddThh:mm:ss」のように変更する
→これでTを含んだかたちで表示されるようになります。
*****************************************************************************************************
### 日時の結合 ###
starttime=${SYEAR}-${SMONTH}-${SDAY}T${SHOUR}:${SMINUTE}:${SSECOND}
endtime=${SYEAR}-${EMONTH}-${EDAY}T${EHOUR}:${EMINUTE}:${ESECOND}
→ここで、psbasemapの横軸範囲指定を行っています。「作図日時指定」を「T」を含んだ形に変換します。ここでも「T」を含んだ形にしましょう。
*****************************************************************************************************
### GMTのデフォルト変更 ###
#gmtdefaults -D
gmtset ANNOT_FONT_SIZE_PRIMARY 28p
gmtset HEADER_FONT_SIZE 35p
gmtset LABEL_FONT_SIZE 30p
gmtset INPUT_DATE_FORMAT yyyy-mm-dd
gmtset INPUT_CLOCK_FORMAT hh:mm:ss
gmtset PLOT_DATE_FORMAT mm-dd
gmtset PLOT_CLOCK_FORMAT hh:mm
→特に橙色で示したところは必ずシェルスクリプト中に書きましょう。
*****************************************************************************************************
#####1軸＜風速＞ #####
psbasemap -R$starttime/$endtime/$minws/$maxws -Jx7.0i/0.5 -Bs1D::/:: -Ba12h:"JST":/2.5:"Wind Speed [m/s]":WS -X5 -Y5 -K > $ofn
→ベースマップの指定です。大事そうなところだけ。
-R$starttime/$endtime/$minws/$maxwsは図の範囲を指定しています。$starttime/$endtimeで横軸(時間)、$minws/$maxwsで縦軸(風速)です。
-Bs1D::/::のs1Dは横軸の月日を表示させるためのものです。この設定だと1日ごとの表記になります。
-Ba12h:"JST":/2.5:"Wind Speed [m/s]":WSのa12hは横軸の目盛を振る間隔です。12時間ごとの場合はこの表記になります。
*****************************************************************************************************
その他のところは、図の体裁にかかわる部分です。各自好きに変えましょう。
おわりに
時系列グラフを作成する場合、「T」を付けて日時を結合することがポイントとなります。ifile、psbasemap両方。
ifleさえ準備できれば、シェルスクリプト完コピで作図されるはずです。お試しください。
