動機
花粉がいつどこで最盛期を迎えているか分かったら面白いかな，と思って初めてみました．
完成形
1フレーム1時間に対応し，2023年3月頭から4月中旬までの様子になります．
追記：たまに見られないです．githubの方で見てください．
（クロップしてあります．実際は横にカラーバーと上に日付時刻が表示され，余白が大きいです）
コードの概要
花粉のデータはWeathernewsが提供しているデータを使用します．
APIが提供されており，取得したい期間と市町村コードを使うことで簡単に取得することができます．
しかし，一回で取得できるデータが市町村1つ，期間は１ヶ月までなので，いくら自動化しても繰り返しの取得には時間がかかります．
そこで，一度取得したデータをローカルに保存するようにし，プログラムの構成を以下のような構成にします．
- 取得したい市町村（都道府県単位），期間を指定
- ローカルのデータベースにアクセスして，なければたりない分（市町村や期間）をオンラインで取得
- 取得したデータをローカルのデータベースに追記
- ローカルにあるデータベースをもとに花粉の分布を描画
データベースの更新と，グラフの描画それぞれにそれなりの時間がかかりそうだったので，これらを別々のスクリプトにわけました．
なお，今回の結果では，花粉飛散量のデータをWeathernewsが提供する一市町村あたり一地点のデータとしています．あくまで観測器が設置された一地点のデータで色分けをしているので，正確な分布を示しているわけではありません．サイエンティフィックな解析ではなく，テクニカルなエンタメとしてお楽しみいただければ幸いです．
使用したライブラリ
- requests: WeathernewsのAPIからデータを取得する用
- geopandas: 地図情報をPandasライクに操作できるようにしたライブラリ（らしい？）
- pandas: データハンドリング用
- h5py: オンラインで取得したデータをローカルにHDF5形式で保存するため
- numpy: 配列操作用
- opencv: 画像から動画を作成
これとは別に，市町村ごとの地形データが必要になります．国土交通省が配布する国土数値情報を利用します．
ソースコードは日本全土のデータを使っています．ファイルの読み取りの箇所で適宜ファイル名を変更すれば各都道府県などの描写も可能です．
工夫した点
データ収集 (collect_data.py)
まずは花粉データの取得を行います．花粉データの取得には，市町村コードとデータ期間（開始日，終了日）が必要なので，市町村コードをどうやって取得するか考えます．
一応，Weathernewsのページに一覧があるのですが，後々のプロットを考えると，プロットに使う地図データから引き出すのがいいかなと思いました．
まずは地図データから市町村コードを抜きだす関数です．引数として地形データファイル名と，対象とする都道府県のリストを受け取り，それをもとに市町村コードのリストを返します．空の都道府県リストが渡された場合は日本全国を対象にします．
def listup_city_code(gdf_jp, prefectures):
'''
list up city code from geometry data
Input:
- gdf_jp : geopandas data to read
- prefectures : list of prefectures to list up city code,
if empty, list up all city codes in japan
Output:
- list of city codes in specified prefectures
'''
### list of city code list to return
cc_list = []
if prefectures:
### if prerectures are specified
### loop for prefectures
for p in prefectures:
### extract geopandas dataframe with one prefecture
gdf = gdf_jp[ gdf_jp['N03_001'] == p ]
for cc in gdf['N03_007']:
### add city code if not in the list yet
if cc not in cc_list:
cc_list.append(cc)
else:
### the same loop above for all prefectures
for cc in gdf_jp['N03_007']:
if cc not in cc_list:
cc_list.append(cc)
return cc_list
あとは市町村コードのリストに則ってデータを集めます．まずはローカルのデータを参照し，ローカルにない分のデータをオンラインからとってきます．
実際にオンラインからデータを取得する関数はこれです．引数に文字列として取得期間の開始と終了，それと市町村コードを受け取ります．APIで取得したデータをPandasのデータフレームに変換して返り値とします．ifの分岐は，指定した市町村コードのデータがWeathernewsのデータベースにないとNot foundのエラーが帰ってくるので，花粉飛散量が-1のダミーのデータフレームを返すようにします．こうしておくと後のプロットが楽になる気がします．
def get_data(start:str, end:str, citycode='04103'):
### get pollen date from wethernews api in one month
url = 'https://wxtech.weathernews.com/opendata/v1/pollen?'
payload = {"citycode": citycode, "start": start, "end": end}
r = requests.get(url, params=payload)
# print(citycode)
if r.ok:
### if request is success, prepare dataframe to return
### convert to dataframe
df_temp = pd.read_csv(
io.StringIO(r.text),
sep=',',
dtype={'citycode': str, 'date': str, 'pollen': int}
)
df = pd.DataFrame({
'Date' : df_temp['date'],
citycode: df_temp['pollen']
})
### convert datetime text to datetime type
df['Date'] = pd.to_datetime(df['Date'], format='%Y-%m-%dT%H:%M:%S%z')
return df
else:
### if request is failed, return '-1' data for all the dates as dummy dataframe
end_point = datetime.datetime.strptime(end + ' 23:00:00', '%Y%m%d %H:%M:%S')
ser_date = pd.Series(
pd.date_range(
start,
end=end_point,
freq='H',
tz='Asia/Tokyo'),
name='Date'
)
ser_pollen = pd.Series(data=np.ones(len(ser_date), dtype=np.int64)*-1, name=citycode)
df = pd.concat([ser_date, ser_pollen], axis=1)
return df
実際は，取得期間の長さに上限（最大31日間）があるので，期間を指定したらその期間を１ヶ月単位に分割して繰り返しデータ取得するような関数を用意しておきます．
def get_city_pollen_data(start_date, end_date, city_code):
'''
collect pollen date during specified term between start_date and end_date in one city
'''
month_list = create_month_list(start_date, end_date)
df = pd.DataFrame()
### get monthly data
for ds in month_list:
print('{0:%Y%m%d} - {1:%Y%m%d}'.format(ds[0], ds[1]))
df_c = get_data(
'{0:%Y%m%d}'.format(ds[0]),
'{0:%Y%m%d}'.format(ds[1]),
citycode=city_code
)
df = pd.concat([df, df_c], axis=0)
return df
なお，関数create_month_listは指定期間を１ヶ月ごとに区切り，１ヶ月の[開始日，終了日]のリストを月数分リストとして返す関数です．
ダウンロードしたデータはローカルにHDF5として保存します．HDF5は個人的に最近気に入って使っているファイル保存形式で，バイナリ形式で保存できて読み書きが早く，階層構造を採用できるのでデータ整理が比較的しやすくて気に入っています．h5pyはHDF5をpythonで扱うためのライブラリです．
ファイル構造は以下のようにします．
year/
├citycode0
│ ├Date
│ └Pollen
├citycode1
│ ├Date
│ └Pollen
:
:
年ごとのメイングループ分けを行い，市町村コードごとのサブグループを用意します．市町村グループ内には取得したデータの日付と花粉情報をそれぞれ格納しておきます．
個人的に工夫したのは，日付情報をUNIX時間として浮動小数点データとして保存したことです．pythonのDatetime型はHDF5にそのまま書き込むことはできないので，別の型に変換する必要があります．書き込みと読み出し時に変換が必要ですが，まあ許容範囲内ですね．
データ収集のメイン関数は少し長いので割愛します．基本は，
- ローカルデータの読み込み
- 指定市町村のデータがローカルにあるかどうか確認
- なければオンラインから取得してローカルデータに追記
の繰り返しです．
分布図描画
続いて，ローカルに保存したデータをもとに花粉分布を描画していきます．
こちらは最初にメインの関数を紹介します．
while文の前では
- 出力画像を保存するディレクトリの用意
- プロットする期間をローカルデータに基づいて決定（この時に花粉データの最大値も取得）
- matplotlibのfigureとカラーマップを用意
- カラーマップはjet
- normalizeは最小値0で規格化（負の値は意味のないデータのため）
をしています．
while文の最初では，対象の日付のプロットがすでにあるかどうか判定しています．プロットにも時間がかかるので，できるだけ描画回数を増やさないようにしています．
if文の中は花粉データを読み込んで，それに応じて色分けをして描画するだけです．
動作的には，毎ループローカルのデータベースにアクセスしているので，あまり賢い方法ではないのですが，多分描画の方が時間がかかるのでいいかな，と思っています．メモリにも優しいかな，と．
def plot_pollen_map(year, gdf, city_codes, hdf5_database):
'''
plot maps based on hdf5 database
'''
### directory to save drawn map
png_dir = 'pollen_map_japan_png'
if not os.path.isdir(png_dir): os.makedirs(png_dir)
### get the datetime of beginning and end and maximum pollen data
dt0, dt1, max_pol = get_plot_duration(year, city_codes, hdf5_database)
print('Plot range {0:%Y-%m-%d} - {1:%Y-%m-%d}'.format(dt0, dt1))
### variable of current datetime of plot
plot_dt = dt0
### prepare campus for plot
fig = plt.figure(figsize=(12, 12))
### use jet for coloring
cmap = mplcm.jet
### norm with maximum pollen data
norm = mplcolors.Normalize(vmin=0, vmax=max_pol)
while plot_dt <= dt1:
### define the png file name
png_filename = png_dir + '/PollenMap_{0:%Y-%m-%d_%H%M}.png'.format(plot_dt)
if not os.path.isfile(png_filename):
### if there is no file, drawing start
print('\rPlot on {0:%Y-%m-%d %H:%M}'.format(plot_dt), end='')
### obtain the pollen data of the current date and hour
ser_pol = get_pollen_data(year, plot_dt, city_codes, hdf5_database)
### prepare list for coloring based on pollen number
fc = [ ser_pol[cc] for cc in gdf['N03_007'] ]
### prepare axis
ax = fig.add_subplot(111)
ax.axis('off')
ax.set_aspect('equal')
### draw map distribution
gdf.plot(
ax = ax,
edgecolor = 'black',
facecolor = mplcm.ScalarMappable(norm=norm, cmap=cmap).to_rgba(fc),
linewidth = 0.1
)
### put colorbar
fig.colorbar(
mplcm.ScalarMappable(norm=norm, cmap=cmap),
orientation = 'vertical',
fraction = 0.05,
aspect = 50
)
### format the final apperance of map
fig.suptitle('{0:%Y-%m-%d %H:%M}'.format(plot_dt))
plt.tight_layout()
# plt.show()
plt.savefig(png_filename)
plt.clf()
### iteration for next datetime (one hour step)
plot_dt += datetime.timedelta(hours=1)
print('')
これで1時間ごとの花粉分布マップが画像として保存されるので，それを動画にするだけです．
一緒にあるmakemovie.pyを実行することで，連続した画像を動画として出力します．\
import cv2
import glob
import tqdm
def make_movie(image_file_path, frame_rate=30.0, video_name='video', image_format='png'):
filelist = glob.glob(image_file_path + '/*.' + image_format)
filelist.sort()
im0 = cv2.imread(filelist[0])
height, width, color = im0.shape
fourcc = cv2.VideoWriter_fourcc('m', 'p', '4', 'v')
video = cv2.VideoWriter('{0}.mp4'.format(video_name), fourcc, frame_rate, (width, height))
for file in tqdm.tqdm(filelist):
img = cv2.imread(file)
video.write(img)
if __name__ == '__main__':
make_movie('pollen_map_japan_png', frame_rate=10.0, video_name='pollen_movie_japan')
最後に
個人的に毎年花粉症に悩まされているので，描画結果を見ているだけでも少しむずむずしました．
南から徐々に北上していくのかな，と思っていたんですが，飛散量の地域差が意外と小さいのかなと思いました．
参考文献
- 地形データのダウンロード先．重いですが全国分入れとくとどの都道府県も描画できるのでおすすめ
https://nlftp.mlit.go.jp/ksj/jpgis/datalist/KsjTmplt-N03.html - Weathernewsの花粉データ．2022年分もあります．なお，オンライン公開の許可は取れています．
https://wxtech.weathernews.com/pollen/index.html
-UNIX時間とpython datetimeの相互変換．その他の記事も毎回非常に参考にしております
https://note.nkmk.me/python-unix-time-datetime/ - Datetimeのタイムゾーンの指定．上記の変換を確実に意図したものにするために．
https://qiita.com/keisuke0508/items/df2594770d63bf124ccd