""" L05: 14日豪雨ヒートマップ — 2024 広島豪雨の空間×時間構造を読み解く 2024-06-29 〜 2024-07-12 の 14 日分・約 280 観測所の 10 分雨量を DoBoX #1275 から取得し、 (時刻 × 観測所) → 日合計 → (観測所 × 日) のマトリクス に再構成して、空間 (観測所/事務所地域) × 時間 (14日) の複合可視化で 「広島豪雨」の構造を多角的に読み解く v2 教材。 実行: cd "2026 DoBoX 教材" py -X utf8 lessons/L05_14days_heavy_rain.py データは存在しなければ DoBoX から自動取得する (ensure_dataset)。 """ from pathlib import Path import pandas as pd import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.colors import LinearSegmentedColormap from _common import (ROOT, ASSETS, LESSONS, render_lesson, code, figure, data_recipe, parse_rain_csv, ensure_dataset) plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False # === 0. データ自動取得 (14日分 × 個別 resource_id) ============================ RAIN_DIR = ROOT / "data" / "rain_2024" RAIN_DIR.mkdir(parents=True, exist_ok=True) RA_PATH = ROOT / "data" / "extras" / "rainfall_annual.csv" # 観測所→事務所のマップ用 RAIN_RID = { "2024-06-29": 94490, "2024-06-30": 94495, "2024-07-01": 94500, "2024-07-02": 94505, "2024-07-03": 94510, "2024-07-04": 94515, "2024-07-05": 94520, "2024-07-06": 94525, "2024-07-07": 94530, "2024-07-08": 94535, "2024-07-09": 94540, "2024-07-10": 94545, "2024-07-11": 94551, "2024-07-12": 94556, } print("=== 0. データ自動取得 ===") for date, rid in RAIN_RID.items(): ensure_dataset(RAIN_DIR / f"rain_{date}.csv", resource_id=rid, min_bytes=500_000, label=f"雨量10分値 {date}") ensure_dataset(RA_PATH, dataset_id=1276, label="雨量年集計 #1276 (事務所メタ)") # === 1. 14日分を縦に連結 → 観測所×日 マトリクスを作る ======================== print("\n=== 1. 14日分を連結 → 観測所×日 マトリクス ===") files = sorted(RAIN_DIR.glob("rain_2024-*.csv")) print(f"ファイル数: {len(files)}") daily_series = [] for f in files: tidy = parse_rain_csv(f) # (時刻 × 観測所) の10分値 daily = tidy.sum(axis=0, min_count=1) # 観測所ごとの日合計 (Series) daily.name = pd.Timestamp(f.stem.replace("rain_", "")) daily_series.append(daily) # 観測所 × 日 マトリクス (NaN は欠測) M = pd.concat(daily_series, axis=1).sort_index(axis=1) M.columns = [c.strftime("%m-%d") for c in M.columns] print(f"マトリクス shape: {M.shape} (観測所 × 14日)") print(f"非欠測セル: {M.notna().sum().sum()} / {M.size}") print(f"全体最大値: {M.max().max():.1f} mm/日") # 全 14 日合計と日数有効値で観測所をランク付け M_filled = M.fillna(0.0) station_total = M_filled.sum(axis=1).sort_values(ascending=False) print(f"観測所 14日合計の上位5件:\n{station_total.head().to_string()}") # === 2. 観測所→事務所メタの構築 (空間グルーピング用) ========================= print("\n=== 2. 観測所→事務所(地域)メタ の構築 ===") ra_raw = pd.read_csv(RA_PATH, header=None, encoding="utf-8-sig") station_office = pd.DataFrame({ "office": ra_raw.iloc[0, 1:].astype(str).tolist(), "station": ra_raw.iloc[4, 1:].astype(str).tolist(), }) # 重複名は最初の値を採用 office_map = {} for _, r in station_office.iterrows(): s, o = r["station"], r["office"] if s not in office_map and s not in ("nan", ""): office_map[s] = o print(f"観測所→事務所 map size: {len(office_map)}") def lookup_office(col): base = col.split("_")[0] if "_" in col else col return office_map.get(col, office_map.get(base, "(不明)")) M_office = pd.Series({c: lookup_office(c) for c in M.index}, name="office") office_count = M_office.value_counts() print(f"事務所別観測所数:\n{office_count.to_string()}") # === 3. 図1: 全観測所×日 ヒートマップ (合計値で並べ替え) ===================== print("\n=== 3. 図1: 全観測所×日 ヒートマップ ===") M_sorted = M_filled.loc[station_total.index] # 上から下へ「降水量の多い観測所順」 fig, ax = plt.subplots(figsize=(11, 13)) # YlOrRd を白背景にした方が情報量が増える cmap = plt.get_cmap("YlOrRd").copy() cmap.set_bad("#eeeeee") disp = M.loc[station_total.index] # NaNは灰色に im = ax.imshow(disp.values, aspect="auto", cmap=cmap, vmin=0, vmax=np.nanpercentile(M.values, 99.5), interpolation="nearest") ax.set_xticks(range(M.shape[1])) ax.set_xticklabels(M.columns, rotation=45, ha="right", fontsize=9) ax.set_yticks([]) # 280観測所はラベル不可、補助線で区切る ax.set_ylabel(f"観測所 (上=14日合計の多い順, n={M.shape[0]})") ax.set_xlabel("日付 (2024年)") ax.set_title("広島県内 全観測所 × 14日 日合計雨量 ヒートマップ\n" "上端ほど14日合計が多い観測所。縦に走る濃い帯=県全域に降った日") cbar = plt.colorbar(im, ax=ax, shrink=0.6, label="日合計雨量 (mm)") # 上位10観測所は名前を出す top10 = station_total.head(10) ypos = np.arange(len(top10)) for i, (s, v) in enumerate(top10.items()): ax.text(-0.6, i, f"{s}", ha="right", va="center", fontsize=8, color="#222") plt.tight_layout() plt.savefig(ASSETS / "L05_heatmap_all.png", dpi=140) plt.close() # === 4. 図2: 上位10観測所 14日重ね描き (small multiples 2x5) ================ print("\n=== 4. 図2: 上位10観測所 small multiples ===") top10_names = station_total.head(10).index.tolist() fig, axes = plt.subplots(2, 5, figsize=(15, 6.2), sharex=True, sharey=True) xs = pd.to_datetime([f"2024-{c}" for c in M.columns]) ymax = M_filled.loc[top10_names].max().max() * 1.08 for ax, name in zip(axes.ravel(), top10_names): series = M.loc[name] bars = ax.bar(xs, series.fillna(0).values, width=0.85, color="#0969da", alpha=0.7, edgecolor="#0a3a72", linewidth=0.4) # ピーク日を赤で強調 peak_i = int(np.nanargmax(series.values)) bars[peak_i].set_color("#cf222e") bars[peak_i].set_alpha(0.95) total = station_total[name] peak_v = series.iloc[peak_i] peak_d = M.columns[peak_i] office = lookup_office(name) ax.set_title(f"{name} ({office})\n14日計 {total:.0f}mm " f"ピーク {peak_d} {peak_v:.0f}mm", fontsize=10) ax.tick_params(axis="x", rotation=45, labelsize=8) ax.set_ylim(0, ymax) ax.grid(alpha=0.25, axis="y") ax.set_axisbelow(True) for ax in axes[:, 0]: ax.set_ylabel("日合計雨量 (mm)") fig.suptitle("14日合計 上位10観測所の 日次雨量 (赤バー = 各観測所のピーク日)", fontsize=13, y=1.0) plt.tight_layout() plt.savefig(ASSETS / "L05_top10_smallmultiples.png", dpi=140) plt.close() # === 5. 図3: 事務所(地域) × 日 ヒートマップ (空間集約マップ) ================= print("\n=== 5. 図3: 事務所×日 ヒートマップ (地域集約) ===") # 事務所単位で平均雨量 (mm/日/観測所) を集計 office_daily_rows = [] for office, group in pd.Series(M_office).groupby(M_office): cols = group.index series = M_filled.loc[cols].mean(axis=0) series.name = f"{office} (n={len(cols)})" office_daily_rows.append(series) office_M = pd.DataFrame(office_daily_rows) # 14日合計でソート (多い→少ない) office_M = office_M.loc[office_M.sum(axis=1).sort_values(ascending=False).index] fig, ax = plt.subplots(figsize=(11, 5.2)) im = ax.imshow(office_M.values, aspect="auto", cmap="YlOrRd", vmin=0, vmax=np.nanpercentile(office_M.values, 99)) ax.set_xticks(range(office_M.shape[1])) ax.set_xticklabels(office_M.columns, rotation=45, ha="right", fontsize=9) ax.set_yticks(range(len(office_M))) ax.set_yticklabels(office_M.index, fontsize=10) ax.set_xlabel("日付 (2024年)") ax.set_title("事務所(=地域)別 平均日雨量 14日推移 — 空間集約マップ\n" "セル数値: 事務所内観測所の平均日合計 (mm)") for i in range(office_M.shape[0]): for j in range(office_M.shape[1]): v = office_M.values[i, j] if not np.isnan(v): txt_color = "white" if v > office_M.values.max() * 0.55 else "black" ax.text(j, i, f"{v:.0f}", ha="center", va="center", color=txt_color, fontsize=8) plt.colorbar(im, ax=ax, label="平均日合計 (mm/観測所)", shrink=0.85) plt.tight_layout() plt.savefig(ASSETS / "L05_office_heatmap.png", dpi=140) plt.close() # === 6. 図4: 日別 県全体雨量ピラミッド ======================================= print("\n=== 6. 図4: 日別 県全体雨量ピラミッド ===") # 日別の代表値: 県内最大, 上位5観測所平均, 中央値, 観測所数 day_stats = pd.DataFrame({ "max": M_filled.max(axis=0), "top5": M_filled.apply(lambda s: s.nlargest(5).mean(), axis=0), "p90": M_filled.quantile(0.90, axis=0), "median": M_filled.median(axis=0), "active_stations": (M_filled > 0).sum(axis=0), }) print(day_stats.round(1).to_string()) fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(11, 7), gridspec_kw={"height_ratios": [2.2, 1]}, sharex=True) xs = pd.to_datetime([f"2024-{c}" for c in M.columns]) # 上: 4 段の代表値 (棒+線) ax1.bar(xs, day_stats["max"].values, width=0.85, color="#cf222e", alpha=0.85, label="県内最大 (1観測所)") ax1.plot(xs, day_stats["top5"].values, "o-", color="#fb8500", linewidth=2, markersize=6, label="上位5観測所平均") ax1.plot(xs, day_stats["p90"].values, "s-", color="#1f883d", linewidth=1.5, markersize=5, label="90パーセンタイル") ax1.plot(xs, day_stats["median"].values, "^-", color="#0969da", linewidth=1.5, markersize=5, label="中央値") ax1.set_ylabel("日合計雨量 (mm)") ax1.set_title("日別 県全体 雨量ピラミッド — 代表値の使い分け") ax1.legend(loc="upper left", fontsize=9) ax1.grid(alpha=0.3, axis="y") # 下: 雨が降った観測所数 (active station count) — 雨域の広がり ax2.bar(xs, day_stats["active_stations"].values, width=0.85, color="#8250df", alpha=0.75) ax2.axhline(M.shape[0], color="gray", linewidth=0.6, linestyle="--", label=f"全観測所数 n={M.shape[0]}") ax2.set_ylabel("有雨観測所数") ax2.set_xlabel("日付 (2024年)") ax2.set_title("雨域の広がり (>0mm の観測所数)") ax2.legend(loc="upper left", fontsize=9) ax2.grid(alpha=0.3, axis="y") ax2.tick_params(axis="x", rotation=30) plt.tight_layout() plt.savefig(ASSETS / "L05_daily_pyramid.png", dpi=140) plt.close() # === 7. 図5: 観測所別 14日合計 ECDF (累積分布関数) ========================== print("\n=== 7. 図5: 観測所別 14日合計 ECDF ===") totals = station_total.values xs_e = np.sort(totals) ys_e = np.arange(1, len(xs_e) + 1) / len(xs_e) fig, ax = plt.subplots(figsize=(9, 5.2)) ax.plot(xs_e, ys_e, "-", color="#0969da", linewidth=2) ax.fill_between(xs_e, 0, ys_e, color="#0969da", alpha=0.10) # 代表的な分位点をマーク for q, color, label in [(0.50, "#1f883d", "中央値 (50%)"), (0.90, "#fb8500", "90%点"), (0.99, "#cf222e", "99%点")]: qv = np.quantile(totals, q) ax.axvline(qv, color=color, linestyle="--", linewidth=1.2) ax.text(qv, 0.04, f" {label}\n {qv:.0f}mm", color=color, fontsize=9, ha="left", va="bottom") # 上位3観測所を散布点で top3 = station_total.head(3) for i, (name, v) in enumerate(top3.items()): rank = (totals >= v).sum() ax.plot(v, 1 - (rank - 1) / len(totals), "o", color="#cf222e", markersize=10, markeredgecolor="black", markeredgewidth=1) ax.annotate(f"{name}\n{v:.0f}mm", xy=(v, 1 - (rank - 1) / len(totals)), xytext=(-50, -22 - i*14), textcoords="offset points", fontsize=9, color="#cf222e", arrowprops=dict(arrowstyle="->", color="#cf222e", lw=0.8)) ax.set_xlabel("14日合計雨量 (mm)") ax.set_ylabel("累積確率 F(x) = P(観測所 ≤ x)") ax.set_title(f"広島県 雨量観測所 14日合計 ECDF (n={len(totals)})\n" "X軸が右に伸びる「重い裾」= 一部の観測所で極端に多雨") ax.grid(alpha=0.3) ax.set_ylim(0, 1.02) plt.tight_layout() plt.savefig(ASSETS / "L05_ecdf.png", dpi=140) plt.close() # === 8. サマリ表 + マトリクスCSV出力 ======================================== print("\n=== 8. サマリ生成 ===") peak_day_idx = day_stats["max"].idxmax() peak_day_max = day_stats["max"].max() peak_station = M_filled[peak_day_idx].idxmax() total_county = M_filled.values.sum() summary_rows = [ ("対象期間", "2024-06-29 〜 2024-07-12 (14日)"), ("観測所数", f"{M.shape[0]} 観測所"), ("事務所数(地域区分)", f"{office_M.shape[0]}"), ("14日 県全体総雨量", f"{total_county:,.0f} mm (観測所×日 の合計)"), ("ピーク日(県内最大)", f"{peak_day_idx} ({peak_day_max:.0f} mm)"), ("ピーク日 最大観測所", f"{peak_station}"), ("14日合計 最多観測所", f"{station_total.index[0]} ({station_total.iloc[0]:.0f} mm)"), ("14日合計 中央値", f"{np.median(totals):.0f} mm"), ("14日合計 90%点", f"{np.quantile(totals, 0.90):.0f} mm"), ("14日合計 99%点", f"{np.quantile(totals, 0.99):.0f} mm"), ] summary_html = "" + "".join( f"" for k, v in summary_rows ) + "
{k}{v}
" print("\n".join(f" {k}: {v}" for k, v in summary_rows)) # 中間ファイル: 観測所×日 マトリクス M.to_csv(ASSETS / "rain_14d_daily_matrix.csv", encoding="utf-8-sig") print(f"saved matrix: {ASSETS / 'rain_14d_daily_matrix.csv'}") # 上位10観測所サマリ表 top10_df = pd.DataFrame({ "観測所": top10_names, "事務所": [lookup_office(n) for n in top10_names], "14日合計(mm)": [round(station_total[n], 1) for n in top10_names], "ピーク日": [M.columns[int(np.nanargmax(M.loc[n].values))] for n in top10_names], "ピーク値(mm)": [round(np.nanmax(M.loc[n].values), 1) for n in top10_names], }) top10_html = top10_df.to_html(index=False) # === 9. HTML レンダリング =================================================== print("\n=== 9. HTML レンダリング ===") sections = [ ("学習目標", """ """), ("使用データ", """

14日分の CSV は 各日が独立した 1 ファイル なので、リソース ID も日ごとに違う。 本スクリプトは下記の RAIN_RID 辞書 (14 リソース) を使って ensure_dataset() で並列に取得する。

"""), ("データ取得手順", data_recipe( [{"label": f"雨量10分値 {d}", "dataset_id": 1275, "resource_id": rid, "file": f"data/rain_2024/rain_{d}.csv", "format": "CSV (5段ヘッダ, 10分値)", "size": "1.0〜1.4 MB"} for d, rid in RAIN_RID.items()] + [{"label": "雨量年集計 (事務所メタの参照)", "dataset_id": 1276, "file": "data/extras/rainfall_annual.csv", "format": "CSV (多段ヘッダ)", "size": "約 500 KB"}] ) + """

取得スクリプト内蔵: 本レッスンの L05_14days_heavy_rain.py は、 RAIN_RID 辞書 (14 日分) を使って ensure_dataset() ヘルパで 14ファイル+1メタを順次自動DLします。 HTML を読んで PowerShell で個別取得しても、fetch_all.py で一括取得しても、結果は同じです。

"""), ("方法", """
  1. 14 リソース ID で 14 ファイル取得: RAIN_RID = {"2024-06-29": 94490, ..., "2024-07-12": 94556} を ensure_dataset(resource_id=...) でループ DL
  2. 各日CSVを tidy 化: parse_rain_csv() (年度差分対応) で (時刻 × 観測所) の 10 分値 DataFrame に
  3. 日合計に集約: tidy.sum(axis=0, min_count=1) で観測所ごとの日合計 Series
  4. 14日分を横結合: pd.concat([...], axis=1) で (観測所 × 日) マトリクス M
  5. 事務所メタの結合: rainfall_annual.csv の 0 行目(事務所名) と 4 行目(観測所名) から {station: office} 辞書を構築
  6. 5 種類の可視化: 全観測所ヒートマップ / 上位10 small multiples / 事務所×日 集約ヒートマップ / 日別代表値ピラミッド+雨域広がり / 14日合計 ECDF
  7. サマリ表: 期間/観測所数/ピーク日/最大観測所/分位点 を 1 表に集約
"""), ("コード解説", code(''' from pathlib import Path import pandas as pd, numpy as np import matplotlib.pyplot as plt from _common import parse_rain_csv, ensure_dataset # === (0) 14リソース ID で 14日分を順次DL === RAIN_RID = { "2024-06-29": 94490, "2024-06-30": 94495, "2024-07-01": 94500, "2024-07-02": 94505, "2024-07-03": 94510, "2024-07-04": 94515, "2024-07-05": 94520, "2024-07-06": 94525, "2024-07-07": 94530, "2024-07-08": 94535, "2024-07-09": 94540, "2024-07-10": 94545, "2024-07-11": 94551, "2024-07-12": 94556, } for date, rid in RAIN_RID.items(): ensure_dataset(f"data/rain_2024/rain_{date}.csv", resource_id=rid, min_bytes=500_000) # === (1) 14日分を 観測所×日 マトリクスに === files = sorted(Path("data/rain_2024").glob("rain_2024-*.csv")) daily_series = [] for f in files: tidy = parse_rain_csv(f) # (時刻 × 観測所) 10分値 daily = tidy.sum(axis=0, min_count=1) # 観測所ごとの日合計 daily.name = pd.Timestamp(f.stem.replace("rain_", "")) daily_series.append(daily) M = pd.concat(daily_series, axis=1) # (観測所 × 14日) M_filled = M.fillna(0.0) # === (2) 観測所→事務所 マップ (空間グルーピング用) === ra = pd.read_csv("data/extras/rainfall_annual.csv", header=None, encoding="utf-8-sig") office_map = dict(zip(ra.iloc[4, 1:].astype(str), # 観測所名 ra.iloc[0, 1:].astype(str))) # 事務所名 # === (3) 14日合計で観測所をランク付け → ヒートマップ === station_total = M_filled.sum(axis=1).sort_values(ascending=False) disp = M.loc[station_total.index] # 上=多雨, 下=少雨 plt.imshow(disp.values, aspect="auto", cmap="YlOrRd", vmin=0, vmax=np.nanpercentile(M.values, 99.5)) # === (4) 上位10観測所 small multiples (2x5) === top10 = station_total.head(10).index fig, axes = plt.subplots(2, 5, figsize=(15, 6.2), sharex=True, sharey=True) for ax, name in zip(axes.ravel(), top10): ax.bar(range(14), M.loc[name].fillna(0).values, color="#0969da") # === (5) 事務所×日 集約ヒートマップ (空間集約) === M_office = pd.Series({c: office_map.get(c, "(不明)") for c in M.index}) office_M = (M_filled.groupby(M_office.values).mean()) # 事務所平均/日 # === (6) 日別 県全体 ピラミッド (代表値の使い分け) === day_stats = pd.DataFrame({ "max": M_filled.max(axis=0), "top5": M_filled.apply(lambda s: s.nlargest(5).mean(), axis=0), "p90": M_filled.quantile(0.90, axis=0), "median": M_filled.median(axis=0), }) # === (7) 観測所別14日合計 ECDF === totals = station_total.values xs = np.sort(totals) ys = np.arange(1, len(xs)+1) / len(xs) plt.plot(xs, ys) # 重い裾を見る plt.axvline(np.quantile(totals, 0.99)) # 99%点 ''')), ("結果", figure("assets/L05_heatmap_all.png", f"全観測所(n={M.shape[0]}) × 14日 ヒートマップ。縦軸は14日合計の多い順。" "縦に走る濃い帯=県全域に降った日、点状の濃さ=局所豪雨を示す") + figure("assets/L05_top10_smallmultiples.png", "14日合計上位10観測所の日次バー (赤=各観測所のピーク日)。" "ピーク日が06-30〜07-01に集中する観測所と07-08〜07-10にずれる観測所がある") + figure("assets/L05_office_heatmap.png", "事務所(=地域区分)×日 の平均日雨量ヒートマップ。" "観測所280点を10程度の事務所に集約することで、雨域の南北・東西の動きが読める") + figure("assets/L05_daily_pyramid.png", "上: 日別の代表値ピラミッド (赤バー=県内最大, 橙=上位5平均, 緑=90%点, 青=中央値)。" "下: 有雨観測所数 (>0mm)。同じ日でも代表値ごとに様相が違う") + figure("assets/L05_ecdf.png", "14日合計の経験累積分布。緑=中央値・橙=90%点・赤=99%点。" "右へ伸びる重い裾=ごく一部の観測所が極端に多雨") + "

期間サマリ

" + summary_html + "

14日合計 上位10観測所

" + top10_html), ("考察", """ """), ("発展課題", """
  1. 10 分粒度のままヒートマップを描く: 日合計に潰さず、(観測所 × 144スロット × 14日) の 3 次元配列を時間軸に並べた長尺ヒートマップで 線状降水帯の通過を可視化する。matplotlib の imshow でも 280 × 2016 = 56 万セル程度なら十分扱える。
  2. 事務所地域 → 緯度経度マップ: 観測所の位置情報 (DoBoX #1275 の観測所諸元 or 国交省 XRAIN メタ) を結合し、folium で 14日合計の色付き円マーカー をプロットする。本レッスンは事務所単位の集約で代用したが、座標があれば真の空間ヒートマップに格上げできる。
  3. 豪雨日の「同時性」検定: 2024-07-01 と 2024-07-10 で多雨の観測所がどれだけ重なるか、Jaccard 係数や Spearman 順位相関で評価。広域前線型と局所型で「上位観測所の入れ替わり」がどれだけ起きるかを定量化する。
  4. ECDF の片対数プロット → 裾指数推定: ECDF の生存関数 1-F(x) を log-log でプロットし、裾部分が直線にのるか検査。直線なら パレート分布、降雨の極値統計 (GEV, GPD) との接続。L210 の極値解析と接続する。
  5. 過去年度との比較: 同じ DoBoX #1275 で 2018年豪雨期間 (2018-07-04 〜 07-08) のリソースを取得し、2024年豪雨と比較する。観測所数・地点・代表値の違いを直接議論できる (時系列の長期トレンド観察)。
"""), ] html = render_lesson( num=5, title="14日豪雨ヒートマップ — 2024 広島豪雨の空間×時間構造", tags=["v2-rewrite", "基礎", "選択", "時系列", "空間集約", "ヒートマップ"], time="100分", level="リテラシ基礎〜応用基礎", data_label='DoBoX #1275 雨量10分値 (14日 × 約280観測所) ' '+ #1276 (事務所メタのみ)', sections=sections, script_filename="lessons/L05_14days_heavy_rain.py", ) out = LESSONS / "L05_14days_heavy_rain.html" out.write_text(html, encoding="utf-8") print(f"\nsaved: {out}") print(f" PNGs: 5 generated in {ASSETS}")