間が空きましたが生きています。
久しぶりの投稿は意外と忘れやすい時系列グラフの描き方の復習です。
慣れていない方はbasemapの記事を最初に見たほうがよいです。
目標
2001年から2021年までの月別地震回数をプロットしたい。
図示までの流れ
- 震源データを取得
- カンマ区切りにでもする
- PandasとかAWKとかで月ごとの回数を抽出
- 図示
震源データの取得
今回は気象庁一元化震源データを使用します。
FAR FIELDの除外は面倒なので省きます。
年別地震回数をさっくり描く
年別でよいなら簡単に作成可能です。
#!/bin/bash
rm year_count.csv
for year in {2001..2022};
do
times=`cat h${year} | wc -l `
echo "${year},${times}" >> year_count.csv
doneでもいいし、
import pandas as pd
import numpy as np
cols = ["year", "count"]
df = pd.DataFrame(index=[], columns=cols)
for i in range(1998, 2022):
filename = "h" + str(i)
count = sum([1 for _ in open(filename)])
print(count)
record = pd.Series([str(i), count], index=df.columns)
df = df.append(record, ignore_index=True)
df.to_csv("year_count.csv", index=None, header=None)でもよいです。
これで年別回数のデータが作成できました。
あとは図示するだけ。これは時系列でなくても行けそうですね。
#!/bin/bash
gmt begin count_year png
gmt basemap -JX16/8 -R2001/2021/0/400000 -Bxa2f1g+l"year" -Byafg+l"count" -BWSne
gmt plot "./year_count.csv" -Sc0.1 -Gred -W0.25 -N
gmt plot "./year_count.csv" -W0.25,red -N
gmt endこの結果が、
となります。
- 2011年は東北地方太平洋沖地震か。
- 2016年は…熊本地震の影響？
- 2021年も地震の回数が増えていたんですね。
月別で表せるようにデータを整理する
さて、本題の月別です。
まずは、気象庁の震源データをカンマ区切りにしなければなりません。
いつまで固定長のみの提供なんですかね、、
自分でも疑わしいですが、上記のスクリプトを正として、
カンマ区切りにしたファイルから月別で回数を取得します。
ついでに、FAR FIELDやUSGSが決定した震源は震源が日本でないことが多いのでサクラエディタなどで、
^U
FAR FIELD
で行削除をしておきましょう。
あとは月別にすればよいです。私だったらPythonで
import numpy as np
import pandas as pd
cols = ["month", "count"]
df = pd.DataFrame(index=[], columns=cols)
cols = ["recode","year","month","day","hour","minute","second","time_SE","latitude_deg","lat_SD","longitude_deg","lon_SD","depth","dep_SD","mag1_1","mag1_type","mag2","mag2_type","travel_time","hypo_judge","hypo_type","max_coefficient","scale_victim","tsunami_victim","region_main","region_sub","hypo_locale","obs_number","determination_way"]
df_hypo = pd.read_csv("hoge.txt", names=cols)
for i in range(2001, 2022):
count_yeardf = df_hypo[df_hypo['year'] == i]
for j in range(1,13):
count_monthdf = count_yeardf[count_yeardf['month'] == j]
output_month = str(i) + "-" + str(j)
record = pd.Series([output_month, len(count_monthdf)], index=df.columns)
df = df.append(record, ignore_index=True)
df.to_csv("month_count.csv", index=None, header=None)とすれば、あとは図示するだけです
図示する
gmt basemap -JX[ヨコの長さ]/[タテの長さ] -R[描画範囲] -B[x|y|z][a][y|o|d|etc.][f][y|o|d|etc.][g][y|o|d|etc.][+l”label”] -B[WSEN][+t”title”]
※ タテの長さを指定しない場合はヨコの長さと同じ
※ 日時は少なくとも日にちまで必要 (例 : 2020-04-01)。
※ 目盛り間隔は数値のあとに”y”を付けると年、”o”を付けると月、”d”を付けると日。
※ "y"は西暦の下2ケタ表記。"Y"と大文字にすると4ケタ表示
今回は月別なので-Rの範囲は12月1日まででよいです。
X軸も主目盛線を年毎、副目盛線を6カ月とすれば分かれば何となくわかるでしょう。
#!/bin/bash
gmt begin count_month png
gmt basemap -JX32/16 -R2001-01-01/2021-12-01/0/60000 -Bxa1Yfg6o+l"year-month" -Byafg+l"count" -BWSne
gmt plot "./month_count.csv" -Sc0.1 -Gred -W0.25
gmt plot "./month_count.csv" -W0.25,red
gmt endとすれば、
となります。
ついでに、
- -Bxa1yfg6oの場合：2ケタ表記
- -Bxa6Ofgの場合
これではメチャクチャなので、、
- -Bxa6Ofg+a45の場合（目盛り文字を45度の角度にする）
これで見やすいですね。
コメント