気象過去データの利用３（豪雨時の降水量を時系列ヒートマップ表示）

気象庁が、2020年3月末まで、気象過去データを無償で提供しています。
(参考) 気象過去データの利用環境
https://www.data.jma.go.jp/developer/past_data/index.html

基本的な気象データは、「利用目的・対象を問わず、どなたでもご利用頂けます。」とのことなので、気象データを使って何か行っていきます。

気象庁の「過去の気象データ・ダウンロード」ページからもダウンロードできますが、まとめて一括してダウンロードできるため非常に便利です。
利用可能なデータは、以下の掲載されています。
https://www.data.jma.go.jp/developer/past_data/data_list_20200114.pdf
もうすぐ期限ですので、必要な方は早めのダウンロードを

今回は、災害時の雨量の推移を地図上にヒートマップで表示してみます。
昨年の台風19号を想定していましたが、2019年7月までしか公開されていないため、2018年の7月豪雨のデータにします。
なお、地図描画用にfolium(0.10.1)を使用しました。

データダウンロード

アメダス各地点の毎時の降水量を取得するため「地域気象観測（アメダス）」－「時・日別値」を利用します。
このファイルには、各アメダス地点の降水量や気温などの観測情報が格納されています。
ファイルの形式については、以下を確認してください。
http://data.wxbc.jp/basic_data/kansoku/amedas/format_amedas.pdf

amedasフォルダにファイルをダウンロードします。約2GBの容量があるため時間がかかります。

import os
import urllib.request

# 地上気象観測　時・日別値ファイルダウンロード
url    = 'http://data.wxbc.jp/basic_data/kansoku/amedas/hourly_daily_1976-2019_v191121.tar'
folder = 'amedas'
path   = 'amedas/hourly_daily_1976-2019_v191121.tar'
# フォルダ作成
os.makedirs(folder, exist_ok=True)
if not os.path.exists(path):
    # ダウンロード
    urllib.request.urlretrieve(url, path)

ファイルは、tar形式のため中身のファイルを確認します。

# ファイル確認
import tarfile

# tarファイルに含まれるファイル確認
with tarfile.open(path, 'r') as tf:
    for tarinfo in tf:
        print(tarinfo.name)

中身は、hourly_daily_1976.tar.gzのようなtar.gzファイルでした。年ごとにtar.gzに固められています。
さらに展開が必要です。tar.gzファイルの内容を確認してみます。

# ファイル確認2
import tarfile

# tarファイルに含まれるファイル確認
with tarfile.open(path, 'r') as tf:
    for tarinfo in tf:
        print(tarinfo.name)
        if tarinfo.isfile():
            # tar.gzファイルに含まれるファイル確認
            with tarfile.open(fileobj=tf.extractfile(tarinfo), mode='r') as tf2:
                for tarinfo2 in tf2:
                    if tarinfo2.isfile():
                        print('\t' + tarinfo2.name)

以下のようなファイルが存在することがわかります。
hourly_daily_1976/area57/ahd1976.57246
ファイル名は、以下の形式とのことです。
ahdYYYY.SSSSS （YYYY：西暦年、SSS：観測所番号）

降水量データの読み込み

2018年7月5日～8日までの全国のアメダス地点の降水量を取得します。
tarファイル、tar.gzを展開せずに読み込みを行います。
各ファイル1日分が1300バイトになります。7月4日まで読み飛ばします。
降雪のみ測定している地点は除外します。また、正常値(RMKが8)以外のデータは、欠損として処理します。

データを格納用のpandasのデータフレームを用意します。
年月日時間を行、観測所番号を列としてデータを格納します。

# データ格納用データフレーム作成
import pandas as pd

prec_df = pd.DataFrame()

# 降水量データ取得
import tarfile
import struct

# 2018年7月5日～8日
year = '2018'
month = 7
start_date = 5
end_date   = 8
head_size = 28
time_size = 48
day_size  = 120
r_size = 1300

# tarファイルに含まれるファイルを取得
with tarfile.open(path, 'r') as tf:
    for tarinfo in tf:
        if tarinfo.isfile():
            # 2018年のファイルを抽出
            if tarinfo.name[13:17] == year:
                # tar.gzファイルに含まれるファイルを取得
                with tarfile.open(fileobj=tf.extractfile(tarinfo), mode='r') as tf2:
                    for tarinfo2 in tf2:
                        if tarinfo2.isfile():
                            print(tarinfo2.name)
                            # ファイルをopen
                            with tf2.extractfile(tarinfo2) as tf3:
                                # 日付まで読み飛ばし
                                buff = tf3.read((month-1)*31*r_size)
                                buff = tf3.read((start_date-1)*r_size)
                                # 日付分読み込み
                                for i in range(start_date, end_date+1):
                                    buff = tf3.read(head_size)
                                    sta_no, sta_type, year, monte, date = struct.unpack('<ih16xhhh', buff)
                                    print(sta_no, sta_type, year, month, date)
                                    # 降雪のみの観測所は対象外
                                    if sta_type != 0:
                                        for time in range(24):
                                            buff = tf3.read(time_size)
                                            prec, rmk = struct.unpack('<hh44x', buff)
                                            # 欠測確認
                                            if rmk != 8:
                                                prec_df.loc['{:04}{:02}{:02}{:02}'.format(year, monte, date, time), str(sta_no)] = None
                                            else:
                                                prec_df.loc['{:04}{:02}{:02}{:02}'.format(year, monte, date, time), str(sta_no)] = prec * 0.1
                                        buff = tf3.read(day_size)
                                    else:
                                        # 残りのデータを読み飛ばす
                                        buff = tf3.read(r_size-head_size)

少し時間がかかりますが、降水量データが取得できました。
データを確認してみましょう。

prec_df

参考までの欠損の数を確認します。

prec_df.isnull().values.sum()

欠損は、492個ありました。

観測所の緯度、経度取得

アメダス地点情報ファイルダウンロード

アメダスの地点情報ファイルをダウンロードします。

import os
import urllib.request

# アメダス地点情報ファイルダウンロード
url    = 'http://data.wxbc.jp/basic_data/kansoku/amedas/amdmaster_201909.tar.gz'
folder = 'amedas'
path   = 'amedas/amdmaster_201909.tar.gz'
# フォルダ作成
os.makedirs(folder, exist_ok=True)
if not os.path.exists(path):
    # ダウンロード
    urllib.request.urlretrieve(url, path)

ファイルを確認します。

# ファイル確認
import tarfile

# tarファイルに含まれるファイル確認
with tarfile.open(path, 'r') as tf:
    for tarinfo in tf:
        print(tarinfo.name)

観測所の緯度、経度の読み込み

amdmaster.indexに観測地点の履歴が格納されています。
各地点の履歴のため、観測日の緯度、経度を取得してもよいのですが、表示上、観測地点の移動は誤差範囲のため最新の緯度、経度を観測点の緯度、経度とします。
pandasのデータフレームに格納します。

# アメダス地点情報格納用データフレーム作成
import pandas as pd

amedas_df = pd.DataFrame()

# 観測所の緯度、経度取得
with tarfile.open(path, 'r') as tf:
    for tarinfo in tf:
        if tarinfo.name == 'amdmaster_201909/amdmaster.index':
            # ファイルをopen
            with tf.extractfile(tarinfo) as tf2:
                # ヘッダ読み飛ばし
                line = tf2.readline()
                line = tf2.readline()
                line = tf2.readline()
                while line:
                    # 最新の緯度経度を取得
                    prec_no = int(line[0:5])
                    amedas_df.loc[prec_no, '観測所名'] = line[6:26].decode('shift_jis').replace("　", "")
                    amedas_df.loc[prec_no, '緯度'] = int(line[74:76]) + float(line[77:81])/60
                    amedas_df.loc[prec_no, '経度'] = int(line[82:85]) + float(line[86:90])/60
                    line = tf2.readline()
            break

観測所の緯度、経度表示

amedas_df

ヒートマップ表示用データ作成

各地点ごとに[緯度, 経度, weight]のリストを作成し追加します。
weightには、降水量を指定します。ただし、0<weight<=1の値を指定する必要があるため、降水量を100で割ることにします。降水量が100mm以上の場合は、100mmとして計算しました。
また、降水量が0または欠損の場合は、地点に含めません。
[緯度, 経度, 降水量]を時間ごとにリストにします。
以下のような３次元のリストになります。
[[[緯度, 経度, 降水量]]]
indexには、時系列表示用の年月日時間を格納します。

import numpy as np

# 時間ごとに、各地点の緯度、経度、降水量データ生成
data = []
index = []
for time in prec_df.index:
    point = []
    index.append(time)
    for sta_no in prec_df.columns:
        #print(sta_no, time)
        # 欠損の場合は対象外
        if np.isnan(prec_df.loc[time, sta_no]): continue
        # 降水量なしは対象外
        if prec_df.loc[time, sta_no] == 0: continue
        # 緯度、経度を取得
        latitude =  None
        longitude = None
        day = time[0:8]
        latitude = amedas_df.loc[int(sta_no), '緯度']
        longitude = amedas_df.loc[int(sta_no), '経度']
        # ヒートマップ用に0～1の範囲にするため、降水量を1/100する。100mm以上は、100mmとする。
        prec = min(prec_df.loc[time, sta_no], 100) / 100
        point.append([latitude, longitude, float(prec)])
    data.append(point)

データを確認します。

data

時系列ヒートマップ表示

foliumを利用し、時系列ヒートマップを表示します。HeatMapWithTime関数を利用します。
パラメータは、以下を参考にしてください。
https://python-visualization.github.io/folium/plugins.html
色は、降水量に応じて、青、ライム、黄、オレンジ、赤に変化します。

中心点は広島にしました。

import folium
from folium import plugins

# アメダス地点を地図に表示
# 中心点設定
sta_no = 67437 # 広島
latitude = amedas_df.loc[int(sta_no), '緯度']
longitude = amedas_df.loc[int(sta_no), '経度']
map = folium.Map(location=[latitude, longitude], zoom_start=7)

# ヒートマップ表示
hm = plugins.HeatMapWithTime(data, index, radius=100, min_opacity=0, max_opacity=1, auto_play=True,
                             gradient={0.1: 'blue', 0.25: 'lime', 0.5:'yellow', 0.75:'orange',1.0: 'red'})
hm.add_to(map)
map

以下のように、地図上に降雨データがヒートマップ表示されます。
1時間ごとに表示が切り替わります。

foliumのヒートマップ表示上、観測地点が近い場合には、データの値が大きく見えます。
アメダス地点をベースにしているため、雨雲レーダーのように正確ではありませんが、大まかな変化は確認できると思います。
速度を10fpsなど早くするとわかりやすいかもしれません。

htmlに保存し表示することもできます。

# 地図を保存
map.save(outfile="prec20180705_08_map.html")

豪雨時の降水量を時系列ヒートマップに表示してみました。

データの公開期間が2020年3月末までのため、必要な方は、早めのダウンロードをお勧めします。