k-means による雨量観測所のクラスタリング — 時間プロファイル × 地理分布

import pandas as pd, numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, adjusted_rand_score
from _common import parse_rain_csv, ensure_dataset

# === (0) データ自動取得 ===
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,
                   min_bytes=500_000)
ensure_dataset("data/extras/stations_master.csv",
               dataset_id=1274, resource_id=39775)

# === (1) 14日分の10分値を日合計に ===
rain_dfs   = [parse_rain_csv(f) for f in sorted(Path("data/rain_2024").glob("rain_2024-*.csv"))]
rain_full  = pd.concat(rain_dfs).sort_index()
rain_daily = rain_full.resample("D").sum(min_count=1)

# === (2) 7次元 特徴量ベクトル ===
feats = pd.DataFrame(index=rain_daily.columns)
feats["sum_mm"]    = rain_daily.sum(axis=0)
feats["max_day"]   = rain_daily.max(axis=0)
feats["std_day"]   = rain_daily.std(axis=0).fillna(0)
feats["peak_hour"] = rain_full.idxmax(axis=0).dt.hour.astype(float)
feats["wet_days"]  = (rain_daily >= 1.0).sum(axis=0)

# 観測所マスタ #1274 から lat/lon を結合
sm_raw = pd.read_csv("data/extras/stations_master.csv",
                     header=None, encoding="utf-8-sig")
sm = sm_raw.iloc[3:].copy()
sm.columns = sm_raw.iloc[2, :].tolist()
sm = sm[sm["データ種別"].astype(str)=="1"]   # 雨量
ll = sm.set_index("観測所名")[["緯度","経度"]].astype(float)
ll.columns = ["lat","lon"]
feats = feats.join(ll, how="inner")

# === (3) 3つの特徴量ビューを用意 ===
RAIN_F = ["sum_mm","max_day","std_day","peak_hour","wet_days"]
GEO_F  = ["lat","lon"]
X_rain = StandardScaler().fit_transform(feats[RAIN_F].values)
X_geo  = StandardScaler().fit_transform(feats[GEO_F].values)
# 統合は雨量を1.0、地理を0.7倍 → 雨量側を主役に保ちつつ地理も効かせる
X_both = np.hstack([X_rain,
                    StandardScaler().fit_transform(feats[GEO_F].values) * 0.7])

# === (4) エルボー + シルエットで k を選ぶ ===
KS = range(2, 11)
def scan(X):
    out = []
    for k in KS:
        km = KMeans(n_clusters=k, n_init=10, random_state=42).fit(X)
        out.append((k, km.inertia_, silhouette_score(X, km.labels_)))
    return pd.DataFrame(out, columns=["k","inertia","silhouette"])
sc_rain = scan(X_rain)
k_opt   = int(sc_rain.loc[sc_rain["silhouette"].idxmax(), "k"])

# === (5) k-means 本番 (雨量特徴ベース) ===
km = KMeans(n_clusters=k_opt, n_init=20, random_state=42).fit(X_rain)
feats["cluster"] = km.labels_

# === (6) クラスタ平均の時間プロファイル ===
for c in range(k_opt):
    members = feats[feats.cluster==c].index
    rain_daily[members].mean(axis=1).plot(label=f"c{c}", linewidth=2)

# === (7) 3パターンを ARI で比較 ===
lab_rain = KMeans(k_opt, random_state=42, n_init=20).fit_predict(X_rain)
lab_geo  = KMeans(k_opt, random_state=42, n_init=20).fit_predict(X_geo)
lab_both = KMeans(k_opt, random_state=42, n_init=20).fit_predict(X_both)
print("ARI rain↔geo  =", adjusted_rand_score(lab_rain, lab_geo))
print("ARI rain↔both =", adjusted_rand_score(lab_rain, lab_both))