【気象データ】「〇〇年に一度の大雨」の算出 (データ取得・ヒストグラム作成編)【統計解析】

【気象データ】「〇〇年に一度の大雨」の算出 (当てはめる分布関数の解説編)【統計解析】

はじめに


ここ最近、大雨のニュースを見るたびに「〇〇年に一度の大雨」という表現を目にします。これは「千年に一人の美少女」や「百年に一度の逸材」といったレトリックではなく、データに基づいて統計的に算出された根拠のある数字です。

気象庁や自治体では観測結果をもとにこのような量を算出し、発表することで防災に役立てています。このシリーズでは、算出方法を確認し、Pythonを使って、実装・計算してみたいと思います。

なお、ここで紹介する手法は、実際に公的機関が発表しているものとは異なる場合が多分にありますので、十分にご注意ください。

動作検証済み環境

macOS Monterey 12.6.2, Python 3.9.15, MetPy 1.3.1

「○○年に一度の大雨」の計算方法とは?


気象庁は観測結果をもとに、各地点について「〇〇年に一度の降水量」を算出し、発表しています。

ある期間内に1回起こると考えられる降水量のことを確率降水量といいます。気象庁の解説ページでは、確率降水量の推定方法が解説されていますので、これに基づいて「〇〇年に一度の大雨」の値を計算してみたいと思います。

大まかには次のような流れになっています。

  1. 年最大日降水量のヒストグラムを作成する
  2. 分布関数を当てはめる
  3. 分布関数の当てはまり具合を確認する
  4. 当てはめた分布関数から確率降水量を算出する

今回は、年日最大降水量 (各年について、24時間に降った雨の最大値のデータ) を取得し、様々なプロットを通してデータを把握したいと思います (「1. 年最大日降水量のヒストグラムを作成する」まで)。

0. データを取得する

今回使用する年最大日降水量は、気象庁の過去の気象データ検索ページから取得することができます。

選択してコピペすることもできますが、非常に使いづらいです。以前の記事で紹介したウェブスクレイピングのPythonのコードを応用して、表データを取得することができます。

年最大日降水量をcsvファイルとして保存するコードを用意したのでお使いください。

1. ヒストグラムを作成する

ヒストグラム (histogram) とは、データの分布を視覚的に把握するためのグラフです。データをある数値の幅を持った階級 (bin) に分割して、その階級内に入るデータの個数を縦軸に表したものです。

また、累積ヒストグラム (cummulative histogram) というものも考えることができます。これは、縦軸にある階級までの累積個数を表したものです (したがって、右上がりの単調増加のグラフになります) 。

データが用意できたら、

  • 横軸を「西暦」、縦軸を「年最大日降水量」とする棒グラフ、
  • 年最大日降水量のヒストグラム
  • 年最大日降水量の累積ヒストグラム

を作成してみましょう。これによって、年最大日降水量の大まかな傾向を把握することができます。

データの取得とヒストグラムの作成


データの取得

まず、気象庁の過去の気象データ検索ページにアクセスし、データを取得するプログラムを以下に示します。

import urllib.request
from bs4 import BeautifulSoup

#---------------------------------------------------------
# 過去の気象データのページにアクセスする
url = f"<https://www.data.jma.go.jp/obd/stats/etrn/view/annually_s.php?prec_no=62&block_no=47772&year=&month=&day=&view=p1>" 
html = urllib.request.urlopen(url).read()

# BeautifulSoupでHTMLを解析
soup = BeautifulSoup(html, "html.parser")

# 空リストを用意
year = [] # 西暦用
annual_max_daily_prcp = [] # 年最大日降水量用

#------------------
# 指定高度面のデータ id='tablefix2'の<table>を抽出
table = soup.find('table', id='tablefix1')

# 各行を解析する
for i , tr in enumerate(table.find_all('tr')):
	# 3行目以降がデータ
	if i >= 3:
		data = tr.find_all('td')
		year.append(data[0].text.strip())
		annual_max_daily_prcp.append(data[4].text.strip())

#------------------
# データの書き出し csvデータとして保存
fout_path = f'data.csv'
with open(fout_path, 'w') as fout:
	# ヘッダ
	fout.write(f'# downloaded from: {url}\\n')
	
	fout.write(f'year,annual maximum daily precipitation (mm)\\n')

	# tablefix1のデータ
	for i in range(len(year)):
		fout.write(f'{year[i]},{annual_max_daily_prcp[i]}\\n')

使用の際には、6行目の url を好きな地点の過去のデータが表示されているページのものに変更してください。例として、大阪のページから取得するようにしています。

このプログラムを get-jma-annual_max_daily_prcp.py として、適当なディレクトリ (ここでは、~/Desktop/labcode/python/probable-precipitationとします。) に保存します。

ヒストグラムの作成

上で述べた3つのグラフを作成してみましょう。以下に実装例を示します。

import matplotlib.pyplot as plt
import numpy as np

#------------------------------------------------------------------------------
# データの読み込み
year = []
annual_max_daily_prcp = []
i = 0
with open('data.csv', 'r') as fin:
    for line in fin:
        if i > 1:
            data = line.split(',')
            year.append(int(data[0]))
            annual_max_daily_prcp.append(float(data[1].replace(']', '')))
        i += 1

# ndarrayに変換
year = np.array(year)
annual_max_daily_prcp = np.array(annual_max_daily_prcp)

#------------------------------------------------------------------------------
# --- プロット
# 全般の設定
fig = plt.figure(figsize=(210/25.4, 294/25.4), dpi=100)
out_fig_path = 'histogram.png'

# プロット枠の設定
ax01 = fig.add_axes([0.10, 0.70, 0.80, 0.20])
ax02 = fig.add_axes([0.10, 0.40, 0.80, 0.20])
ax03 = fig.add_axes([0.10, 0.10, 0.80, 0.20])

# --- plot 1
# 年最大日降水量を西暦毎にバーグラフで表示
ax01.bar(year, annual_max_daily_prcp)

ax01.set_xlabel('year')
ax01.set_ylabel('annual maximum\\ndaily precipitation (mm)')

# グラフ上に最大値、最小値、平均値を表示
fig.text(0.01, 0.9, f'max: {annual_max_daily_prcp.max()} mm ({year[np.argmax(annual_max_daily_prcp)]})', transform=ax01.transAxes)
fig.text(0.01, 0.8, f'min: {annual_max_daily_prcp.min()} mm ({year[np.argmin(annual_max_daily_prcp)]})', transform=ax01.transAxes)
fig.text(0.01, 0.7, f'mean: {annual_max_daily_prcp.mean():.1f} mm', transform=ax01.transAxes)

# --- plot 2 と plot 3
# bin のためのリストを作成する
bins = [10*(i+1) for i in range(int(annual_max_daily_prcp.max()/10) + 1)]

# ヒストグラムを表示
ax02.hist(annual_max_daily_prcp, bins=bins, density=True)

# 累積ヒストグラムを表示
ax03.hist(annual_max_daily_prcp, bins=bins, cumulative=True, density=True)
ax03.set_xlabel('annual maximum daily precipitation (mm)')

#------------------------------------------------------------------------------
# --- 図の保存
plt.savefig(out_fig_path, transparent=False)

このコードを plot-histogram.pyとして get-jma-annual_max_daily_prcp.py を保存したのと同じ適当なディレクトリ (ここでは、~/Desktop/labcode/python/probable-precipitationとします。) に保存します。

プログラムを実行する

ターミナルを開き、

$ cd ~/Desktop/labcode/python/probable-precipitation

と入力し、ディレクトリを移動します( $マークは無視してください)。まず、データを取得しましょう。

$ python get-jma-annual_max_daily_prcp.py

とします。data.csv というファイルが、現在のディレクトリ (つまり、~/Desktop/labcode/python/probable-precipitation ) にできているはずです。内容は、西暦とその年の年最大日降水量がコンマ区切りで2列に並んだデータです。

つぎに、

$ python plot-histogram.py

とします。現在のディレクトリに histogram.pngというファイルができて、次のようになっていれば成功です。

一番上は入力データの1883年から2023年までの年最大日降水量が棒グラフとして表示されています。最大値が1957年の250.7 mm、最小値が1981年の42.0 mm、平均値が93.0 mmであることも表示されています。

二番目はヒストグラムで、bin幅 (階級幅)を10 mmとして、正規化して (すべて足すと1になるようにしてある) 表示されています。60 mm から70 mm の bin に入るデータの個数と 70 mmから80 mmに入るデータの個数が同数で最も多いことや、右に裾野が広い分布であることがわかります。また、250 mmを超えたイベントががいかに大きかったかがわかります。

一番下は累積ヒストグラムです。同じく、bin幅を10 mmとして、正規化して表示されています。累積ヒストグラムはヒストグラムを小さい方から積み上げていったものです。ほぼ、統計を取った年のうち、80%程度の年が年最大日降水量が120 mm以下に収まることがわかります。

コードの解説


ヒストグラムを作成したプログラムのソースコードの解説をしていきます。データ取得のプログラムについては以前の記事を参照してください。

import matplotlib.pyplot as plt
import numpy as np

まず、必要なモジュールをインポートします。

#------------------------------------------------------------------------------
# データの読み込み
year = []
annual_max_daily_prcp = []
i = 0
with open('data.csv', 'r') as fin:
    for line in fin:
        if i > 1:
            data = line.split(',')
            year.append(int(data[0]))
            annual_max_daily_prcp.append(float(data[1].replace(']', '')))
        i += 1

# ndarrayに変換
year = np.array(year)
annual_max_daily_prcp = np.array(annual_max_daily_prcp)

get-jma-annual_max_daily_prcp.pyが作成したデータファイルを読み込みます。データには「 ] 」という記号が入っていますが、これは、そのデータがそれ以上であることを示すものです (気象庁「値欄記号の説明」参照) 。今回 (2023/08/21執筆時現在) は2023年のデータが135.5 ] となっていますが、これは2023年がまだ続いているからです。普通はこのようなデータは使いませんが、今回は含めるようにしています。

#------------------------------------------------------------------------------
# --- プロット
# 全般の設定
fig = plt.figure(figsize=(210/25.4, 294/25.4), dpi=100)
out_fig_path = 'histogram.png'

# プロット枠の設定
ax01 = fig.add_axes([0.10, 0.70, 0.80, 0.20])
ax02 = fig.add_axes([0.10, 0.40, 0.80, 0.20])
ax03 = fig.add_axes([0.10, 0.10, 0.80, 0.20])

# --- plot 1
# 年最大日降水量を西暦毎にバーグラフで表示
ax01.bar(year, annual_max_daily_prcp)

ax01.set_xlabel('year')
ax01.set_ylabel('annual maximum\\ndaily precipitation (mm)')

# グラフ上に最大値、最小値、平均値を表示
fig.text(0.01, 0.9, f'max: {annual_max_daily_prcp.max()} mm ({year[np.argmax(annual_max_daily_prcp)]})', transform=ax01.transAxes)
fig.text(0.01, 0.8, f'min: {annual_max_daily_prcp.min()} mm ({year[np.argmin(annual_max_daily_prcp)]})', transform=ax01.transAxes)
fig.text(0.01, 0.7, f'mean: {annual_max_daily_prcp.mean():.1f} mm', transform=ax01.transAxes)

# --- plot 2 と plot 3
# bin のためのリストを作成する
bins = [10*(i+1) for i in range(int(annual_max_daily_prcp.max()/10) + 1)]

# ヒストグラムを表示
ax02.hist(annual_max_daily_prcp, bins=bins, density=True)

# 累積ヒストグラムを表示
ax03.hist(annual_max_daily_prcp, bins=bins, cumulative=True, density=True)
ax03.set_xlabel('annual maximum daily precipitation (mm)')

この部分はプロットをしています。

棒グラフは .bar(x, height, width=0.8, bottom=None, *, align='center', data=None, **kwargs) でできます。xyearheightannual_max_daily_prcp を入れています。

グラフ上に文字を表示するには、..text(x, y, s, fontdict=None, withdash=<deprecated parameter>, **kwargs) を使います。transform=ax01.transAxes とすることで、xy で指定した位置はグラフ上の座標を指定するもの (グラフ上で左下が原点 (0, 0) で右上が (1, 1) ) とみなされます。

bins は 階級 (bin) を指定するためのリストです。0mm 以上 10 mm 未満、10 mm 以上 20 mm 未満、20 mm 以上 30mm 未満、… というように10 mm刻みで階級を作ります。階級の区切りの値のリストを指定すればよいです。したがって、0から年最大日降水量の最大値が含まれるような数のリストを作成しています。

ヒストグラムを作成するには.hist(x, bins=None, range=None, density=None, weights=None, cumulative=False, bottom=None, histtype='bar', align='mid', orientation='vertical', rwidth=None, log=False, color=None, label=None, stacked=False, normed=None, *, data=None, **kwargs) を使います。ここでは、基本的な使い方として、binsを指定し、densityTrueとして、正規化しています。

bins を指定しない場合、データの範囲を10分割したものになるようです。bin幅は上で述べたように恣意的に決めていますが、データの個数等から最適と思われるbin数や幅を決定する手法がいくつか知られています。サポートされているものとしては、numpyのページをご覧ください。

また、cumulativeTrueにすると、三番目で示したように、累積ヒストグラムが得られます。

#------------------------------------------------------------------------------
# --- 図の保存
plt.savefig(out_fig_path, transparent=False)

最後に図を保存して終了します。

最後に


今回は、「〇〇年に一度の大雨」について理解するために、その基礎資料となる年最大日降水量のデータ取得と、ヒストグラムの作成方法について紹介しました。

次回は、このヒストグラムに分布を当てはめる方法や当てはめるべき分布の種類について紹介します。

参考文献


本記事を執筆するに当たり、次の資料・文献を参考にしました。

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です