PCA で雨域構造を読み解く — 時間モード vs 空間モード

学習目標

• 多変量データを行列化する 視点を持つ: 観測所×時間 と 観測所×日 のどちらの軸を取るかで何が見えるかが変わる
• PCA を生雨量フィールドに適用し、寄与率・loading・スコアを 独立した3つの可視化 で読み解く
• 主成分が 「時間モード」(全県共通の降雨タイミング) と 「空間モード」(地域差) のどちらに寄っているかを、η²(分散比)で定量化する
• 再構成誤差曲線から「上位 k 個で何 % のシグナルを再現できるか」を測り、有効次元を判断できる
• NMF (非負行列因子分解) と比較し、符号制約の違いが解釈にどう影響するかを理解する

使用データ

• 1) 雨量 10分値 (14日分) — DoBoX #1275 (観測情報_雨量日集計)。2024-06-29 〜 07-12, 約 280 観測所 × 144 時刻/日 × 14日 = 約 2016 列。
• 2) 雨量 年集計 — DoBoX #1276。本レッスンでは観測所→事務所/水系の名前マッピングのみ利用。

14日窓には 2024年7月豪雨の前駆〜本格降雨 (7/1, 7/10〜7/12) が含まれており、PCA で時間モードが立ちやすい好題材。

データ取得手順

一括取得（全レッスン共通, 推奨）:

cd "2026 DoBoX 教材"
py -X utf8 data\fetch_all.py

fetch_all.py はカタログ・追加データを data/ と data/extras/ に再現可能ダウンロード。DoBoX のオープンデータは申請不要、商用・非商用とも利用可。本レッスンの .py スクリプトは、データが無ければ自動取得してから処理を始めるよう実装されています（ensure_dataset() ヘルパ）。

スクリプト（全体ソースコード）

⬇ L08_pca_rain.py

cd "2026 DoBoX 教材"
python -X utf8 lessons/L08_pca_rain.py

方法 — 行列化の選択

14日分の10分雨量を PCA で解析するために、行と列をどう割り当てるかがまず最大の意思決定になる。

本教材では (A) 観測所 × 時間 (2016列) を採用する。理由:

1. 降雨イベントの形 (10分単位の立ち上がり/減衰) を保存できる → 第1主成分が「7月豪雨イベントそのもの」として現れることを確認できる
2. n=280 と p=2016 は p>n だが、PCA は SVD で安定計算可能 (rank ≤ min(n,p)−1 = 279)
3. 各観測所を row-center + row-scale で標準化 → 「降雨の総量」ではなく「降雨パターンの形」を比較する

1. 14ファイル統合: parse_rain_csv() で 14 日分を tidy 化し転置
2. 欠損 5% 超 / 全期間ゼロの観測所を除外
3. 各行 (観測所) を z-score 化: パターン比較のため
4. PCA(n_components=10): scree plot, loadings, scores, 再構成誤差 を計算
5. 事務所名を空間プロキシとして scatter を色分け、η²(分散比) で空間構造の強さを定量化
6. NMF(n_components=3): 非負入力で同じデータを分解し、PCA との解釈差を比較

コード解説

↑ L08_pca_rain.py 行 416–536

 1
 2
 3
 4
 5
 6
 7
 8
 9
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475

import re, numpy as np, pandas as pd
from pathlib import Path
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA, NMF
from _common import parse_rain_csv, ensure_dataset

# === (0) 14日分データを取得 ===
RAIN_RID = {"2024-06-29": 94490, ..., "2024-07-12": 94556}
for date, rid in RAIN_RID.items():
    ensure_dataset(f"data/rain_2024/rain_{date}.csv", resource_id=rid)

# === (1) 観測所×時間 行列を作る ===
files = sorted(Path("data/rain_2024").glob("rain_2024-*.csv"))
rain  = pd.concat([parse_rain_csv(f) for f in files]).sort_index()  # (時刻, 観測所)
M     = rain.T                                                       # (観測所, 時刻)
M     = M.loc[M.notna().mean(axis=1) >= 0.95].fillna(0.0)
M     = M.loc[M.sum(axis=1) > 0.5]                                   # 全期間ゼロ除外
print(M.shape)   # 例: (271, 2016)

# === (2) 行ごとに z-score 標準化 (観測所間の "降雨スケール" を揃える) ===
Xc      = M.values - M.values.mean(axis=1, keepdims=True)
row_std = Xc.std(axis=1, keepdims=True);  row_std[row_std==0] = 1.0
X       = Xc / row_std

# === (3) PCA(k=10) ===
pca = PCA(n_components=10, random_state=42)
Z   = pca.fit_transform(X)            # (n_station, 10)  観測所スコア
EVR = pca.explained_variance_ratio_   # 寄与率
print("PC1 =", EVR[0], "  cum top3 =", EVR[:3].sum())

# === (4) loadings は時間方向のベクトル → そのまま時系列としてプロット ===
times = pd.to_datetime(M.columns)
for i in range(3):
    plt.plot(times, pca.components_[i], label=f"PC{i+1}")
# loadings の "ピーク時刻" を見れば、その PC が表す時間モードが分かる
# 例: PC1 のピークが 7/2 早朝 → PC1 = 「2024-07-02 豪雨イベント」

# === (5) 観測所スコア散布: 事務所別カラーで空間構造を見る ===
office_for_station = lookup_office_from_rainfall_annual(...)
plt.scatter(Z[:,0], Z[:,1], c=office_for_station)
# 事務所が PC スコアでどれだけ説明されるかを η² で測る
def eta2(scores, groups):
    grand = scores.mean()
    sst = ((scores - grand)**2).sum()
    ssb = sum(m.sum()*(scores[m].mean()-grand)**2
              for g in np.unique(groups) for m in [groups==g])
    return ssb / sst

# === (6) 再構成誤差 vs k ===
rmse = [np.sqrt(((X - Z[:,:k] @ pca.components_[:k])**2).mean())
        for k in range(1, 11)]
# k=3 で何 % 再現できるかを RMSE で見る

# === (7) NMF 比較: 非負入力 → 解釈しやすい "降雨イベント基底" になる ===
X_pos_n = M.values / M.values.sum(axis=1, keepdims=True)   # 各観測所で正規化
nmf = NMF(n_components=3, init="nndsvda", max_iter=500, random_state=42)
W = nmf.fit_transform(X_pos_n)        # 観測所 × 3 (非負スコア)
H = nmf.components_                   # 3 × 時間 (非負基底=降雨イベント形状)
# H[i] が「i 番目の典型的な降雨イベント時系列」、W[s,i] が「観測所 s が i に
# どれだけ似ているか」 — PCA の符号反転問題が無く解釈が一意

結果

scree plot — PC1=13.5%, 上位3で 28.9%, 上位5で 37.3%, 上位10で 50.4%。1モードに集中する豪雨期間ながら、空間多様性ゆえ累積80%には10成分前後を要する。

PC1〜PC3 の loadings を時間軸でプロット。点線は日付境界。loading が強くピークする時刻 = その PC が表す時間モード。

PC1×PC2 観測所スコアの散布図。色 = 事務所 (空間プロキシ)。同じ事務所が固まって見えれば PCA は『空間構造』も捉えている。

(左) RMSE/MAE と k の関係。(右) 最大降雨観測所の波形を k=1,3,10 で再構成。k=3 でピークの位置を、k=10 で細部の波形を概ね再現。

(上) NMF が抽出した3つの非負基底 = 異なるタイミングの降雨イベント。(下) 各観測所をその主導 NMF 基底でラベル付けし PCA 散布上に重ね描き。

PC ごとのサマリ (寄与率 + 空間構造の強さ η²)

再構成誤差 vs k

考察

• 主成分は「時間モード」と「空間モード」を 同時に 捉える: PC1 の loadings は 2024-07-10 20:00 付近に鋭くピークし (=「7月10日夜の豪雨イベント」)、 さらに観測所スコアの事務所別 η² は PC1=0.55, PC2=0.47 と非常に大きい。 つまり「どの観測所が 7/10 イベントで最も濡れたか」が事務所単位 (=地理的にひとかたまりの地域) で揃っている。 PCA は "いつ降ったか" と "どこが濡れたか" を同時に説明する単一の軸として PC1 を選んだ。 これは降雨が 時空間で分離不可能 (non-separable) なフィールドであることの数学的反映である。
• PC1 と PC2 の役割は「異なる豪雨イベント」: PC1 ピークは 2024-07-10、PC2 ピークは 2024-07-01。 loading パターンを見ると、PC2 が正の観測所は「7/1 派 (前駆豪雨)」、負は「7/10 派 (本格豪雨)」に分かれる傾向。 つまり PCA は 「どちらのイベントで強く濡れたか」を観測所ごとに対比させる軸を 2 番目に立てている。 PC2 の事務所 η²=0.47 → これも空間構造を強く反映 (前線の通過位置・地形効果による地域差)。
• PC3 以降は時空間が「混ざる」: PC3 (η²=0.07) は空間構造がほぼ消え、ピークも 7/10 20:50 (PC1 と近接)。 PC3 は 同イベント内の細部 (10分単位の波形差)を表す可能性が高い。「主成分が立つほど解釈が難しくなる」典型例。
• 累積寄与率は緩やか — n << p の効果: 上位3で 28.9%、上位10で 50.4%。 RMSE は k=1 で 0.95 → k=10 で 0.81 と緩やかに減る。 これは p=2016 の自由度に対し n=392 観測所の振る舞いが多様なため。 気象フィールドでは 「上位数モードで主要イベントは説明できるが、ローカルな細部は高次成分にしか現れない」のが普通。
• NMF はイベント分解としてより直接的: NMF の各基底はそのまま「7/1 豪雨」「7/10 豪雨」などの 非負の時間集中パターンとして現れ、 観測所ごとの重み W[s,i] は「観測所 s が i 番目イベントでどれだけ濡れたか」と直読できる。 PCA の符号反転 (PC2 で正/負に分かれる) は数学的に必然だが解釈を一段複雑にする。「データが本来非負なら NMF を試す」は応用基礎の鉄則。
• 行列化の選択が結果を決める: もし観測所×日 (14列) で PCA したら、PC1 は単に「総雨量」を捉えるだけで時間構造は潰れる (旧 v1 教材の限界)。 今回 2016 列を保ったことで「PC1 = 7/10 イベント × 西部・廿日市の地域強雨」という 時空間が結びついた強い解釈ができた。 「集約のタイミングが情報構造を決める」はデータサイエンスの普遍則。

発展課題

1. 観測所緯度経度を結合 (DoBoX 別データセットや観測所マスタを別途取得) して、PC2/PC3 のスコアを地図にカラーマップする。 本教材では事務所名を空間プロキシとしたが、緯度経度があれば 北東-南西の降雨グラデーション が可視化できる。
2. 列を中心化のみ・行は中心化のみ・両方の3条件で PCA を再実行し、結果がどう変わるかを比較する (前処理の影響学習)。
3. EOF (Empirical Orthogonal Function) 解析: 気象学で使われる時空間PCAの呼称。本教材の手順と全く同じだが、地球流体力学の文脈での解釈を学ぶ。
4. ICA (Independent Component Analysis) を試す: PCA の「直交性」より「統計的独立性」を優先する分解。降雨イベントが独立な 'ソース' から来ていると仮定するなら ICA の方が物理的に合う。
5. テンソル分解 (Tucker / CP): (観測所, 日, 時刻) の3階テンソルとして直接分解すれば、日内パターンと日間パターンを同時に分離できる。
6. PCA で外れ値検出: 再構成誤差が他観測所より大きい観測所は「特異な雨域」「センサ異常」「地形効果」のいずれか。本データで誰が外れ値かを調べる。