こんにちは、とりりんです。
今回は、GMTを使って2021年の台風データの可視化をしてみたいと思います。
こんな方向けに書いています！
- GMT初心者の方
- データ解析をしている方
- データの可視化をしたい方
是非最後まで読んでいただけると嬉しいです！
台風データ
台風の有名なデータとして、ベストトラックデータというデータがあります。
ベストトラックデータには、各台風の位置や中心気圧、最大風速、暴風域・強風域の半径などがあり、テキスト形式で記録されています。
データはこちらから取得できます。
今回は中心気圧(データの６列目)と最大風速（同７列目）を使って時系列（時刻はデータ１列目）変化を見ていこうと思います。
基本的な話ですが、強い台風に発達するほど、中心気圧が低く、最大風速が強くなる傾向になります。
今回は、発達や衰退指定く様子がグラフで読み取れたら合格というところですね。
実際の処理
では、GMTのpsxyコマンドをベースに作成していきましょう。
psxyコマンドについて知りたい方は、こちらをご覧ください。
また、今回紹介するソースコードはGitHubに置いてあるので、実際に動かしたい方は使ってみてください。
スポンサーリンク
台風ごとにデータを分割
まず、ダウンロードしたベストトラックデータは2021年すべての台風が載ってありました。
今回使用する上で台風ごとに分かれていて欲しかったので、まずは台風ごとに分解していきます。
ソースコードはGitHub上のseparate.shです。
#!/bin/bash
data="./data/typhoon2021.txt"
stdDir="./data"
while read line; do
time=`echo $line | awk '{printf $1}'`
# ヘッダ行の場合
if [ $time -eq 66666 ]; then
# Ex. 2021年１号 -> 2101
tyNum=`echo $line | awk '{printf $2}'`
# 保存先ディレクトリがなければ新規作成
outDir=$stdDir/$tyNum
if [ ! -d $outDir ]; then
mkdir $outDir
fi
# 保存先パスに既にデータがあれば削除
outPath=$outDir/$tyNum".txt"
if [ -e $outPath ]; then
rm -f $outPath
fi
echo $outPath
fi
echo $line >> $outPath
done < $data
このソースコードは、ダウンロードしたデータ（今回は./data/typhoon2021.txt）を使って各台風ごとのデータ（１号の場合は./data/2021/2021.txt）を作成しています。
まず、while文でデータを１行ずつ読み込み、その１行は変数$lineに格納されます。
echo $line | awk ‘{printf}’がよく使われていますが、これはline指定した列を取得できます。
また、１列目に時刻があるのですが、ヘッダ行は66666になっています（11行目）。
これを利用して、ヘッダ行の場合にファイルの保存パス（./data/2021/2021.txtなど）を作って、１行ずつ保存していきます（31行目）。
ヘッダ行のif文を無視すれば、ファイルを読み込んで書き込む作業をしているだけです。
スポンサーリンク
グラフを作成
下処理が終わったので後は書くだけ！と思いたいですが、下処理をまだしています。
後、めんどくさいので色々自動化してしまいました。
#!/bin/bash
tyYear=21
for n in `seq -f %02g 8`; do
tyNum=$tyYear$n
dataDir=./data/$tyNum
tyData=$dataDir/${tyNum}.txt
presTxt=$dataDir/${tyNum}_Pressure.txt
mswTxt=$dataDir/${tyNum}_MaximamSustainedWind.txt
imgDir=./img/$tyNum
if [ ! -d $imgDir ]; then
mkdir $imgDir
fi
savePs=$imgDir/${tyNum}.ps
savePng=$imgDir/${tyNum}.png
if [ -e $mswTxt ]; then
rm -f $mswTxt
fi
if [ -e $presTxt ]; then
rm -f $presTxt
fi
if [ -e $savePs ]; then
rm -f $savePs
fi
if [ -e $savePng ]; then
rm -f $savePng
fi
count=0
while read line; do
# ヘッダ行読み飛ばし
if [ $count -ne 0 ]; then
# GMT用の日付表示に変更（Ex. 21010100 -> 2021-01-01T00:00:00）
date=`echo $line | awk '{printf $1}'`
year=20${date:0:2}
month=${date:2:2}
day=${date:4:2}
hour=${date:6:2}
gmtDate=${year}-${month}-${day}T${hour}:00:00
# 描画用にデータ開始時刻を保存
if [ $count -eq 1 ]; then
startDate=${gmtDate:0:11}00:00:00
fi
hpa=`echo $line | awk '{printf $6}'`
msw=`echo $line | awk '{printf $7}'`
echo $gmtDate $hpa >> $presTxt
echo $gmtDate $msw >> $mswTxt
fi
count=$((count+1))
done < $tyData
# データ表示間隔調整(描画用)
diffDay=$((${gmtDate:8:2}-${startDate:8:2}))
if [ $diffDay -lt 0 ]; then
diffDay=$((30-${startDate:8:2}+${gmtDate:8:2}))
fi
timeInterval=`echo "scale=5; 0.00031 / ${diffDay}" | bc`
# 表示間隔変更（描画用）
if [ $diffDay -ge 10 ]; then
dispTime=a3Df12Hg1D
elif [ $diffDay -lt 2 ]; then
dispTime=a1Df1Hg3H
else
dispTime=a1Df3Hg6H
fi
gmt psxy $presTxt -B$dispTime:"UTC":/a20g10:"Pressure[hPa]"::."20${tyNum:0:2} No.${tyNum:2:2}":WSn -R$startDate/$gmtDate/940/1020 -Jx${timeInterval}T/0.15 -Sc0.2 -Gblue -K > $savePs
gmt psxy $presTxt -B -R -J -W1,blue -P -K -O >> $savePs
gmt psxy $mswTxt -B/a10:"MSW[kt]":E -R$startDate/$gmtDate/0.1/100 -Jx${timeInterval}T/0.12 -Ss0.2 -Ggreen -P -K -O >> $savePs
gmt psxy $mswTxt -B -R -J -W1,green -P -K -O >> $savePs
gmt pslegend -R -J -Dx1.8/11.5/2.5/1.2 -F+p1,black+gwhite -P -O <> $savePs
S 0.2 c 0.2 blue black 0.5 Pressure
S 0.2 c 0.2 green black 0.5 MSW
EOF
gmt ps2raster $savePs -E100 -Tg
echo $tyNum Save
done
描画は終盤のとしているところだけです。
それまでは下処理です。
実際の処理を説明すると、台風の数だけ繰り返し（for文）とファイルの読み込み（while文）をしています。
６行目から32行目はファイル指定や準備をしているだけなので、気にしないでください。
34行目からのwhile文では、時刻, 中心気圧と時刻, 最大風速のデータを作成しています。
この形式にしないとGMTさんが働いてくれないので…。
また、時刻を2021-01-01T00:00:00などの形式にしないといけないので、その形に直しています（40~46行目）。
その後、中心気圧と最大風速を取得（53, 54行目）し、ファイルに追記しています（55, 56行目）。
>>は追記し続けるため、while文に入る前にrm -fでファイルを消したりしているわけです。
ようやくここからが本番ですが、今まで作成したファイルを使ってGMTで描画していきます。
例として、書き込んだ気圧データはこんな感じになっています。
2021-02-16T06:00:00 1004
2021-02-16T12:00:00 1006
2021-02-16T18:00:00 1004
…
x, yデータとなっていますね（今回は時系列なのでxではなくtが適切ですが）。
前回学習したpsxyコマンドを使用して、中心気圧と最大風速について、線と丸を描画します（79~82行目）。
ついでに凡例をつけて完成です。
まとめ
今回はGMTのpsxyコマンドを使ってグラフを作成してみました。
下処理がすごく長かったですが、GMT自体は10行足らずでできました。
GMTは良い描画ツールなので、是非マスターしてください。
最後までご覧いただき、ありがとうございました！
Python
↓とりりんを応援して頂けると嬉しいです！↓
にほんブログ村
コメント