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

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')

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)

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

$ python get-jma-annual_max_daily_prcp.py

$ python plot-histogram.py

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)