""" L03: 避難所4065件 JSON→DataFrame — 構造探索・災害種別共起・収容力・バリアフリー DoBoX #42 「避難所情報」(4,065 施設・35+ フィールド) の JSON を ``pd.json_normalize`` で平坦化し、以下を多角的に提示する v2 教材: 1. 全フィールドの dtype/null率/ユニーク数 を表化 (構造探索) 2. 5 災害種別フラグの **共起行列 (Jaccard)** ヒートマップ 3. capacity (収容人数) の対数正規フィット + Q-Q プロット 4. 市町村別カバレッジ — capacity 合計 vs 件数 散布 + 上位15バー 5. folium による災害種別カラーマーカー (代表災害ごとの色分け) 6. バリアフリー総合スコア (11 フラグ合計) の市町村別箱ひげ 実行: cd "2026 DoBoX 教材" py -X utf8 lessons/L03_shelter_analysis.py """ from pathlib import Path import json import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.colors import LinearSegmentedColormap import folium from folium.plugins import MarkerCluster from scipy import stats from _common import (ROOT, ASSETS, LESSONS, render_lesson, code, figure, data_recipe, ensure_dataset) plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False # === 0. データ自動取得 ======================================================= SHELTER_PATH = ROOT / "data" / "shelters.json" print("=== 0. データ自動取得 ===") ensure_dataset(SHELTER_PATH, dataset_id=42, label="避難所情報 #42") # === 1. JSON → DataFrame ===================================================== print("\n=== 1. JSON → DataFrame (json_normalize) ===") with open(SHELTER_PATH, encoding="utf-8") as f: raw = json.load(f) df = pd.json_normalize(raw["items"]) print(f"shape: {df.shape} (レコード={len(df)}, フィールド={df.shape[1]})") # 数値化 for c in ["latitude", "longitude", "capacity", "westernNum", "japaneseNum", "parkingNum"]: if c in df.columns: df[c] = pd.to_numeric(df[c], errors="coerce") # 災害種別 5 フラグ DISASTER = { "洪水": "floodShFlg", "土砂": "sedimentDisasterShFlg", "高潮": "stormSurgeShFlg", "地震": "earthquakeShFlg", "津波": "tsunamiShFlg", } for k, c in DISASTER.items(): df[k] = (df[c].astype(str) == "1").astype(int) # バリアフリー / 居住性 11 フラグ (>=1 を positive とする;'2'=「一部」も含む) BF_COLS = ["handicappedFlg", "ostomateFlg", "petFlg", "parkingFlg", "internetFlg", "bathFlg", "breastfeedingFlg", "powerFlg", "cookingFlg", "heatingFlg", "coolongFlg"] for c in BF_COLS: df[c + "_b"] = (pd.to_numeric(df[c], errors="coerce").fillna(0) >= 1).astype(int) df["bf_score"] = df[[c + "_b" for c in BF_COLS]].sum(axis=1) # === 2. 構造探索: 全フィールドの dtype / null率 / unique 数 =================== print("\n=== 2. 全フィールド プロファイル ===") profile_rows = [] for c in df.columns: if c in DISASTER or c.endswith("_b") or c == "bf_score": continue s = df[c] profile_rows.append({ "field": c, "dtype": str(s.dtype), "null率(%)": round(s.isna().mean() * 100, 1), "ユニーク数": s.nunique(dropna=True), "例": str(s.dropna().iloc[0])[:30] if s.notna().any() else "-", }) profile = pd.DataFrame(profile_rows) print(profile.to_string(index=False)) profile_html = profile.to_html(index=False, classes="profile") # === 3. 災害種別共起行列 (Jaccard) ============================================ print("\n=== 3. 災害種別 Jaccard 共起行列 ===") keys = list(DISASTER) J = np.zeros((len(keys), len(keys))) for i, a in enumerate(keys): for j, b in enumerate(keys): inter = ((df[a] == 1) & (df[b] == 1)).sum() union = ((df[a] == 1) | (df[b] == 1)).sum() J[i, j] = inter / union if union else 0 jaccard_df = pd.DataFrame(J, index=keys, columns=keys) print(jaccard_df.round(3).to_string()) fig, ax = plt.subplots(figsize=(6.4, 5.4)) cmap = LinearSegmentedColormap.from_list( "jaccard", ["#f0f8ff", "#0969da", "#0a3069"]) im = ax.imshow(J, cmap=cmap, vmin=0, vmax=1, aspect="equal") ax.set_xticks(range(len(keys))) ax.set_yticks(range(len(keys))) ax.set_xticklabels(keys, fontsize=11) ax.set_yticklabels(keys, fontsize=11) for i in range(len(keys)): for j in range(len(keys)): ax.text(j, i, f"{J[i, j]:.2f}", ha="center", va="center", color="white" if J[i, j] > 0.55 else "#1f2328", fontsize=10, fontweight="bold") plt.colorbar(im, ax=ax, label="Jaccard 係数 |A∩B|/|A∪B|", shrink=0.85) ax.set_title("災害種別フラグの共起 (Jaccard)\n対角=自己一致, 非対角=同時対応の重なり度") plt.tight_layout() plt.savefig(ASSETS / "L03_jaccard.png", dpi=140) plt.close() # 災害種別カウントの棒も並走 (どの災害が何件か) counts = pd.Series({k: int(df[k].sum()) for k in keys}).sort_values() fig, ax = plt.subplots(figsize=(7, 3.8)) bars = ax.barh(counts.index, counts.values, color=["#cf222e", "#fb8500", "#0969da", "#1f883d", "#8250df"]) ax.set_xlabel("対応避難所数") ax.set_title(f"災害種別カバレッジ (全 {len(df):,} 件中)") for i, v in enumerate(counts.values): ax.text(v + 30, i, f"{v:,} ({v/len(df)*100:.1f}%)", va="center", fontsize=10) ax.set_xlim(0, counts.max() * 1.18) plt.tight_layout() plt.savefig(ASSETS / "L03_disaster_counts.png", dpi=140) plt.close() # === 4. capacity の分布: ヒスト + 対数正規フィット + Q-Q ===================== print("\n=== 4. capacity 分布解析 ===") cap = df["capacity"].dropna() cap = cap[cap > 0] print(f" n={len(cap)}, mean={cap.mean():.1f}, median={cap.median():.0f}, " f"max={cap.max():.0f}, p95={cap.quantile(0.95):.0f}") # 対数正規フィット (shape=sigma, scale=exp(mu)) shape, loc, scale = stats.lognorm.fit(cap, floc=0) mu, sigma = np.log(scale), shape print(f" lognorm fit: μ_log={mu:.3f}, σ_log={sigma:.3f}, " f"中央値={np.exp(mu):.0f}人") fig, axes = plt.subplots(1, 2, figsize=(11.5, 4.4)) # (左) ヒスト + フィット曲線 ax = axes[0] log_cap = np.log10(cap) n_bins = 40 ax.hist(log_cap, bins=n_bins, color="#85e89d", edgecolor="#1a7f37", density=True, label="観測 (log10)") xs = np.linspace(log_cap.min(), log_cap.max(), 300) # log10(x) のヒストに重ねるには、x = 10**xs での lognorm pdf を変換 pdf = stats.lognorm.pdf(10 ** xs, shape, loc, scale) * (10 ** xs) * np.log(10) ax.plot(xs, pdf, color="#cf222e", linewidth=2.2, label=f"対数正規 fit (μ={mu:.2f}, σ={sigma:.2f})") ax.axvline(np.log10(cap.median()), color="#0969da", linestyle="--", label=f"中央値 {cap.median():.0f}人") ax.set_xlabel("log10(収容人数)") ax.set_ylabel("密度") ax.set_title(f"収容人数の分布 (n={len(cap):,}) と対数正規 fit") ax.legend(fontsize=9, loc="upper right") ax.grid(alpha=0.3) # (右) Q-Q プロット (理論分位 vs 観測分位 on log scale) ax = axes[1] log_cap_sorted = np.sort(np.log(cap.values)) n = len(log_cap_sorted) theor_q = stats.norm.ppf((np.arange(1, n + 1) - 0.5) / n, loc=mu, scale=sigma) ax.scatter(theor_q, log_cap_sorted, s=6, alpha=0.45, color="#0969da") lo, hi = min(theor_q.min(), log_cap_sorted.min()), max(theor_q.max(), log_cap_sorted.max()) ax.plot([lo, hi], [lo, hi], color="#cf222e", linewidth=1.5, label="y=x (完全一致)") ax.set_xlabel("理論分位 (対数正規 ln-scale)") ax.set_ylabel("観測分位 ln(capacity)") ax.set_title("Q-Q プロット: capacity ~ Lognormal?") ax.legend(fontsize=9) ax.grid(alpha=0.3) plt.tight_layout() plt.savefig(ASSETS / "L03_capacity_dist.png", dpi=140) plt.close() # Kolmogorov-Smirnov 検定 (参考) ks_stat, ks_p = stats.kstest(cap.values, "lognorm", args=(shape, loc, scale)) print(f" KS test: D={ks_stat:.4f}, p={ks_p:.2e}") # === 5. 市町村別カバレッジ =================================================== print("\n=== 5. 市町村別カバレッジ ===") muni = df.groupby("municipalityName").agg( n=("name", "size"), cap_total=("capacity", "sum"), cap_med=("capacity", "median"), flood=("洪水", "sum"), sediment=("土砂", "sum"), tsunami=("津波", "sum"), bf_med=("bf_score", "median"), ).sort_values("n", ascending=False) print(muni.head(15).to_string()) muni_summary_html = muni.head(15).round(1).to_html() # (図) 散布: 件数 vs capacity 合計 fig, axes = plt.subplots(1, 2, figsize=(13, 5.0)) ax = axes[0] ax.scatter(muni["n"], muni["cap_total"] / 1000, s=40, alpha=0.55, color="#0969da", edgecolor="white") # 上位 8 市町村を注釈 for name, row in muni.head(8).iterrows(): ax.annotate(name, (row["n"], row["cap_total"] / 1000), fontsize=9, color="#1f2328", xytext=(4, 3), textcoords="offset points") ax.set_xlabel("避難所件数") ax.set_ylabel("収容人数合計 (千人)") ax.set_title(f"市町村別: 件数 × 収容力合計 ({len(muni)} 市町村)") ax.grid(alpha=0.3) # (図) 上位 15 市町村バー ax = axes[1] top15 = muni.head(15) ypos = np.arange(len(top15)) ax.barh(ypos, top15["n"], color="#0969da", label="件数") ax.barh(ypos, top15["cap_total"] / 1000, left=top15["n"], color="#fb8500", alpha=0.7, label="収容人数 (千人)") ax.set_yticks(ypos) ax.set_yticklabels(top15.index, fontsize=9) ax.invert_yaxis() ax.set_xlabel("件数 / 収容人数 (千人)") ax.set_title("上位 15 市町村") ax.legend(loc="lower right") ax.grid(alpha=0.3, axis="x") plt.tight_layout() plt.savefig(ASSETS / "L03_muni_coverage.png", dpi=140) plt.close() # === 6. バリアフリー総合スコア — 市町村別 箱ひげ ============================== print("\n=== 6. バリアフリー総合スコア (11 flags) ===") print(f" 全体: mean={df['bf_score'].mean():.2f}, " f"median={df['bf_score'].median():.0f}, max={df['bf_score'].max()}") # 件数上位 12 市町村のみ箱ひげ (見やすさのため) target_munis = muni.head(12).index.tolist() groups = [df.loc[df["municipalityName"] == m, "bf_score"].values for m in target_munis] fig, ax = plt.subplots(figsize=(11, 5.0)) bp = ax.boxplot(groups, labels=target_munis, patch_artist=True, showmeans=True, meanprops={"marker": "D", "markerfacecolor": "#cf222e", "markeredgecolor": "white"}) palette = plt.cm.tab20(np.linspace(0, 1, len(target_munis))) for patch, c in zip(bp["boxes"], palette): patch.set_facecolor(c) patch.set_alpha(0.7) ax.set_ylabel("バリアフリー総合スコア (0〜11)") ax.set_title("市町村別 バリアフリー総合スコア分布 (赤◇=平均, 件数上位12市町村)") ax.grid(axis="y", alpha=0.3) plt.setp(ax.get_xticklabels(), rotation=30, ha="right", fontsize=9) plt.tight_layout() plt.savefig(ASSETS / "L03_bf_boxplot.png", dpi=140) plt.close() # === 7. folium 地図 — 災害種別カラーマーカー ================================== print("\n=== 7. folium 地図 (災害種別カラー) ===") df_map = df.dropna(subset=["latitude", "longitude"]).copy() df_map = df_map[(df_map["latitude"].between(33, 36)) & (df_map["longitude"].between(131, 134))] print(f" 地理表示対象: {len(df_map):,} 件") # 「最も特徴的な災害」を 1 つ選んで色分け (優先順: 津波→土砂→洪水→高潮→地震) PRIORITY = [("津波", "#cf222e"), ("土砂", "#fb8500"), ("洪水", "#0969da"), ("高潮", "#8250df"), ("地震", "#1f883d")] def primary_color(row): for k, c in PRIORITY: if row[k] == 1: return k, c return "対応なし", "#999999" df_map[["primary", "color"]] = df_map.apply( lambda r: pd.Series(primary_color(r)), axis=1) m = folium.Map([df_map["latitude"].mean(), df_map["longitude"].mean()], zoom_start=9, tiles="cartodbpositron") # 凡例 legend_html = """
主災害カテゴリ (優先順)
● 津波対応
● 土砂対応
● 洪水対応
● 高潮対応
● 地震のみ
● 5種別いずれも非対応
""" m.get_root().html.add_child(folium.Element(legend_html)) mc = MarkerCluster(name="避難所", disable_clustering_at_zoom=13).add_to(m) for _, r in df_map.iterrows(): badges = " ".join( [f"{k}" for k in keys if r[k]]) pop = folium.Popup( f"{r['name']}
" f"{r['address01']}{r['address02']}
" f"収容: {int(r['capacity']) if pd.notna(r['capacity']) else '不明'}人
" f"BFスコア: {int(r['bf_score'])}/11
" f"主区分: {r['primary']}
{badges}", max_width=320) folium.CircleMarker([r["latitude"], r["longitude"]], radius=3.5, color=r["color"], weight=1, fill=True, fill_color=r["color"], fill_opacity=0.78, popup=pop).add_to(mc) map_path = ASSETS / "L03_map.html" m.save(str(map_path)) print(f" saved: {map_path}") # === 8. HTML レンダリング =================================================== print("\n=== 8. HTML レンダリング ===") sections = [ ("学習目標", """ """), ("使用データ", """ """), ("データ取得手順", data_recipe([ {"label": "避難所情報 (4,065 施設)", "dataset_id": 42, "file": "data/shelters.json", "format": "JSON (items 配列, UTF-8)", "size": "約 4 MB"}, ])), ("方法", """
  1. JSON 読込 → 平坦化: json.load → pd.json_normalize(raw["items"]) で 4,065 × 36 の DataFrame を得る。latitude / longitude / capacity は文字列で来るため pd.to_numeric(errors="coerce")
  2. 構造プロファイル: 各列の dtype・null率・ユニーク数を 1 表にまとめる ('field profiling' = データ理解の最初の一歩)
  3. 5 災害種別フラグを {0,1} に整数化し、件数・カバレッジ率を集計
  4. Jaccard 共起行列: |A∩B| / |A∪B| を 5×5 で計算 → ヒートマップ。同時対応のしやすさが定量化される
  5. capacity 分布: scipy.stats.lognorm.fit で μ_log, σ_log を推定 → ヒストに pdf を重ね、Q-Q プロットで末尾の乖離を可視化
  6. 市町村集約: groupby("municipalityName") で件数・capacity合計・BFスコア中央値 → 散布 + 上位 15 バーで二側面提示
  7. バリアフリー総合スコア: 11 フラグの >=1 を集計 (0〜11) → 件数上位 12 市町村で箱ひげ比較
  8. folium 地図: 主災害カテゴリ (優先順 津波→土砂→洪水→高潮→地震) で 5 色に塗り分け、MarkerCluster で 4,065 点を破綻なく描画
"""), ("コード解説", code(''' import json import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy import stats import folium from folium.plugins import MarkerCluster from _common import ensure_dataset # === (0) データ自動取得 === ensure_dataset("data/shelters.json", dataset_id=42, label="避難所情報 #42") # === (1) JSON → DataFrame: json_normalize で平坦化 === with open("data/shelters.json", encoding="utf-8") as f: raw = json.load(f) df = pd.json_normalize(raw["items"]) # (4065, 36) for c in ["latitude", "longitude", "capacity"]: df[c] = pd.to_numeric(df[c], errors="coerce") # === (2) 災害種別 5 フラグを整数化 === DISASTER = {"洪水":"floodShFlg", "土砂":"sedimentDisasterShFlg", "高潮":"stormSurgeShFlg", "地震":"earthquakeShFlg", "津波":"tsunamiShFlg"} for k, c in DISASTER.items(): df[k] = (df[c].astype(str) == "1").astype(int) # === (3) フィールド プロファイリング === profile = pd.DataFrame([{ "field": c, "dtype": str(df[c].dtype), "null率(%)": round(df[c].isna().mean()*100, 1), "ユニーク数": df[c].nunique(), } for c in df.columns]) # === (4) Jaccard 共起行列: |A∩B| / |A∪B| === keys = list(DISASTER) J = np.zeros((5, 5)) for i, a in enumerate(keys): for j, b in enumerate(keys): inter = ((df[a]==1) & (df[b]==1)).sum() union = ((df[a]==1) | (df[b]==1)).sum() J[i, j] = inter / union if union else 0 # imshow で 5x5 ヒートマップ → ファイル保存 # === (5) capacity の対数正規フィットと Q-Q プロット === cap = df["capacity"].dropna(); cap = cap[cap > 0] shape, loc, scale = stats.lognorm.fit(cap, floc=0) mu, sigma = np.log(scale), shape # log-平均と log-SD # 観測 vs 理論の Q-Q log_sorted = np.sort(np.log(cap.values)) n = len(log_sorted) theor = stats.norm.ppf((np.arange(1, n+1) - 0.5) / n, loc=mu, scale=sigma) plt.scatter(theor, log_sorted, s=6) # 直線に乗れば lognormal が良 fit # === (6) 市町村集約: 件数 × capacity合計 × BFスコア === muni = df.groupby("municipalityName").agg( n=("name", "size"), cap_total=("capacity", "sum"), bf_med=("bf_score", "median"), ).sort_values("n", ascending=False) # === (7) folium: 主災害カテゴリで 5 色に塗り分け === PRIORITY = [("津波","#cf222e"), ("土砂","#fb8500"), ("洪水","#0969da"), ("高潮","#8250df"), ("地震","#1f883d")] def primary(row): for k, c in PRIORITY: if row[k] == 1: return c return "#999" df["color"] = df.apply(primary, axis=1) m = folium.Map([34.4, 132.5], zoom_start=9, tiles="cartodbpositron") mc = MarkerCluster(disable_clustering_at_zoom=13).add_to(m) for _, r in df.dropna(subset=["latitude","longitude"]).iterrows(): folium.CircleMarker([r.latitude, r.longitude], radius=3.5, color=r.color, fill=True, fill_color=r.color, fill_opacity=0.78).add_to(mc) m.save("lessons/assets/L03_map.html") ''')), ("結果", "

(1) 全フィールド プロファイル

" "

36 列のうち shelterId / shelterStartTimestamp / " "shelterEndTimestamp は 100% null (開設状況の予約フィールド)。" "capacity の null は 535 件 (13%) — 「収容力不明」の施設が一定数ある。

" + profile_html + figure("assets/L03_disaster_counts.png", f"5 災害種別カバレッジ。最多は洪水 {int(df['洪水'].sum()):,} 件 ({df['洪水'].mean()*100:.1f}%)、" f"最少は津波 {int(df['津波'].sum()):,} 件 ({df['津波'].mean()*100:.1f}%)。" "津波が極端に少ないのは沿岸市町に偏在しているため") + figure("assets/L03_jaccard.png", "災害種別間 Jaccard 共起行列。洪水×土砂が約 0.59 と最も重なり大、" "津波は他カテゴリと小さい (海岸線の地理的制約)") + figure("assets/L03_capacity_dist.png", f"(左) 収容人数の log10 ヒストに対数正規 fit (μ={mu:.2f}, σ={sigma:.2f}, " f"中央値 {np.exp(mu):.0f} 人) を重ね描き／" "(右) Q-Q プロット — 中央域は直線にほぼ乗るが裾 (大規模施設) は乖離。" "「全体は対数正規だが極大はべき的」という現象を可視化") + figure("assets/L03_muni_coverage.png", "(左) 市町村別: 件数 × 収容人数合計の散布。広島市・福山市・東広島市が圧倒的／" "(右) 上位 15 市町村のスタックバー (青=件数, オレンジ=収容人数(千人))") + figure("assets/L03_bf_boxplot.png", "件数上位 12 市町村のバリアフリー総合スコア (0〜11) 箱ひげ。" "中央値は概ね 1〜2 と低く、市町村間で分布の散らばり方が違う点に注目") + '

(6) 上位 15 市町村 集計表

' + muni_summary_html + '

(7) 避難所マップ — 主災害カテゴリ別カラー

' '

画面拡大で個別マーカーが現れる (zoom 13 以上で MarkerCluster 解除)。' 'クリックで施設名・住所・収容人数・BFスコア・対応災害種別を表示。

' ''), ("考察", """ """), ("発展課題", """
  1. 人口あたり指標: e-Stat で広島県下市町村の人口を取り、「人口千人あたり収容力」を市町村別に計算 → 散布図上での「過剰／不足」を可視化。
  2. クラスタリング: 災害種別 5 フラグ + BF スコア + log(capacity) を特徴量として k-means / HDBSCAN で施設を類型化し、地図上で色分け。「沿岸特化型」「内陸大規模型」などの自然なグループが浮かぶか検証。
  3. ハザードマップとの突合: DoBoX #46 (津波浸水想定) と緯度経度で空間結合し、「津波対応フラグが立っているのに浸水想定域にある」避難所を抽出。
  4. capacity の欠損補完: 535 件の null を、同一市町村・同一施設種別の中央値で補完して合計収容力を再評価し、補完前後で結論が変わるかを見る。
  5. バリアフリースコアの主成分分析: 11 フラグの相関構造を PCA で抽出し、第 1 主成分が「設備充実度」、第 2 主成分が「想定災害弱者種別」を分離するかを検証。
"""), ] html = render_lesson( num=3, title="避難所4,065件 JSON→DataFrame — 構造探索・共起・分布・地理", tags=["導入", "基礎", "心得", "v2-rewrite"], time="90分", level="リテラシ基礎", data_label='避難所情報 (DoBoX #42)', sections=sections, script_filename="lessons/L03_shelter_analysis.py", ) out = LESSONS / "L03_shelter_analysis.html" out.write_text(html, encoding="utf-8") print(f"\nsaved: {out}") print(f" PNGs: 5 (L03_disaster_counts, L03_jaccard, L03_capacity_dist, " f"L03_muni_coverage, L03_bf_boxplot)") print(f" HTML map: assets/L03_map.html")