""" L04 (v2): 雨量10分値CSV — データ前処理の解剖 DoBoX #1275 の雨量10分値CSV (5段ヘッダ × 1605列, 401観測所 × 4列構造) を tidy data に整形する過程を「解剖」する教材。単に整形するだけでなく: (1) ヘッダ位置の動的検出ロジック (年度違いに堅牢) (2) 観測所メタの抽出 (事務所/水系/河川 の階層構造) (3) 完備性: 観測所×時刻 の欠測ヒートマップ (4) フラグ列 (0x...) の bit パターン解読 (5) サンプル観測所の10分時系列重ね描き (6) 整形品質指標 (raw shape → tidy shape, 欠測率, ユニーク観測所数) 実行: cd "2026 DoBoX 教材" py -X utf8 lessons/L04_tidy_rainfall.py """ from pathlib import Path import re import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap 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. データ自動取得 ======================================================= SRC = ROOT / "data" / "rain_2024" / "rain_2024-07-01.csv" print("=== 0. データ自動取得 ===") ensure_dataset(SRC, resource_id=94500, min_bytes=500_000, label="雨量10分値 2024-07-01") # === 1. 生CSVを読み込み — header=None でフラットに ============================ print("\n=== 1. 生CSVを header=None で全部「データ」として読込 ===") raw = pd.read_csv(SRC, header=None, encoding="utf-8-sig") print(f"raw shape: {raw.shape} (行 × 列)") print(f" → 1605 列 = 401 観測所 × 4 列 (10分雨量/フラグ/累計/フラグ) + 時刻列") # === 2. ヘッダ位置の動的検出 ================================================== print("\n=== 2. ヘッダ位置の動的検出 ===") # step-by-step: 各行の "10分雨量" 出現件数を数える header_scan = [] for r in range(min(15, len(raw))): row_str = raw.iloc[r, :].astype(str).tolist() n_10min = sum("10分雨量" in v for v in row_str) n_obsname = (str(raw.iloc[r, 0]).strip() == "観測所名") n_riverlbl = (str(raw.iloc[r, 0]).strip() == "河川名") header_scan.append({ "row": r, "first_cell": str(raw.iloc[r, 0])[:18], "n_10分雨量": n_10min, "is_観測所名行": n_obsname, "is_河川名行": n_riverlbl, }) scan_df = pd.DataFrame(header_scan) print(scan_df.to_string(index=False)) header_idx = next(r for r in range(15) if any("10分雨量" in str(v) for v in raw.iloc[r, :])) obs_idx = next(r for r in range(header_idx) if str(raw.iloc[r, 0]).strip() == "観測所名") print(f"\n→ ヘッダ行 (10分雨量を含む): row={header_idx}") print(f"→ 観測所名行: row={obs_idx}") print(f"→ メタ行 (上から順): {[str(raw.iloc[r, 0]) for r in range(obs_idx + 1)]}") # === 3. tidy DataFrame の生成 ================================================ print("\n=== 3. tidy DataFrame (観測所×時刻) を生成 ===") tidy = parse_rain_csv(SRC) print(f"tidy shape: {tidy.shape} (時刻 × 観測所)") print(f" raw 列数 1605 → tidy 列数 {tidy.shape[1]} (≒ 1605/4)") print(f" raw データ行 144 → tidy 行数 {tidy.shape[0]} (10分×24時間=144)") miss_rate = tidy.isna().mean().mean() * 100 print(f" 全セルの欠測率: {miss_rate:.2f}%") print(f" ユニーク観測所名: {tidy.columns.nunique()} (重複名は連番化)") # === 4. 観測所メタ抽出 (事務所/水系/河川の階層) ================================ print("\n=== 4. 観測所メタ抽出 ===") meta_rows = { "office": raw.iloc[0, 1:].astype(str).tolist(), # 事務所名 "owner": raw.iloc[1, 1:].astype(str).tolist(), # データ所管 "system": raw.iloc[2, 1:].astype(str).tolist(), # 水系名 "river": raw.iloc[3, 1:].astype(str).tolist(), # 河川名 "sid": raw.iloc[4, 1:].astype(str).tolist(), # 観測所番号 "station": raw.iloc[5, 1:].astype(str).tolist(), # 観測所名 } # 観測所ごとに 4 列繰り返しているので、最初の 1 列だけ拾う (10分雨量列) cols_10min = [i for i, v in enumerate(raw.iloc[header_idx, :]) if "10分雨量" in str(v)] meta = pd.DataFrame({k: [v[c - 1] for c in cols_10min] for k, v in meta_rows.items()}) # NaN 文字列を整理 for c in meta.columns: meta[c] = meta[c].replace({"nan": "(不明)"}) print(f"meta shape: {meta.shape}") print("水系×事務所 観測所数:") ws_off = meta.groupby(["system", "office"]).size().unstack(fill_value=0) print(ws_off.to_string()) # === 5. 図1: 観測所階層構造のヒートマップ ==================================== print("\n=== 5. 図1: 水系×事務所 観測所数ヒートマップ ===") top_systems = meta["system"].value_counts().head(15).index.tolist() top_offices = meta["office"].value_counts().head(8).index.tolist() hm = (meta[meta["system"].isin(top_systems) & meta["office"].isin(top_offices)] .groupby(["system", "office"]).size().unstack(fill_value=0) .reindex(index=top_systems, columns=top_offices, fill_value=0)) fig, ax = plt.subplots(figsize=(9, 0.45 * len(hm) + 1.5)) im = ax.imshow(hm.values, cmap="YlOrRd", aspect="auto") ax.set_xticks(range(len(hm.columns))) ax.set_xticklabels(hm.columns, rotation=30, ha="right") ax.set_yticks(range(len(hm.index))) ax.set_yticklabels(hm.index) for i in range(hm.shape[0]): for j in range(hm.shape[1]): v = hm.values[i, j] if v > 0: ax.text(j, i, str(v), ha="center", va="center", color="white" if v > hm.values.max() * 0.55 else "black", fontsize=9) plt.colorbar(im, ax=ax, label="観測所数", shrink=0.85) ax.set_title("観測所階層: 水系 × 事務所 の観測所数 (上位15水系×8事務所)") ax.set_xlabel("事務所") ax.set_ylabel("水系") plt.tight_layout() plt.savefig(ASSETS / "L04_meta_hierarchy.png", dpi=140) plt.close() # === 6. 図2: 完備性 — 観測所×時刻 の欠測ヒートマップ ========================= print("\n=== 6. 図2: 観測所×時刻 欠測ヒートマップ ===") miss_mat = tidy.isna().astype(int).T # rows=観測所, cols=時刻 n_missing_per_station = miss_mat.sum(axis=1) print(f"欠測のある観測所数: {(n_missing_per_station > 0).sum()} / {len(miss_mat)}") print(f"完全欠測 (1日通して NaN) の観測所数: {(n_missing_per_station == tidy.shape[0]).sum()}") # 欠測の多い順に並べる order = n_missing_per_station.sort_values(ascending=False).index miss_mat_sorted = miss_mat.loc[order] fig, ax = plt.subplots(figsize=(11, 6)) cmap = ListedColormap(["#ffffff", "#111111"]) ax.imshow(miss_mat_sorted.values, cmap=cmap, aspect="auto", interpolation="nearest") ax.set_xlabel("時刻 (10分刻み, 0:00 → 23:50)") ax.set_ylabel(f"観測所 (n={len(miss_mat_sorted)}, 欠測多い順)") n_xt = 12 xticks = np.linspace(0, tidy.shape[0] - 1, n_xt).astype(int) ax.set_xticks(xticks) ax.set_xticklabels([tidy.index[i].strftime("%H:%M") for i in xticks]) ax.set_yticks([]) ax.set_title(f"観測所×時刻 欠測パターン (黒=NaN) " f"全体欠測率 {miss_rate:.2f}% / 欠測のある観測所 " f"{(n_missing_per_station > 0).sum()}/{len(miss_mat)}") plt.tight_layout() plt.savefig(ASSETS / "L04_missing_heatmap.png", dpi=140) plt.close() # === 7. フラグ列の解読 ======================================================= print("\n=== 7. フラグ列 (0x...) の解読 ===") # フラグ列は header_idx 行で "フラグ" と書かれている列の中で、'10分雨量'の隣 (= cols_10min + 1) flag_cols_10min = [c + 1 for c in cols_10min] flag_data = raw.iloc[header_idx + 1:, flag_cols_10min].copy() flag_flat = flag_data.values.ravel() flag_series = pd.Series([str(v) for v in flag_flat if pd.notna(v)]) flag_freq = flag_series.value_counts() print(f"ユニークなフラグ値の数: {flag_freq.size}") print("フラグ出現頻度 (上位10):") print(flag_freq.head(10).to_string()) def parse_flag(hex_str): """0x00010000 形式を整数化。失敗したら NaN。""" try: return int(str(hex_str), 16) except (ValueError, TypeError): return np.nan flag_int = flag_series.map(parse_flag).dropna().astype(np.int64).to_numpy() # bit ごとに立っている回数を集計 (上位 32 bit を見る) N_BITS = 32 bit_counts = np.array([int(((flag_int >> b) & 1).sum()) for b in range(N_BITS)]) total_records = len(flag_int) bit_rate = bit_counts / max(total_records, 1) * 100 print(f"\n総レコード数: {total_records}, 各 bit が 1 の割合 (上位は出現多い):") for b in range(N_BITS): if bit_counts[b] > 0: print(f" bit{b:2d} (= 0x{1 << b:08X}): {bit_counts[b]:>6}件 " f"({bit_rate[b]:5.1f}%)") # === 8. 図3: フラグ bit 出現頻度 + 値別出現頻度 =============================== print("\n=== 8. 図3: フラグ bit 出現頻度 ===") fig, axes = plt.subplots(1, 2, figsize=(13, 4.6), gridspec_kw={"width_ratios": [1.1, 1]}) # (左) bit ごとの出現割合 ax = axes[0] active_bits = [b for b in range(N_BITS) if bit_counts[b] > 0] bars = ax.bar(active_bits, bit_rate[active_bits], color="#0969da", edgecolor="black", linewidth=0.6) for b, bar in zip(active_bits, bars): ax.text(bar.get_x() + bar.get_width() / 2, bar.get_height() + 0.5, f"0x{1 << b:X}", ha="center", va="bottom", fontsize=8, color="#57606a") ax.set_xticks(active_bits) ax.set_xticklabels([f"bit{b}" for b in active_bits], rotation=30, ha="right") ax.set_ylabel("出現割合 (%, 全フラグ値中)") ax.set_title(f"フラグ列の bit パターン解読 ({total_records:,}レコード)") ax.grid(alpha=0.3, axis="y") ax.set_axisbelow(True) # (右) 値別出現頻度 (上位8) ax = axes[1] top_flags = flag_freq.head(8) yp = np.arange(len(top_flags))[::-1] ax.barh(yp, top_flags.values, color="#cf222e", edgecolor="black", linewidth=0.6) ax.set_yticks(yp) ax.set_yticklabels(top_flags.index) for y, v in zip(yp, top_flags.values): ax.text(v, y, f" {v:,} ({v/total_records*100:.1f}%)", va="center", fontsize=9, color="#24292f") ax.set_xlabel("出現件数") ax.set_title(f"フラグ値 出現頻度 上位8 ({flag_freq.size} 種類中)") ax.set_axisbelow(True) ax.grid(alpha=0.3, axis="x") plt.tight_layout() plt.savefig(ASSETS / "L04_flag_decoded.png", dpi=140) plt.close() # === 9. 図4: サンプル観測所の10分時系列 重ね描き ============================== print("\n=== 9. 図4: サンプル5観測所の10分時系列重ね描き ===") # 日合計が大きい順に 5 観測所を選ぶ (= 主役) day_sum = tidy.sum(axis=0, skipna=True).sort_values(ascending=False) sample_stations = day_sum.head(5).index.tolist() print(f"日合計上位5観測所: {sample_stations}") fig, axes = plt.subplots(2, 1, figsize=(11, 6.5), sharex=True, gridspec_kw={"height_ratios": [1, 1]}) PALETTE = ["#0969da", "#cf222e", "#fb8500", "#1f883d", "#8250df"] # 上: 10分雨量 ax = axes[0] for st, c in zip(sample_stations, PALETTE): ax.plot(tidy.index, tidy[st].values, color=c, linewidth=1.4, label=f"{st} (日合計 {day_sum[st]:.1f}mm)") ax.set_ylabel("10分雨量 (mm)") ax.set_title("2024-07-01 日合計 上位 5 観測所の 10 分時系列") ax.legend(loc="upper right", fontsize=9, ncol=1) ax.grid(alpha=0.3) # 下: 累積雨量 (cumsum) ax = axes[1] for st, c in zip(sample_stations, PALETTE): ax.plot(tidy.index, tidy[st].cumsum().values, color=c, linewidth=1.6) ax.set_ylabel("累積 (10分雨量を積算, mm)") ax.set_xlabel("時刻 (2024-07-01)") ax.set_title("同 観測所の累積雨量カーブ — 「降った時間帯」が形で見える") ax.grid(alpha=0.3) plt.tight_layout() plt.savefig(ASSETS / "L04_sample_timeseries.png", dpi=140) plt.close() # === 10. 整形品質サマリ表 ==================================================== print("\n=== 10. 整形品質サマリ ===") quality_rows = [ ("生 CSV shape", f"{raw.shape[0]} × {raw.shape[1]}"), ("生 CSV 列数の内訳", f"1 (時刻) + {len(cols_10min)*4} (= {len(cols_10min)} 観測所 × 4 列)"), ("ヘッダ行の動的検出", f"row={header_idx} で '10分雨量' を発見"), ("観測所名行の検出", f"row={obs_idx}"), ("メタ行 (事務所/所管/水系/河川/番号/観測所)", f"row 0〜{obs_idx}"), ("tidy 後 shape", f"{tidy.shape[0]} 時刻 × {tidy.shape[1]} 観測所"), ("ユニーク観測所数 (重複名連番化前)", f"{meta['station'].nunique()} → 連番化後 {tidy.shape[1]}"), ("水系数 / 事務所数 / 河川数", f"{meta['system'].nunique()} / {meta['office'].nunique()} / {meta['river'].nunique()}"), ("セル全体の欠測率", f"{miss_rate:.3f}%"), ("欠測のある観測所数", f"{(n_missing_per_station > 0).sum()} / {len(miss_mat)}"), ("フラグ列のユニーク値数", f"{flag_freq.size}"), ("フラグ最多値", f"{flag_freq.index[0]} ({flag_freq.iloc[0]:,}件)"), ] quality_html = "" + "".join( f"" for k, v in quality_rows) + "
項目値
{k}{v}
" # raw / tidy / meta のプレビュー HTML preview_raw_df = raw.iloc[:8, :7].fillna("").astype(str) preview_raw_html = preview_raw_df.to_html(index=False, header=False, classes="rawtbl") preview_tidy = tidy.iloc[:6, :6].copy() preview_tidy.index = preview_tidy.index.strftime("%Y-%m-%d %H:%M") preview_tidy_html = preview_tidy.to_html() meta_preview_html = meta.head(8).to_html(index=False) scan_html = scan_df.to_html(index=False) # === 11. HTML レンダリング =================================================== print("\n=== 11. HTML レンダリング ===") sections = [ ("学習目標", """ """), ("使用データ", """ """), ("データ取得手順", data_recipe([ {"label": "雨量10分値 2024-07-01", "dataset_id": 1275, "resource_id": 94500, "file": "data/rain_2024/rain_2024-07-01.csv", "format": "CSV (5段ヘッダ, 1605列, UTF-8 BOM)", "size": "約 1.8 MB"}, ])), ("方法", """
  1. header=None で全部「データ」として読込: 多段ヘッダがある CSV は素直に read_csv せず、まず生表として確保
  2. ヘッダ位置の動的検出: 各行の "10分雨量" 出現件数を数え、最大の行をヘッダとする (位置で取らず内容で取る)
  3. 観測所メタの抽出: ヘッダより上の行を meta DataFrame に分離 (事務所/所管/水系/河川/番号/観測所)。観測所ごとの4列繰り返しの中で 10分雨量列のみを残す
  4. 値部分の数値化と tidy 化: 10分雨量の列だけ抜き、観測所名を列名に、時刻を index に
  5. 完備性チェック: 観測所×時刻 の欠測ヒートマップ (黒=NaN) で、観測停止 / センサ故障 / 欠落を可視化
  6. フラグ列の解読: 0x... を整数に変換し、bit ごとに立っている件数を集計 → 各 bit が「観測あり/欠測/異常」のどれかを推測
  7. 整形品質指標: raw shape, tidy shape, 欠測率, ユニーク観測所数 をサマリ表に
"""), ("コード解説", code(''' import re import numpy as np import pandas as pd from _common import parse_rain_csv, ensure_dataset # === (0) データ自動取得 === SRC = "data/rain_2024/rain_2024-07-01.csv" ensure_dataset(SRC, resource_id=94500, min_bytes=500_000) # === (1) header=None で全部「データ」として読込 === raw = pd.read_csv(SRC, header=None, encoding="utf-8-sig") print(raw.shape) # (150行, 1605列) ← 401観測所 × 4列 + 時刻 # === (2) ヘッダ行の動的検出 === # 年度で位置が違う (2023年式: 5行目, 2024年式: 6行目) ので # 「'10分雨量' を含む行」をスキャンして見つける。 header_idx = next(r for r in range(15) if any("10分雨量" in str(v) for v in raw.iloc[r, :])) obs_idx = next(r for r in range(header_idx) if str(raw.iloc[r, 0]).strip() == "観測所名") # === (3) 各観測所は 4 列構造。10分雨量の列だけ拾う === cols_10min = [i for i, v in enumerate(raw.iloc[header_idx, :]) if "10分雨量" in str(v)] # 401 列 stations = [str(raw.iloc[obs_idx, c]).strip() for c in cols_10min] # === (4) 観測所メタを別 DataFrame に分離 === meta = pd.DataFrame({ "office": [str(raw.iloc[0, c]) for c in cols_10min], # 事務所 "system": [str(raw.iloc[2, c]) for c in cols_10min], # 水系 "river": [str(raw.iloc[3, c]) for c in cols_10min], # 河川 "station": stations, }) # === (5) 値部分を数値化し tidy 化 === data = raw.iloc[header_idx + 1:, :].reset_index(drop=True) ts = pd.to_datetime(data.iloc[:, 0]) vals = data.iloc[:, cols_10min].apply(pd.to_numeric, errors="coerce") tidy = pd.DataFrame(vals.values, index=ts, columns=stations) # 重複観測所名は連番化 (`西部建設` が複数あるなど) # === (6) 欠測ヒートマップ用に観測所×時刻のマスクを作る === miss_mat = tidy.isna().astype(int).T miss_rate = tidy.isna().mean().mean() * 100 print(f"欠測率: {miss_rate:.2f}%") # === (7) フラグ列 (0x...) の bit 解読 === flag_cols = [c + 1 for c in cols_10min] # 10分雨量 の隣がフラグ flag_data = raw.iloc[header_idx + 1:, flag_cols] flag_int = pd.Series([int(str(v), 16) for v in flag_data.values.ravel() if pd.notna(v)], dtype=np.int64) # bit ごとに 1 が立っている件数をカウント bit_counts = [int(((flag_int >> b) & 1).sum()) for b in range(24)] # 一番よく立つ bit が「観測あり」を表す可能性が高い ''')), ("結果", f"""

① 生 CSV のスナップショット (先頭8行 × 先頭7列)

{preview_raw_html}

5 段ヘッダ + 観測時刻ヘッダ。各観測所は 4 列を占有 (10分雨量/フラグ/累計/フラグ)。空セル (NaN) が観測所間の繰り返しに見える。

② ヘッダ位置の動的検出スキャン

{scan_html}

"10分雨量" 出現件数が最大の行 = ヘッダ行。 "観測所名" で始まる行 = 観測所名行。位置で取らず内容で取るので、 2023 年式 (5 段ヘッダ) でも 2024 年式 (6 段ヘッダ) でも壊れない。

③ tidy DataFrame (観測所×時刻, 先頭6行 × 先頭6観測所)

{preview_tidy_html}

shape: {tidy.shape[0]} 時刻 × {tidy.shape[1]} 観測所。 1 行 = 1 時刻、1 列 = 1 観測所、値 = 10 分雨量 (mm)。

④ 観測所メタ DataFrame (先頭8件)

{meta_preview_html}

事務所 → 水系 → 河川 → 観測所 という 4 段の階層構造。観測所メタを別 DataFrame に切り出すと、 水系単位の集約 (L080) や河川単位の地図 (L02) と連携できる。

{figure("assets/L04_meta_hierarchy.png", "観測所階層: 水系 × 事務所 の観測所数ヒートマップ。同じ水系内に複数の事務所が観測網を持つ重層的な体制が見える")} {figure("assets/L04_missing_heatmap.png", "観測所×時刻 の欠測パターン (黒=NaN)。1 日通して欠測している観測所と、部分的に欠測する観測所がある。横筋が出ているのは「観測停止 (センサ故障)」のサイン")} {figure("assets/L04_flag_decoded.png", "フラグ列の解読。(左) 各 bit が立っている割合 — 大半が 0x00000000 で全 bit=0 (=正常観測)。立つ bit はごく一部のレコードのみで、欠測/異常状態に対応すると推測できる。(右) フラグ値別の出現頻度 上位8")} {figure("assets/L04_sample_timeseries.png", "日合計上位5観測所の10分雨量 (上) と累積雨量 (下)。tidy DataFrame からは 1 行で tidy.cumsum() と書けて即座に可視化できる")}

⑤ 整形品質サマリ

{quality_html} """), ("考察", """ """), ("発展課題", """
  1. 14日分を結合した縦長 tidy: data/rain_2024/*.csv 全14ファイルを parse_rain_csv で読み込み、pd.concat で時系列を 14 日分に伸ばす。L080 と同じ前処理スタイル。
  2. 累計雨量との突合: tidy.cumsum() と CSV 内の累計雨量列を観測所ごとに比べて、仕様書なしで「累計雨量 = 当日 0:00 からの積算」を検証する。
  3. フラグ列の意味推定 → ラベル化: bit ごとの頻度から「観測ビット (常時1)」「欠測ビット (値が NaN の時のみ1)」を統計的に推定し、flag_meaning 列を tidy に追加した拡張版を作る。
  4. 観測所メタを使った空間集約: meta DataFrame を観測所名で結合して、水系別・河川別の日合計を算出。L080 のクロス相関分析の入力になる。
  5. ヘッダ仕様の自動診断レポート: 任意の DoBoX 系現場 CSV を投げると「ヘッダ何段か」「メタ列の意味推定」「データ開始行」を出力する 汎用パーサジェネレータを実装する。
"""), ] html = render_lesson( num=4, title="雨量CSV — データ前処理の解剖 (tidy 化)", tags=["基礎", "前処理", "現場データ", "v2-rewrite"], time="100分", level="リテラシ基礎〜応用基礎", data_label='' '観測情報_雨量日集計 (DoBoX #1275) 2024-07-01分 ' '(resource_id=94500)', sections=sections, script_filename="lessons/L04_tidy_rainfall.py", ) out = LESSONS / "L04_tidy_rainfall.html" out.write_text(html, encoding="utf-8") print(f"\nsaved: {out}") print(f" PNGs: 4 generated in {ASSETS}")