2024-07-01 豪雨日 — 雨域セルの空間 - 時間進行を読む

データ取得手順

⚠️ このスクリプトは自動取得に対応していません。以下のデータセットを DoBoX から手動でダウンロードし、data/extras/ 以下に保存してください。

実行コマンド:

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

DoBoX のオープンデータは申請不要・商用/非商用とも利用可。 data/extras/ は .gitignore 対象（約 57 GB のキャッシュ）。 スクリプト実行で自動再生成されます。

学習目標と問い

このレッスンで答えたい問い

2024-07-01 の豪雨日、雨域は県内をどう動き、どこに集中し、いつピークを迎えたか？ 10分間隔の生雨量データだけから、時間と空間の両方で起きていたことを読み解けるか。

用語の定義（このレッスン独自）

立てた仮説

1. H1（雨域は移動する）: 雨域は梅雨前線の進行に従って 数時間で県内を抜けるはず。 重心の移動軌跡（24 個の点を線で結んだ図）として可視化できる
2. H2（空間クラスタ性は時刻で変動）: 強雨時間帯ほど雨は 1〜2 か所に局所集中し （Moran's I が高い）、合間は分散して降るはず
3. H3（県全体での同時ピーク）: 観測所ごとのピーク時刻は 1〜2 の時間帯に モードを持つはず （梅雨前線が県を一斉に通過するため同時多発）
4. H4（属性結合で座標欠落を代替できる）: 雨量観測所には緯経度が無いが、 河川名をキーにカメラと結合すれば、KDE/重心/Moran 計算に十分な代理座標が得られるはず
5. H5（上位観測所の波形は短時間集中型）: 日合計上位 10 観測所は、24 時間に分散ではなく 1〜2 時間にピークが集中する波形を持つはず（局所豪雨 = 短時間に降りきる）

到達点

• 10 分値 → 1 時間集約と空間 KDE を実装し、雨域の 時間進行を可視化できる
• 重心・Moran's I・ピーク時刻分布を 「ツール」として使い分けられる （数式の細部ではなく、入出力と限界を理解する）
• 属性結合による 座標欠落の代替推定のパターンを習得する（実務頻発）
• 仮説 H1〜H5 を結果と照合して 支持／反証／部分支持を判定できる

使用データ

本レッスンは 2 つのデータセットを結合して使う。雨量データには公式座標が無いため、 カメラの座標を 河川名をキーに代用する（仮説 H4 の検証対象）。

• 1) 雨量 10 分値 2024-07-01 — DoBoX #1275 (resource_id 94500)
  144 時刻 × 約 400 観測所。各観測所の 水系名 / 河川名 / 事務所名は CSV ヘッダから取れるが、緯経度は無い。
• 2) 監視カメラ一覧 — DoBoX #1279
  351 件、緯経度・路河川名等を含む。「同じ河川にあるカメラの平均緯経度」を 雨量観測所の代理座標として用いる。

なぜ 2024-07-01 か: 2018 年 7 月の西日本豪雨と同じ梅雨末期に当たる時期。 14 日窓（DoBoX #1275 の公開単位）の中で県内日合計が最大級だった日のひとつで、 時空間進行が読みやすい。

ダウンロード（再現用データ・中間データ・図）

本レッスンの全成果物に直リンクを置いた。途中ステップから再現したい学習者向け。

1. 生データ（DoBoX 由来）

2. プログラムで生成される中間データ（本レッスンの実行成果物）

3. 図 PNG

• L06_kde_smallmult.png — 図 1 時間別空間 KDE small multiples (24 枚)
• L06_com_trajectory.png — 図 2 雨域重心の軌跡と移動速度
• L06_moran_hourly.png — 図 3 時間別 Moran's I
• L06_peak_hist.png — 図 4 ピーク時刻分布と CDF
• L06_top10_timeseries.png — 図 5 Top10 観測所の 10 分値時系列

4. 再現スクリプト

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

スクリプト本体: lessons/L06_july1_heatmap.py。 データが無ければ ensure_dataset() が DoBoX から自動取得。

分析1: 雨域の空間 KDE small multiples (24時間)

狙い

「1 日の中で雨域がどこからどこへ移動したか」を、24 枚の地図を時系列で並べて一望する。 仮説 H1（雨域は移動する）を視覚的に検証するための分析。

このツールで何ができて、何ができないか（要件 J）

1 日のデータがどう変換されるか（Before/After 表、要件 K）

実装（狙いと要点）

狙い: 加重ガウシアンを各観測所に置いてグリッドに足し算するだけのシンプル実装。 要点: w > 0 の観測所だけループに入れて高速化。バンド幅 h = 0.07。

↑ L06_july1_heatmap.py 行 700–758

 1
 2
 3
 4
 5
 6
 7
 8
 9
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726

def kde_grid(weights, bandwidth=0.07, grid=80):
    """加重 KDE を等間隔グリッド上で評価。
    weights: 各観測所のその時刻の雨量 (mm)
    bandwidth: ガウシアンの「幅」(度)。0.07° ≒ 7 km
    """
    xs = np.linspace(LON_MIN, LON_MAX, grid)
    ys = np.linspace(LAT_MIN, LAT_MAX, grid)
    XX, YY = np.meshgrid(xs, ys)
    Z = np.zeros_like(XX)
    h2 = bandwidth ** 2
    for la, lo, w in zip(lats, lons, weights):
        if w <= 0:
            continue
        # 各観測所を中心とした「滑らかな小山」を足し合わせる
        Z += w * np.exp(-((XX - lo) ** 2 + (YY - la) ** 2) / (2 * h2))
    return XX, YY, Z

# 24 コマを 6 × 4 で並べて時間進化を一望 (small multiples)
fig, axes = plt.subplots(6, 4, figsize=(13, 16), sharex=True, sharey=True)
for h in range(24):
    XX, YY, Z = kde_grid(hourly.iloc[h].values, bandwidth=0.07)
    ax = axes[h // 4, h % 4]
    if Z.max() > 0:
        ax.contourf(XX, YY, Z, levels=12, cmap="Blues")
    ax.scatter(lons, lats, s=3, c="#666", alpha=0.35)  # 観測所点
    total = float(np.nansum(hourly.iloc[h].values))
    ax.set_title(f"{h:02d}:00-{h:02d}:59  Σ={total:.0f}mm", fontsize=9)

結果（図と読み取り）

なぜこの図か: 「雨域がどう移動したか」を見るには、同じ地図枠を時刻順に 24 枚並べるのが最も直感的。 1 枚の動画にすれば滑らかに見えるが、紙面では並列表示が読み取りやすい。 バンド幅 7 km 程度に揃えると、強雨セルが「色のかたまり」として現れる。

図 1: 1 時間ごとの空間 KDE 24 枚。塗り色の濃さが雨量の強さ、灰色点が観測所

この図から読み取れること:

• 時刻 03:00 前後がピーク: 県内総雨量が最大 (Σ ≈ 4686 mm)。 塗り色のかたまりが最も濃い時間帯
• 未明〜朝に強雨域が現れ、午後にかけて弱まる: 24 枚を時系列で追うと、 雨域が県の一部から始まり時間進行に伴って縮小していく
• 雨域は 1 つではなく 2 〜 3 セル同時の時刻もある — 重心 (分析 2) は これらのセルの「平均位置」を出すため注意
• 仮説 H1（雨域は移動する）への支持: 雨域は静止せず、時刻ごとに位置と強度が変化することが視覚的に確認できる

分析2: 雨域重心 (Center of Mass) の軌跡

狙い

「24 枚の KDE をたった 1 点に要約する」。各時刻の雨量を「重さ」と見なした つり合いの中心を計算し、24 個の点を線で結べば「雨域がどっちへ動いたか」が一目で分かる。 仮説 H1 を 定量化する分析。

このツールで何ができて、何ができないか（要件 J）

1 日のデータがどう変換されるか（Before/After 表、要件 K）

実装（狙いと要点）

狙い: 雨量を重みにした座標の加重平均で各時刻 1 点を出す。 要点: 雨量ゼロ時刻は NaN にして矢印を切る。差分 → km 換算で速度化。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

KM_LON = 111.0 * np.cos(np.deg2rad(34.5))   # 経度 1 度 ≒ 91 km @広島緯度
KM_LAT = 111.0                              # 緯度 1 度 ≒ 111 km

com = []
for h in range(24):
    w = np.where(np.isnan(hourly.iloc[h].values), 0, hourly.iloc[h].values)
    if w.sum() > 0:
        # 重み付き平均 = 雨量を「重さ」とした重心
        clon = float(np.average(lons, weights=w))
        clat = float(np.average(lats, weights=w))
    else:
        clon, clat = np.nan, np.nan
    com.append({"hour": h, "lon": clon, "lat": clat, "total": float(w.sum())})

com_df = pd.DataFrame(com)
# 移動速度 (km/h)
dlon = com_df["lon"].diff().values * KM_LON
dlat = com_df["lat"].diff().values * KM_LAT
com_df["speed_kmh"] = np.sqrt(dlon ** 2 + dlat ** 2)

結果（図と読み取り）

なぜこの図か: 雨域の 進行方向と速度を見るには「軌跡（散布 + 矢印）」と 「速度の時刻別棒」を並列表示するのが最も読みやすい。色は時刻 (0→23 h) で、 矢印が連続時刻の移動を示す。

図 2: 雨域重心の軌跡 (左, 色 = 時刻 0→23h) と時刻別移動速度 (右)

この図から読み取れること:

• 1 日の総移動距離 ≒ 334 km (差分の和)
• 強雨時間帯ほど重心は安定: 雨が 1 か所に集まるので重心が動かない
• 降りはじめ・降りやみ時に重心が大きく跳ぶ: 雨量が小さい時刻は分母が小さく、 わずかな分布変化で重心が振れる（重心の弱点）
• 仮説 H1 の判定: 軌跡は西→東への単調進行ではなく、複数のセルが入れ替わる動き。 H1 は 部分支持（移動はするが単純な前線通過ではない）

結果（表と読み取り）

表 1: 24 時刻の重心位置と移動速度

この表から読み取れること: 強雨時刻 ( 03:00 前後) は速度が小さく重心が安定する一方、 弱雨時刻は速度が跳ねる。重心は 「雨量の多い時刻だけ」信頼してよい指標。

分析3: 時間別 Moran's I (空間自己相関)

狙い

「近くの観測所同士で雨量が似ているか」を 1 数値に要約し、24 時刻でどう変化するかを見る。 雨が 1〜2 か所に集中している時刻は I が高く、分散して降る時刻は I が低い。 仮説 H2 を検証する。

このツールで何ができて、何ができないか（要件 J）

1 日のデータがどう変換されるか（Before/After 表、要件 K）

実装（狙いと要点）

狙い: 観測点の k 近傍 inverse-distance 重み行列で I を計算。 要点: 距離は緯経度を km に換算してから計算。重みは行標準化。

↑ L06_july1_heatmap.py 行 881–943

 1
 2
 3
 4
 5
 6
 7
 8
 9
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909

def moran_i(values, lats_arr, lons_arr, k=8):
    """k 近傍 (inverse-distance) で Global Moran's I を計算"""
    v = np.asarray(values, dtype=float)
    mask = ~np.isnan(v)
    if mask.sum() < 10 or np.nanstd(v) == 0:
        return np.nan
    v = v[mask]
    la, lo = lats_arr[mask], lons_arr[mask]
    n = len(v)
    # 距離行列 (km)
    dla = (la[:, None] - la[None, :]) * KM_LAT
    dlo = (lo[:, None] - lo[None, :]) * KM_LON
    d = np.sqrt(dla ** 2 + dlo ** 2)
    np.fill_diagonal(d, np.inf)
    # 各点の k 近傍だけ inverse-distance 重みを残す
    W = np.zeros_like(d)
    for i in range(n):
        idx = np.argsort(d[i])[:k]
        W[i, idx] = 1.0 / np.maximum(d[i, idx], 0.1)
    # 行標準化 (各行の合計を 1 に)
    W = W / W.sum(axis=1, keepdims=True)
    z = v - v.mean()
    s2 = (z ** 2).sum() / n
    if s2 == 0:
        return np.nan
    num = (W * np.outer(z, z)).sum()
    return float(num / (s2 * n) * (n / W.sum()))

moran_per_hour = [moran_i(hourly.iloc[h].values, lats, lons, k=8) for h in range(24)]

結果（図と読み取り）

なぜこの図か: I の時刻変化と 時間総雨量を 同じ時間軸で並列に見たい。 左軸 = I の棒（赤=強・青=中・灰=弱）、右軸 = 総雨量の折れ線。 両者が同期するか／ずれるかが仮説 H2 の核心。

図 3: 時間別 Moran's I (棒) と時間総雨量 (橙線)。同期すれば「雨多い時刻 = 集中」

この図から読み取れること:

• 最大 I = 0.77 @ 03:00: その時刻に雨域が最も局所集中していた
• 日合計 I = 0.65: 1 日全体で見ても クラスタ的な配置 (>0)
• I が高い時刻 ≒ 強雨時間帯: 仮説 H2（強雨時間帯ほど集中）に 支持の傾向
• 低 I 時刻 = 雨域が拡散・併合 or 弱雨で観測ノイズ的: 雨域が消滅していく時間帯にも見られる

結果（表と読み取り）

表 2: 24 時刻の Moran's I と総雨量

この表から読み取れること: 「強い局所集中」と判定された時刻は雨量も多い傾向。 仮説 H2 への部分支持（同期はするが完全に一致するわけではない）。

分析4: 観測所別ピーク時刻ヒストグラム

狙い

「県内の観測所たちがいつピーク雨量を記録したか」を 1 つのヒストグラムに集約する。 1〜2 の時間帯に観測所が集中していれば 梅雨前線の同時通過を示唆（仮説 H3）。 広く分散していれば 独立した複数のセルが時間差で降ったことを示す。

狙いに対する手法選択（要件 H）

• 「県全体での集中時間帯」を見たい → 観測所別 1 時刻に縮約してヒストグラム
• 累積分布 (CDF) も並列表示すると「P25/P50/P75 の時刻」が読める → ピーク到達の早さ・遅さを定量化
• str.contains や K-Means のような複雑手法は不要 — idxmax 一発で OK

1 日のデータがどう変換されるか（Before/After 表、要件 K）

実装（狙いと要点）

狙い: 各観測所のピーク時刻を集計して分布を見る。 要点: 日合計が小さい観測所はピーク時刻に意味がないので除外する。

↑ L06_july1_heatmap.py 行 347–372

 1
 2
 3
 4
 5
 6
 7
 8
 9
356
357

peak_hours = []
for st in geo_cols:
    s = hourly[st]
    if s.sum() > 5:                          # 日合計 5 mm 超だけ採用
        peak_hours.append(int(s.idxmax().hour))
peak_hours = np.array(peak_hours)

# 24 ビンのヒストグラム + 累積分布 (CDF)
counts, edges, _ = ax.hist(peak_hours, bins=np.arange(-0.5, 24.5, 1))
mode_h = int(np.argmax(counts))              # 最頻時間帯
ax2.plot(np.sort(peak_hours), np.arange(1, len(peak_hours)+1)/len(peak_hours)*100)

結果（図と読み取り）

なぜこの図か: 同時多発性を見るには ヒストグラムが定石。 モード（最頻時間帯）を縦線で示し、CDF を併置すると「P25/P50/P75 はいつか」が読める。

図 4: 観測所別ピーク時刻のヒストグラム (左) と累積分布 CDF (右)

この図から読み取れること:

• モード時刻 = 03:00 (133 観測所): 県内の最も多くの観測所が同じ時間帯にピークを迎えた
• ピーク時刻の P25-P75 = 3.0h 〜 4.0h: この幅が狭いほど「短時間集中型」、広いほど「長時間持続型」
• 仮説 H3 の判定: モードに観測所が集中 → 支持。ただし裾も長く、雨域が単純な 1 セルではないことを示す
• 避難判断への含意: モード時刻に県全体で同時に避難判断が必要になるシナリオ

結果（表と読み取り）

表 3: ピーク時刻 上位 5 時間帯（観測所数の多い順）

この表から読み取れること: 上位 5 時間帯だけで全観測所の過半が説明できる。 雨域の主要なピークはこの数時間に集中していたことが定量的に確認できる。

分析5: Top10 観測所の 10 分値時系列と座標推定

狙い

「日合計上位 10 観測所の波形」を 10 分値で重ね描きし、ピーク集中度を観測所別に見る。 仮説 H5（上位観測所は短時間集中型）を検証。 同じ分析の中で、仮説 H4（属性結合で座標欠落を代替できる）の検証も行う （10 分値は座標を使わないが、「Top10 を地図上に置きたい」発展課題で必要になる）。

STEP 1: 河川名による属性結合 — 座標欠落の代替推定（仮説 H4）

表 4: 座標解決ソースの内訳（実データ）

この表から読み取れること: 河川一致 (227 観測所) が大半を占め、 未解決は 88 観測所。仮説 H4 への支持: 属性結合だけで KDE/重心/Moran 計算に十分な座標が得られた。

STEP 2: Top10 観測所の 10 分値波形

狙い: 日合計上位 10 観測所を抽出し、10 分値の波形を重ね描き。 ピークがどれくらい鋭いか、複数ピークか単峰かを見る。

1 日のデータがどう変換されるか（Before/After 表、要件 K）

実装（狙いと要点）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

# Top10 抽出と 10 分値の重ね描き
top10_names = day_total[geo_cols].sort_values(ascending=False).head(10).index.tolist()

fig, ax = plt.subplots(figsize=(12, 5))
cmap = plt.get_cmap("tab10")
for i, st in enumerate(top10_names):
    s = tidy[st]
    ax.plot(s.index, s.values, color=cmap(i), alpha=0.85, lw=1.3,
            label=f"{st} (Σ={day_total[st]:.0f}mm)")

# 集中度 = ピーク 1 時間雨量 / 日合計 (1 に近いほど短時間集中型)
for st in top10_names:
    peak_share = hourly[st].max() / day_total[st]

結果（図と読み取り）

なぜこの図か: 上位観測所の 波形の鋭さを見るには 10 分値を 重ね描きするのが最も直感的。 1 観測所 1 線で色分けし、凡例に日合計を出す。波形が短時間に立ち上がる観測所が短時間集中型。

図 5: 日合計上位 10 観測所の 10 分値時系列重ね描き

この図から読み取れること:

• ほとんどの観測所が同じ時間帯にピーク: 重ね描きでピーク位置が揃う — 仮説 H3 と整合的
• 1 観測所内では数十分のスパイク状ピークが中心 — 1 時間以下のスケールの強雨
• 仮説 H5 の判定: 表 5 の「最大時間/日合計」がおおよそ 20〜37% → 部分支持 （1 時間でその日の 2〜3 割が降る = 2〜4 時間に降りきる短時間集中型）

結果（表と読み取り）

表 5: 日合計上位 10 観測所のサマリ（集中度付き）

この表から読み取れること: 「最大時間/日合計」が 20〜37% に分布 （3 時間で日合計の 6〜9 割が降った計算）。 水系列を見ると上位は 「単独河川」「その他（沿岸部）」「小瀬川」に偏っており、 県西部沿岸〜山地のラインに集中する地形性降水のサイン。 発展課題（地形・標高との照合）への伏線になる。

仮説検証と考察

仮説と結果の照合

考察

• 雨は「動く」だけでなく「同期する」: 図 1〜3 の組み合わせから、雨域は 時間移動しつつ 複数地点で同時にピークを迎える複雑な振る舞いを示す。 単純な前線通過モデルでは説明できない部分が残り、L05/L06 の発展（多日比較・気象レーダ照合）が必要。
• 「同じ日合計」でも空間クラスタ性で被害の出方が違う: Moran's I が高い時刻 = 1〜2 か所に集中（狭い谷でフラッシュフラッドのリスク）、 低い時刻 = 広域に分散（多くの河川で同時に水位上昇のリスク）。 合計だけでなく クラスタ性も防災判断に必要。
• 属性結合の威力と限界: 仮説 H4 の支持から、河川名のような属性キーは 座標欠落データを救う有効な手段だと分かった。一方で「同じ河川でも上流/下流で雨量が違う」という 誤差は残り、点状の精度を要する分析（特定地点の警報判断など）には不十分。 fallback 階層を持つ設計は実務でも有用。
• KDE のバンド幅は「見たい現象のスケール」: h=0.07°≒7 km は強雨セル 1 つの典型サイズ。 これより細かくすれば点状（観測所配置のバイアスが見える）、太くすれば全県塗り潰し（時間進行が消える）。 「どの空間スケールで雨を見たいか」を決めることがバンド幅選びと同義。
• 方法の限界と次の一歩: (1) Global Moran's I は どこに集中するかを答えない → Local Moran's I (発展課題)。 (2) 重心は 2 セル併存に弱い → 複数モード抽出 (k-means 重心) が次の選択肢。 (3) 観測点ベース KDE は 山陰の観測薄で過小評価 → 気象レーダ照合 (発展課題)。

発展課題（結果から導かれる新たな問い）

各課題は、上の 結果と 新しい仮説に裏打ちされている。 「結果 X → 新仮説 Y → 課題 Z」の 3 段で書く。

1. 多日比較で前線の動きの再現性を見る
   • 結果 X: 図 2 の重心軌跡は 2024-07-01 に固有のパターン。前線通過のたびに同じ動きをするのか不明
   • 新仮説 Y: 14 日窓の他の豪雨日 (06-30, 07-02 など) の重心軌跡を重ねると、 梅雨前線の典型的進行方向が浮かび上がるはず（ベクトル場として可視化）
   • 課題 Z: rain_2024 ディレクトリの 14 日分について同じ重心計算を行い、 14 本の軌跡を 1 枚に重ねて方向の主成分を取る (PCA)
2. Local Moran's I でホットスポットを特定
   • 結果 X: 図 3 の Global I は 「県全体での集中度」1 個しか示さず、 具体的な集中地点が分からない
   • 新仮説 Y: 観測所ごとに Local Moran's I (LISA) を計算すれば、 HH (高×高クラスタ) / LH (Outlier) / LL (低×低クラスタ) として 地図上で どこが豪雨セルの中心かが直接読めるはず
   • 課題 Z: PySAL の esda.Moran_Local で計算 → 4 タイプを地図に色分け。 最強雨時刻 (03:00) で確認
3. KDE バンド幅の感度分析
   • 結果 X: 図 1 は h=0.07° で固定だが、強雨セルのサイズは時刻で変わる可能性がある
   • 新仮説 Y: 各時刻で 適応的バンド幅 (Silverman's rule や交差検証) を使えば、 強雨時刻ほど自動的に細かい解像度になり、KDE の見かけが変わるはず
   • 課題 Z: scipy.stats.gaussian_kde の bw_method を変えて h=0.03/0.07/0.15 の 3 通りを並列描画。差分から「現象スケール」が読めるかを検証
4. 観測所→流域水位ピーク時刻のラグ
   • 結果 X: 図 4 で観測所のピーク時刻は 03:00 にモード。 これに対し河川水位ピークは何時間遅れるのか不明
   • 新仮説 Y: 観測所ピーク → 河川水位ピークの 時間ラグは流域面積でスケールする （大きい流域ほど遅れて応答）。観測所と水位観測所をペアにして散布すると単調関係が見える
   • 課題 Z: 水位 10 分値データ (DoBoX 1276 系) を取得し、流域別にピーク時刻を取って 雨量ピーク時刻との差分を散布。流域面積 (河川シェイプ) と重ね合わせる
5. Top10 観測所の地形配置と短時間集中度の関係
   • 結果 X: 表 5 で上位観測所は特定の水系・河川に偏り、最大時間/日合計が高い (短時間集中型)
   • 新仮説 Y: 短時間集中型の観測所は 急峻な谷に位置するはず（地形性降水）。 標高や傾斜と相関するはず
   • 課題 Z: 標高 DEM (DoBoX) を取得し、Top10 観測所の代理座標で 標高・傾斜を抽出。最大時間/日合計と散布図で相関を取る