""" X02: 埋蔵文化財 × 観光地化潜在性 — 史跡分布から都市プロファイルを読む 広島県の埋蔵文化財包蔵地データ（都城・官衙跡 13件 + その他 198件 = 計211件）と 観光対象施設データ（インフラツーリズム 40件）を用いて、 - 文化財種別ごとの市町分布 - 種別×市町のクロス集計を入力にした **階層クラスタリング**（市町プロファイル） - 同じクロス集計を入力にした **PCA**（主成分分析、2次元圧縮） - 観光化された地域と文化財種別の対応 を分析する **L01 水準** の研究記事用スクリプト。 実行: cd "2026 DoBoX 教材" py -X utf8 lessons/X02_burial_tourism_potential.py """ from pathlib import Path import sys import re import math import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from scipy.cluster.hierarchy import linkage, dendrogram, fcluster from sklearn.decomposition import PCA from sklearn.preprocessing import StandardScaler from _common import ROOT, ASSETS, LESSONS, render_lesson, code, figure plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False # === 0. データ読み込み =================================================== PATH_GOVT = ROOT / "data" / "extras" / "burial_castle_govt.csv" PATH_OTHER = ROOT / "data" / "extras" / "burial_other.csv" PATH_TOUR = ROOT / "data" / "extras" / "infra_tourism.csv" print("=== 0. データ読み込み ===") gv = pd.read_csv(PATH_GOVT, encoding="utf-8-sig") ot = pd.read_csv(PATH_OTHER, encoding="utf-8-sig") tr = pd.read_csv(PATH_TOUR, encoding="utf-8-sig") gv["source_series"] = "都城・官衙跡(#1663)" ot["source_series"] = "その他(#1669)" print(f" 都城・官衙跡: {len(gv)}件 / その他: {len(ot)}件 / 観光対象: {len(tr)}件") # === 1. 埋蔵文化財の統合と種別正規化 ===================================== print("\n=== 1. 埋蔵文化財統合と種別正規化 ===") burial = pd.concat([gv, ot], ignore_index=True) def normalize_kind(s: str) -> str: """生の「種別」列を 11 カテゴリに丸める（複数併記は主要カテゴリ優先）。""" s = str(s).replace("\r", "").replace("\n", "").strip() if "官衙" in s: return "官衙跡" if "寺院" in s: return "寺院跡" if "集落" in s: return "集落跡" if "祭祀" in s: return "祭祀遺跡" if "経塚" in s: return "経塚" if "砂留" in s: return "砂留" if "庭園" in s: return "庭園跡" if "建物" in s or "居宅" in s or "牢" in s or "番所" in s or "市跡" in s: return "建物跡" if "石造" in s: return "石造物" if "埋納" in s: return "埋納地" if "塚" in s and "経" not in s and "里" not in s: return "塚" if "一里塚" in s or "道" in s or "船渠" in s or "船着" in s or "台場" in s or "交通" in s: return "交通遺跡" return "その他" burial["kind11"] = burial["種別"].apply(normalize_kind) burial["city"] = burial["市町名"].fillna("不明").astype(str).str.strip() # 不明・空白行はクラスタリングから外す mask_valid = (burial["city"] != "不明") & (burial["city"] != "") burial_v = burial[mask_valid].copy() print(f" 統合後: {len(burial)}件、市町確定: {len(burial_v)}件") print(f" 種別カテゴリ: {burial_v['kind11'].nunique()} 種") print(f" 市町数: {burial_v['city'].nunique()} 市町") # 中間CSV1: 統合後の生データ（種別正規化付き） burial_v[["番号", "名称", "種別", "kind11", "時代", "city", "所在地", "緯度", "経度", "source_series"]].to_csv( ASSETS / "X02_burial_unified.csv", index=False, encoding="utf-8-sig") # === 2. 種別×市町 クロス集計 ============================================= print("\n=== 2. 種別×市町 クロス集計 ===") crosstab = pd.crosstab(burial_v["city"], burial_v["kind11"]) # 列順を件数降順で安定化 col_order = crosstab.sum(axis=0).sort_values(ascending=False).index.tolist() crosstab = crosstab[col_order] # 行は合計件数昇順（後段デンドログラム可読性） crosstab = crosstab.loc[crosstab.sum(axis=1).sort_values(ascending=False).index] print(f" クロス表サイズ: {crosstab.shape} (市町 × 種別)") print(crosstab.to_string()) # 中間CSV2 crosstab.to_csv(ASSETS / "X02_kind_by_city_crosstab.csv", encoding="utf-8-sig") # === 3. 図1: 種別別件数 + 市町別件数 ==================================== print("\n=== 3. 図1: 種別と市町の単独分布 ===") kind_total = crosstab.sum(axis=0).sort_values(ascending=True) city_total = crosstab.sum(axis=1).sort_values(ascending=True) fig, axes = plt.subplots(1, 2, figsize=(13, 5.4), gridspec_kw={"width_ratios": [1, 1]}) ax = axes[0] COLORS_KIND = { "祭祀遺跡": "#8250df", "経塚": "#bf3989", "その他": "#888888", "砂留": "#0a7e8c", "官衙跡": "#cf222e", "交通遺跡": "#fb8500", "石造物": "#8a5a00", "埋納地": "#a40e26", "庭園跡": "#1f883d", "建物跡": "#0969da", "船渠": "#0550ae", "塚": "#d4a72c", "台場跡": "#888", "寺院跡": "#cf2230", "集落跡": "#1f8838", "寺院跡": "#cf2230", } colors_k = [COLORS_KIND.get(k, "#888") for k in kind_total.index] ax.barh(range(len(kind_total)), kind_total.values, color=colors_k, edgecolor="black", linewidth=0.4) ax.set_yticks(range(len(kind_total))) ax.set_yticklabels(kind_total.index, fontsize=10) for i, v in enumerate(kind_total.values): ax.text(v + 0.5, i, str(int(v)), va="center", fontsize=9) ax.set_xlabel("件数") ax.set_title(f"種別別件数（n={int(kind_total.sum())}件、{len(kind_total)}種別）") ax.grid(axis="x", alpha=0.3) ax = axes[1] ax.barh(range(len(city_total)), city_total.values, color="#0969da", edgecolor="black", linewidth=0.4) ax.set_yticks(range(len(city_total))) ax.set_yticklabels(city_total.index, fontsize=10) for i, v in enumerate(city_total.values): ax.text(v + 0.5, i, str(int(v)), va="center", fontsize=9) ax.set_xlabel("件数") ax.set_title(f"市町別件数（n={len(city_total)}市町）") ax.grid(axis="x", alpha=0.3) plt.tight_layout() plt.savefig(ASSETS / "X02_kind_city_bars.png", dpi=140) plt.close() # === 4. 図2: 緯度経度散布（種別色分け） ================================== print("\n=== 4. 図2: 地理散布（種別色分け） ===") geo = burial_v.dropna(subset=["緯度", "経度"]).copy() geo["lat"] = pd.to_numeric(geo["緯度"], errors="coerce") geo["lon"] = pd.to_numeric(geo["経度"], errors="coerce") geo = geo.dropna(subset=["lat", "lon"]) print(f" 座標確定: {len(geo)}件") fig, ax = plt.subplots(figsize=(10.5, 7.0)) top_kinds = kind_total.sort_values(ascending=False).head(7).index.tolist() for k in top_kinds: sub = geo[geo["kind11"] == k] ax.scatter(sub["lon"], sub["lat"], s=42, alpha=0.78, edgecolor="white", linewidth=0.7, color=COLORS_KIND.get(k, "#888"), label=f"{k} (n={len(sub)})") # その他種別を薄い灰色で other_mask = ~geo["kind11"].isin(top_kinds) if other_mask.sum() > 0: sub = geo[other_mask] ax.scatter(sub["lon"], sub["lat"], s=22, alpha=0.4, color="#bbb", label=f"その他種別 (n={len(sub)})") # 観光施設を×印で重ね描き tr_geo = tr.dropna(subset=["緯度", "経度"]).copy() tr_geo["lat"] = pd.to_numeric(tr_geo["緯度"], errors="coerce") tr_geo["lon"] = pd.to_numeric(tr_geo["経度"], errors="coerce") tr_geo = tr_geo.dropna(subset=["lat", "lon"]) ax.scatter(tr_geo["lon"], tr_geo["lat"], s=80, marker="x", color="#cf222e", linewidth=1.6, label=f"観光対象施設 (n={len(tr_geo)})") ax.set_xlabel("経度") ax.set_ylabel("緯度") ax.set_title("広島県 埋蔵文化財 × 観光対象施設 — 地理分布") ax.legend(loc="lower right", fontsize=9, framealpha=0.95) ax.grid(alpha=0.3) plt.tight_layout() plt.savefig(ASSETS / "X02_geo_scatter.png", dpi=140) plt.close() # === 5. 階層クラスタリング（市町×種別 の行ベクトル）====================== print("\n=== 5. 階層クラスタリング ===") # 入力: クロス集計表（行=市町、列=種別） # 各市町の文化財プロファイルを 1 本のベクトル（13次元）として扱う X = crosstab.values.astype(float) # 行内総和で正規化（市町の規模差を消し、種別構成比に揃える） row_sum = X.sum(axis=1, keepdims=True) row_sum[row_sum == 0] = 1.0 X_norm = X / row_sum # Ward 法で階層クラスタリング（距離行列はユークリッド） Z = linkage(X_norm, method="ward") print(f" 入力: {X_norm.shape} (市町 × 種別構成比)") print(f" リンク行列: {Z.shape}") # 図3: デンドログラム fig, ax = plt.subplots(figsize=(11.5, 5.8)) ddata = dendrogram(Z, labels=list(crosstab.index), ax=ax, color_threshold=0.6 * max(Z[:, 2]), leaf_font_size=10, leaf_rotation=45) ax.set_title("市町の文化財プロファイル — Ward 法による階層クラスタリング") ax.set_ylabel("併合時のばらつき（Ward 距離）") ax.grid(axis="y", alpha=0.3) plt.tight_layout() plt.savefig(ASSETS / "X02_dendrogram.png", dpi=140) plt.close() # 3群に分割（仮説 H4 に従い 3 グループ） K_CLUSTERS = 3 cluster_id = fcluster(Z, t=K_CLUSTERS, criterion="maxclust") city_cluster = pd.Series(cluster_id, index=crosstab.index, name="cluster") print(f" クラスタ件数: {city_cluster.value_counts().sort_index().to_dict()}") print(f" クラスタ別市町:") for c in sorted(city_cluster.unique()): members = city_cluster[city_cluster == c].index.tolist() print(f" cluster{c}: {members}") # 中間CSV3: クラスタ割当 cluster_assign = pd.DataFrame({ "city": crosstab.index, "cluster": cluster_id, "total_burial": crosstab.sum(axis=1).values, }) cluster_assign.to_csv(ASSETS / "X02_city_clusters.csv", index=False, encoding="utf-8-sig") # === 6. PCA（主成分分析） ================================================ print("\n=== 6. PCA（主成分分析） ===") # 標準化（種別ごとの分散を揃える） scaler = StandardScaler() X_std = scaler.fit_transform(X_norm) pca = PCA(n_components=2) PCs = pca.fit_transform(X_std) expl = pca.explained_variance_ratio_ print(f" 寄与率 PC1={expl[0]:.3f}, PC2={expl[1]:.3f}, 累積={expl.sum():.3f}") # ローディング: 各種別が PC1/PC2 に効く強さ loadings = pd.DataFrame(pca.components_.T, index=crosstab.columns, columns=["PC1", "PC2"]) print("ローディング:") print(loadings.round(3).to_string()) # 中間CSV4: PCA座標とローディング pd.DataFrame({"city": crosstab.index, "PC1": PCs[:, 0], "PC2": PCs[:, 1], "cluster": cluster_id}).to_csv( ASSETS / "X02_pca_scores.csv", index=False, encoding="utf-8-sig") loadings.round(4).to_csv(ASSETS / "X02_pca_loadings.csv", encoding="utf-8-sig") # 図4: PCA 散布図（クラスタ色分け + ローディング矢印） fig, ax = plt.subplots(figsize=(10.5, 7.5)) CLUSTER_COLORS = {1: "#0969da", 2: "#cf222e", 3: "#1f883d"} for c in sorted(set(cluster_id)): m = (cluster_id == c) ax.scatter(PCs[m, 0], PCs[m, 1], s=140, alpha=0.85, edgecolor="black", linewidth=0.7, color=CLUSTER_COLORS.get(c, "#888"), label=f"クラスタ{c} (n={int(m.sum())})", zorder=3) for i, name in enumerate(crosstab.index): ax.annotate(name, (PCs[i, 0], PCs[i, 1]), fontsize=9, xytext=(6, 4), textcoords="offset points", zorder=4) # ローディング矢印（スケール調整） arrow_scale = max(np.abs(PCs).max() * 0.55, 1.5) for kind, row in loadings.iterrows(): if abs(row["PC1"]) < 0.15 and abs(row["PC2"]) < 0.15: continue ax.arrow(0, 0, row["PC1"] * arrow_scale, row["PC2"] * arrow_scale, color="#666", alpha=0.55, head_width=0.08, length_includes_head=True, linewidth=1.0, zorder=2) ax.text(row["PC1"] * arrow_scale * 1.15, row["PC2"] * arrow_scale * 1.15, kind, fontsize=9, color="#444", ha="center", bbox=dict(boxstyle="round,pad=0.15", fc="#fff", ec="#bbb", alpha=0.85), zorder=5) ax.axhline(0, color="#888", linewidth=0.5) ax.axvline(0, color="#888", linewidth=0.5) ax.set_xlabel(f"PC1 ({expl[0]*100:.1f}%)") ax.set_ylabel(f"PC2 ({expl[1]*100:.1f}%)") ax.set_title("市町の文化財プロファイル — PCA 2 次元散布 + 階層クラスタ色分け") ax.legend(loc="lower right", fontsize=10) ax.grid(alpha=0.3) plt.tight_layout() plt.savefig(ASSETS / "X02_pca_scatter.png", dpi=140) plt.close() # === 7. 観光対象施設の市町割り当て（lat/lon 最近傍）====================== print("\n=== 7. 観光対象施設の市町割り当て（最近傍埋蔵文化財） ===") # 観光対象施設の lat/lon → 最寄りの埋蔵文化財の市町を割り当て def haversine_km(lat1, lon1, lat2, lon2): R = 6371.0 p1, p2 = math.radians(lat1), math.radians(lat2) dp = math.radians(lat2 - lat1) dl = math.radians(lon2 - lon1) a = (math.sin(dp/2)**2 + math.cos(p1) * math.cos(p2) * math.sin(dl/2)**2) return 2 * R * math.asin(min(1.0, math.sqrt(a))) def assign_city_by_text(row): cities = list(crosstab.index) for col in ["施設名称", "施設の詳細説明", "施設概要"]: s = str(row.get(col, "")) for c in cities: if c in s: return c return None tr_geo["city_text"] = tr_geo.apply(assign_city_by_text, axis=1) geo_arr = geo[["lat", "lon", "city"]].dropna().reset_index(drop=True) def assign_city_by_geo(row): if pd.notnull(row.get("city_text")): return row["city_text"] # 最近傍の埋蔵文化財の市町 dists = geo_arr.apply( lambda r: haversine_km(row["lat"], row["lon"], r["lat"], r["lon"]), axis=1) return geo_arr.loc[dists.idxmin(), "city"] tr_geo["city"] = tr_geo.apply(assign_city_by_geo, axis=1) tour_count = tr_geo.groupby("city").size().rename("tour_count") print(f" 観光施設市町割り当て:") print(tour_count.sort_values(ascending=False).to_string()) # 中間CSV5: 観光施設の市町割り当て結果 tr_geo[["施設名称", "lat", "lon", "分類", "city", "city_text"]].to_csv( ASSETS / "X02_tourism_city_assigned.csv", index=False, encoding="utf-8-sig") # === 8. 文化財密度 vs 観光資源数（市町別の散布 + 単回帰）================= print("\n=== 8. 文化財密度 vs 観光資源数 ===") city_burial = crosstab.sum(axis=1).rename("burial_count") joint = pd.concat([city_burial, tour_count], axis=1).fillna(0).astype(int) joint["cluster"] = city_cluster.reindex(joint.index).fillna(0).astype(int) print(joint.sort_values("burial_count", ascending=False).to_string()) # 単回帰 (numpy.polyfit, deg=1) xs = joint["burial_count"].values.astype(float) ys = joint["tour_count"].values.astype(float) slope, intercept = np.polyfit(xs, ys, 1) ss_tot = ((ys - ys.mean()) ** 2).sum() ss_res = ((ys - (slope * xs + intercept)) ** 2).sum() r2 = 1 - ss_res / ss_tot if ss_tot > 0 else 0.0 r_pearson = np.corrcoef(xs, ys)[0, 1] if len(xs) > 1 else 0.0 print(f" 単回帰: y = {slope:.3f}x + {intercept:.3f}, R² = {r2:.3f}, Pearson r = {r_pearson:.3f}") # 図5: 散布 + 回帰直線 fig, ax = plt.subplots(figsize=(9.5, 6.0)) for c in sorted(joint["cluster"].unique()): m = joint["cluster"] == c if c == 0: continue ax.scatter(joint.loc[m, "burial_count"], joint.loc[m, "tour_count"], s=140, alpha=0.85, edgecolor="black", linewidth=0.7, color=CLUSTER_COLORS.get(int(c), "#888"), label=f"クラスタ{int(c)}", zorder=3) for city in joint.index: ax.annotate(city, (joint.loc[city, "burial_count"], joint.loc[city, "tour_count"]), fontsize=9, xytext=(5, 4), textcoords="offset points") xfit = np.linspace(0, xs.max() * 1.05, 30) ax.plot(xfit, slope * xfit + intercept, "--", color="#cf222e", lw=1.4, label=f"回帰直線 y={slope:.2f}x+{intercept:.2f}\nR²={r2:.3f}, r={r_pearson:.3f}") ax.set_xlabel(f"市町別 埋蔵文化財件数（{int(crosstab.values.sum())}件 / {len(crosstab)}市町）") ax.set_ylabel("市町別 観光対象施設件数（40件）") ax.set_title("文化財密度 vs 観光資源数 — 市町散布 + 単回帰") ax.legend(loc="upper left", fontsize=10) ax.grid(alpha=0.3) plt.tight_layout() plt.savefig(ASSETS / "X02_burial_vs_tourism.png", dpi=140) plt.close() # 中間CSV6: 市町別 文化財数 × 観光数 結合 joint.reset_index().to_csv( ASSETS / "X02_burial_vs_tour_by_city.csv", index=False, encoding="utf-8-sig") # === 9. クラスタ別 観光化率 + 種別構成 =================================== print("\n=== 9. クラスタ別 観光化率と種別構成 ===") cluster_summary_rows = [] for c in sorted(set(cluster_id)): cities_in_c = city_cluster[city_cluster == c].index.tolist() sub_burial = crosstab.loc[cities_in_c] n_burial = int(sub_burial.values.sum()) n_tour = int(tour_count.reindex(cities_in_c).fillna(0).sum()) # 「その他」を除いた上位2 種別を主種別とする sums = sub_burial.sum(axis=0) sums_no_other = sums[sums.index != "その他"] top_kinds = sums_no_other.sort_values(ascending=False).head(2) top_kind_str = "、".join(f"{k}({int(v)})" for k, v in top_kinds.items() if v > 0) cluster_summary_rows.append({ "cluster": int(c), "n_cities": len(cities_in_c), "cities": "、".join(cities_in_c), "n_burial": n_burial, "n_tourism": n_tour, "tour_per_burial": round(n_tour / max(n_burial, 1), 3), "dominant_kinds": top_kind_str, }) cluster_summary = pd.DataFrame(cluster_summary_rows) print(cluster_summary.to_string(index=False)) # 図6: クラスタ別 種別構成 + 観光化率 fig, axes = plt.subplots(1, 2, figsize=(13, 5.4), gridspec_kw={"width_ratios": [1.4, 1]}) ax = axes[0] # クラスタ別 種別構成（積み上げ棒） ratio = pd.DataFrame() for c in sorted(set(cluster_id)): cities_in_c = city_cluster[city_cluster == c].index.tolist() s = crosstab.loc[cities_in_c].sum(axis=0) ratio[f"C{c}"] = s / s.sum() ratio = ratio.T # 行=クラスタ、列=種別 bottom = np.zeros(len(ratio)) xpos = np.arange(len(ratio)) for kind in ratio.columns: vals = ratio[kind].values ax.bar(xpos, vals, bottom=bottom, label=kind, color=COLORS_KIND.get(kind, "#888"), edgecolor="white", linewidth=0.5) bottom += vals ax.set_xticks(xpos) ax.set_xticklabels(ratio.index, fontsize=11) ax.set_ylabel("種別構成比") ax.set_ylim(0, 1) ax.set_title("クラスタ別 文化財 種別構成（割合）") ax.legend(loc="center left", bbox_to_anchor=(1.0, 0.5), fontsize=8, ncol=1, framealpha=0.95) ax.grid(axis="y", alpha=0.3) # 観光化率（単位文化財件数あたりの観光施設数） ax = axes[1] tpb = cluster_summary["tour_per_burial"].values ax.bar(range(len(cluster_summary)), tpb, color=[CLUSTER_COLORS[int(c)] for c in cluster_summary["cluster"]], edgecolor="black", linewidth=0.5) ax.set_xticks(range(len(cluster_summary))) ax.set_xticklabels([f"C{int(c)}" for c in cluster_summary["cluster"]], fontsize=11) for i, v in enumerate(tpb): ax.text(i, v + 0.005, f"{v:.3f}", ha="center", fontsize=10) ax.set_ylabel("観光対象施設 ÷ 埋蔵文化財件数") ax.set_title("クラスタ別 観光化率") ax.grid(axis="y", alpha=0.3) plt.tight_layout() plt.savefig(ASSETS / "X02_cluster_compare.png", dpi=140) plt.close() # 中間CSV7 cluster_summary.to_csv(ASSETS / "X02_cluster_summary.csv", index=False, encoding="utf-8-sig") # === 10. HTML 用テーブル ================================================ print("\n=== 10. HTML テーブル生成 ===") crosstab_top = crosstab.head(8) # 上位8市町だけ表示 crosstab_html = crosstab_top.to_html() kind_total_html = (kind_total.sort_values(ascending=False) .rename("件数").to_frame().to_html()) city_total_html = (city_total.sort_values(ascending=False) .rename("件数").to_frame().to_html()) cluster_summary_html = cluster_summary.to_html(index=False) loadings_html = loadings.round(3).to_html() joint_html = (joint.sort_values("burial_count", ascending=False) .head(10).to_html()) # 実測値で本文に埋め込む数値を確定 N_TOTAL = int(len(burial)) N_VALID = int(len(burial_v)) N_KINDS = int(crosstab.shape[1]) N_CITIES = int(crosstab.shape[0]) PC1_PCT = float(expl[0] * 100) PC2_PCT = float(expl[1] * 100) PC_CUM_PCT = float(expl.sum() * 100) R2_VALUE = float(r2) R_VALUE = float(r_pearson) SLOPE_VALUE = float(slope) INTERCEPT_VALUE = float(intercept) # 1件追跡表（要件K）: 「地蔵河原一里塚」を例にする trace_rows = [ ("0. 元データ", "burial_other.csv の 1 行（番号=167）", "名称='地蔵河原一里塚', 種別='交通遺跡', 時代='近世', " "市町名='広島市', 緯度=34.530, 経度=132.498", "1×13 列"), ("1. 種別正規化", "normalize_kind('交通遺跡') が '交通遺跡' に丸める", f"kind11='交通遺跡'（{N_KINDS}カテゴリの1つ）", "1×1"), ("2. クロス集計", "(city='広島市', kind11='交通遺跡') のセルに +1", f"crosstab.loc['広島市','交通遺跡'] = {int(crosstab.loc['広島市','交通遺跡'])}", f"{N_CITIES}行×{N_KINDS}列の表"), ("3. 行内正規化", "広島市の行を「広島市の総件数」で割って構成比に", f"X_norm['広島市'] の '交通遺跡' 比率 = {int(crosstab.loc['広島市','交通遺跡'])}/{int(crosstab.loc['広島市'].sum())} = {int(crosstab.loc['広島市','交通遺跡'])/max(int(crosstab.loc['広島市'].sum()),1):.3f}", f"{N_CITIES}×{N_KINDS} の比率行列"), ("4. 階層クラスタ", f"Ward 法で {N_CITIES} 市町を 3 群に併合", f"広島市 → cluster={int(city_cluster.loc['広島市'])}", "1市町 → 1 クラスタ番号"), ("5. PCA 圧縮", "種別ごとの分散を揃え、2 主成分に圧縮", f"広島市の (PC1, PC2) = ({PCs[list(crosstab.index).index('広島市'),0]:.3f}, {PCs[list(crosstab.index).index('広島市'),1]:.3f})", "1市町 → (PC1, PC2)"), ("6. 観光化率", f"広島市は観光対象 {int(tour_count.get('広島市', 0))} 件 / 文化財 {int(crosstab.loc['広島市'].sum())} 件", f"tour_per_burial = {int(tour_count.get('広島市', 0))/max(int(crosstab.loc['広島市'].sum()),1):.3f}", "1 市町 → 比率"), ] trace_html = ( "" "" + "".join(f"" for a, b, c, d in trace_rows) + "
段階このデータで何が起きるか結果サイズ
{a}{b}{c}{d}
" ) # === 11. HTML レンダリング =============================================== print("\n=== 11. HTML レンダリング ===") sections = [ ("学習目標と問い", f"""

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

「広島県の埋蔵文化財は、どの市町に・どんな種別で偏在し、観光資源化されている／されていない地域を地図と数値で見分けられるか？」

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

立てた仮説

  1. H1（種別×地理）: 文化財種別ごとに地理分布が違う。都城・官衙跡は古代の行政中心地（広島市・福山市・府中市）に集中、 集落跡・経塚は山間部に広く分布する。
  2. H2（観光化の偏り）: 観光対象施設は 都城・古墳・建造物系の大型遺構が多く、経塚・砂留などの小規模文化財は観光化されにくい。
  3. H3（密度相関）: 市町別の 文化財件数と 観光資源数は正の相関を持つ（多い所は両方多い）。
  4. H4（市町プロファイルの3群）: {N_KINDS} 種別の市町別件数をベクトル化して階層クラスタリングすると、 市町は 都市型（古代行政＋近代観光）／農村型（祭祀・経塚分散）／遺跡密集型（特定種別が突出）の 3 群に分かれる。
  5. H5（PCA で2軸抽出）: 同じデータを PCA で 2 次元に圧縮すると、 PC1 が 「古代行政 vs 近世生活」、PC2 が 「祭祀・宗教 vs 交通・産業」のような主題軸を表すはず。

到達点

"""), ("使用データ", f"""

※ 11 シリーズのうち 2 シリーズ（都城・官衙跡、その他）の合計 {N_TOTAL} 件を使用。 他 9 シリーズ（古墳・横穴／貝塚／集落跡・散布地／城館跡／社寺跡／生産遺跡／その他の墳墓／近代以降の単独遺跡／水中遺跡）は 発展課題で網羅予定。本レッスンでは種別正規化により {N_KINDS} カテゴリ に丸めるため、 データ範囲 ＝ 11 シリーズの直交分類ではなく、種別文字列ベースの分類になる点に注意。

"""), ("ダウンロード（再現用データ・中間データ・図）", f"""

本レッスンの全成果物に直リンク。学習者は途中ステップから再現できる。

1. 生データ（DoBoX 由来）

ファイル形式行数取得元
data/extras/burial_castle_govt.csv CSV13DoBoX #1663
data/extras/burial_other.csv CSV198DoBoX #1669
data/extras/infra_tourism.csv CSV40DoBoX #1280
data/dataset_index.csv CSV551DoBoX 全件カタログ

2. プログラムが生成する中間データ

ファイル内容使う分析
X02_burial_unified.csv 埋蔵文化財統合（{N_TOTAL}件）+ 種別正規化列 kind11分析2 以降の入力
X02_kind_by_city_crosstab.csv 市町×種別 クロス集計（{N_CITIES}×{N_KINDS}）分析3 / 分析4 / 分析5
X02_city_clusters.csv 各市町のクラスタ番号（Ward法 / k=3）分析4
X02_pca_scores.csv 各市町の PC1 / PC2 座標分析5
X02_pca_loadings.csv {N_KINDS} 種別の PC1 / PC2 ローディング分析5
X02_tourism_city_assigned.csv 観光施設 40 件の市町割り当て結果分析6
X02_burial_vs_tour_by_city.csv 市町別 文化財数 × 観光数 結合表分析6
X02_cluster_summary.csv クラスタ別 件数・観光化率・主種別分析6

3. 図 PNG

4. 再現スクリプト

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

スクリプト本体: lessons/X02_burial_tourism_potential.py

5. 1 件のデータが分析を通る過程（要件K）

例: 「地蔵河原一里塚」（広島市安佐北区可部、近世の交通遺跡、burial_other.csv #167 行）。 このレッスンの全 6 段階を 1 件で追う:

""" + trace_html), ("分析1: 種別と市町の単独分布", f"""

狙い

埋蔵文化財 {N_TOTAL} 件を「種別ごと」「市町ごと」に集計し、それぞれの分布の偏りを見る。 クロス集計に進む前の入口として、各次元の サイズ感 を掴む（仮説 H1 の前提確認）。

手法（種別文字列の正規化 + 件数集計）

実装

""" + code(''' gv = pd.read_csv("data/extras/burial_castle_govt.csv", encoding="utf-8-sig") ot = pd.read_csv("data/extras/burial_other.csv", encoding="utf-8-sig") burial = pd.concat([gv, ot], ignore_index=True) def normalize_kind(s): s = str(s).replace("\\r","").replace("\\n","").strip() if "官衙" in s: return "官衙跡" if "祭祀" in s: return "祭祀遺跡" if "経塚" in s: return "経塚" if "砂留" in s: return "砂留" if "庭園" in s: return "庭園跡" if "石造" in s: return "石造物" if "埋納" in s: return "埋納地" if "塚" in s and "経" not in s and "里" not in s: return "塚" if "一里塚" in s or "道" in s or "船" in s or "台場" in s or "交通" in s: return "交通遺跡" return "その他" burial["kind11"] = burial["種別"].apply(normalize_kind) burial["city"] = burial["市町名"].fillna("不明").astype(str).str.strip() ''') + f"""

結果（図と読み取り）

なぜこの図か: 種別と市町の 件数の相対比較 を一目で見たい。 両方を 1 枚の左右並びにすると 「どの種別が多いか／どの市町が多いか」 が同時に対比できる。

""" + figure("assets/X02_kind_city_bars.png", "図1: 種別別件数（左）と市町別件数（右）") + f"""

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

結果（表と読み取り）

表1: 種別別件数（{N_KINDS}カテゴリ）

""" + kind_total_html + f"""

下位の「塚」「庭園跡」は数件しかなく — クラスタリングや PCA で 影響度が低い 列となる。

表2: 市町別件数（{N_CITIES}市町）

""" + city_total_html + """

福山市が断トツ。中山間部の市町は件数が少なく、後段クラスタリングで 独立クラスタ になりやすい。

"""), ("分析2: 種別 × 市町 クロス集計", f"""

狙い

「どの市町に・どの種別が何件あるか」を 1 枚の表にする。 これが以降の階層クラスタリング・PCA の 共通入力 になる重要ステップ。

手法（pandas crosstab）

このレッスンの中心テーブル: ここで作る {N_CITIES}×{N_KINDS} の表が、分析4（階層クラスタリング）と分析5（PCA）の入力。 この表を「市町ごとの{N_KINDS}個の数字の並び（=プロファイル）」と読むと、 以降の手法は 「{N_CITIES} 個のプロファイルが似ているもの同士をグループ化」「{N_CITIES} 個を 2 次元に圧縮」と言える。

実装

""" + code(''' crosstab = pd.crosstab(burial["city"], burial["kind11"]) # 列順を全体件数の降順で安定化 crosstab = crosstab[crosstab.sum(axis=0).sort_values(ascending=False).index] # 行も総件数降順 crosstab = crosstab.loc[crosstab.sum(axis=1).sort_values(ascending=False).index] print(crosstab.shape) # (16, 13) ''') + f"""

結果（表と読み取り）

なぜ表だけか: ここでは「数字を並べた表そのもの」が成果物。 ヒートマップ化すると後の分析と重複するので、まず数値を確認する。

表3: 種別×市町 クロス集計（上位8市町を抜粋, 行=市町、列=種別）

""" + crosstab_html + f"""

※ 表は上位 8 市町だけ抜粋。実際は {N_CITIES} 市町。完全版は X02_kind_by_city_crosstab.csv

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

"""), ("分析3: 地理散布で見る種別の空間分布", """

狙い

クロス集計（数字の表）だけでは「西部寄り／東部寄り」「沿岸／山間」が見えない。 緯度経度の散布図に種別色分けで重ね描きし、空間的偏りを目で確認する（仮説 H1 の検証）。

手法（緯度経度散布、種別色分け、観光対象重ね描き）

実装

""" + code(''' geo = burial.dropna(subset=["緯度","経度"]).copy() geo["lat"] = pd.to_numeric(geo["緯度"], errors="coerce") geo["lon"] = pd.to_numeric(geo["経度"], errors="coerce") geo = geo.dropna(subset=["lat","lon"]) fig, ax = plt.subplots(figsize=(10.5, 7)) top_kinds = kind_total.sort_values(ascending=False).head(7).index.tolist() for k in top_kinds: sub = geo[geo["kind11"] == k] ax.scatter(sub["lon"], sub["lat"], s=42, alpha=0.78, label=k) # 観光対象を ×印で重ね描き ax.scatter(tr["経度"], tr["緯度"], s=80, marker="x", color="red", label="観光対象施設") ax.set_xlabel("経度"); ax.set_ylabel("緯度") ''') + """

結果（図と読み取り）

なぜこの図か: 種別ごとの空間配置と観光対象の位置関係を 1 枚で同時に見たい。 種別を色、観光を ×印にすると2系列が干渉せずに重ねられる。

""" + figure("assets/X02_geo_scatter.png", "図2: 緯度経度散布（種別色分け + 観光対象 ×印）") + """

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

"""), ("分析4: 階層クラスタリング — 市町を文化財プロファイルでグループ化", f"""

狙い

「{N_CITIES} 市町を、文化財の種別構成が似ている者同士でまとめると、何個のグループに分かれるか？」 仮説 H4（3 群仮説）を検証する。

STEP の全体像（複数手法を組み合わせる）

STEP担当入力出力
STEP 1行内正規化（各市町の合計件数で割る） {N_CITIES}×{N_KINDS} の整数表{N_CITIES}×{N_KINDS} の比率表（行合計=1.0）
STEP 2Ward 法で階層クラスタリング {N_CITIES}×{N_KINDS} の比率表系統樹（デンドログラム） + 各市町のクラスタ番号

STEP 1: なぜ「件数そのもの」ではなく「比率」を入力にするのか

福山市は {int(city_total.get('福山市',0))} 件、安芸太田町は {int(city_total.get('安芸太田町',0))} 件。件数を直接使うと、福山市が他の市町と何もかも違う「特異点」として孤立してしまう。 本当に見たいのは 「文化財の種別の混ざり具合」であって、件数の絶対値ではない。 そこで各市町の行をその市町の合計で割って比率に変換する（市町の規模差をならす）。

STEP 1 の Before/After 例（広島市の行）:

段階祭祀遺跡その他経塚砂留官衙跡交通遺跡...合計
Before（件数）{int(crosstab.loc['広島市'].get('祭祀遺跡', 0))}{int(crosstab.loc['広島市'].get('その他', 0))}{int(crosstab.loc['広島市'].get('経塚', 0))}{int(crosstab.loc['広島市'].get('砂留', 0))}{int(crosstab.loc['広島市'].get('官衙跡', 0))}{int(crosstab.loc['広島市'].get('交通遺跡', 0))}...{int(crosstab.loc['広島市'].sum())}
After（比率）{crosstab.loc['広島市'].get('祭祀遺跡',0)/max(crosstab.loc['広島市'].sum(),1):.3f}{crosstab.loc['広島市'].get('その他',0)/max(crosstab.loc['広島市'].sum(),1):.3f}{crosstab.loc['広島市'].get('経塚',0)/max(crosstab.loc['広島市'].sum(),1):.3f}{crosstab.loc['広島市'].get('砂留',0)/max(crosstab.loc['広島市'].sum(),1):.3f}{crosstab.loc['広島市'].get('官衙跡',0)/max(crosstab.loc['広島市'].sum(),1):.3f}{crosstab.loc['広島市'].get('交通遺跡',0)/max(crosstab.loc['広島市'].sum(),1):.3f}...1.000

STEP 2: 階層クラスタリング（Ward 法）はツールとして何をするか

役割: {N_CITIES} 個の比率の並び（={N_CITIES} 個の市町プロファイル）を、似ている者同士で順に併合し、系統樹を作る。

動作のイメージ:

  1. 最初は {N_CITIES} 市町がそれぞれ単独グループ（={N_CITIES} グループ）
  2. 「最も似ている 2 グループ」を 1 つに併合する → {N_CITIES-1} グループ
  3. これを繰り返して、最後は 1 グループ（全部まとまる）
  4. 併合の履歴を 系統樹 として描くと、枝の高さが「違いの大きさ」

このツールの入出力:

STEP 1 + STEP 2 の実装

""" + code(f''' from scipy.cluster.hierarchy import linkage, dendrogram, fcluster X = crosstab.values.astype(float) # {N_CITIES}×{N_KINDS} 整数表 row_sum = X.sum(axis=1, keepdims=True); row_sum[row_sum==0] = 1.0 X_norm = X / row_sum # STEP 1: 行内比率に正規化 Z = linkage(X_norm, method="ward") # STEP 2-A: 系統樹（リンク行列） dendrogram(Z, labels=list(crosstab.index)) # STEP 2-B: 描画 cluster_id = fcluster(Z, t=3, criterion="maxclust") # 3 群に切る ''') + """

結果（図と読み取り）

なぜこの図か: 「何市町が・どの段階で・どれくらい違うのか」を 系統樹で見たい。 枝の高さが「ばらつき＝違いの大きさ」を表すので、どこで切れば自然な群分けになるかが目で判断できる。

""" + figure("assets/X02_dendrogram.png", "図3: 市町プロファイルの系統樹（Ward 法）") + f"""

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

結果（表と読み取り）

表4: 各クラスタの構成市町と特徴

""" + cluster_summary_html + f"""

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

"""), (f"分析5: PCA — {N_KINDS} 種別を 2 次元に圧縮", f"""

狙い

「{N_KINDS} 個の数字の並び（市町プロファイル）を、2 つの数字に要約できないか？」 仮説 H5 を検証。圧縮した 2 次元の意味を ローディングで読む。

このツールは何をするか（黒箱でOK）

PCA = 主成分分析（Principal Component Analysis）。 たくさんの列を持つ表（=多次元データ）を「いちばん情報を残す方向」を新しい軸として 2 本だけ選び、表を 2 列に圧縮するツール。

入出力:

用語の平易な言い換え（要件P）

専門用語平易な言い方このレッスンでの具体
主成分（PC1, PC2）新しく作った 2 本の軸 {N_CITIES} 市町を散布図に置くための X 軸 / Y 軸
ローディング各種別がその新しい軸にどれくらい効くか（影響の重み） 「祭祀遺跡は PC1 にプラス {float(loadings.loc['祭祀遺跡','PC1']):+.2f} 効く」など
寄与率（explained variance ratio）その軸が元データの情報のうち何割を残すか PC1={PC1_PCT:.0f}%、PC2={PC2_PCT:.0f}% → 2 軸合計で {PC_CUM_PCT:.0f}% 保持
標準化（StandardScaler）列ごとに「平均0・分散1」に揃える前処理 件数が多い種別が PCA を支配しないようにする

パラメータと限界

実装

""" + code(f''' from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA # X_norm = 分析4 で作った {N_CITIES}×{N_KINDS} 比率表 X_std = StandardScaler().fit_transform(X_norm) # 列ごとに平均0,分散1 pca = PCA(n_components=2) PCs = pca.fit_transform(X_std) # {N_CITIES}×2 の座標 expl = pca.explained_variance_ratio_ # [PC1寄与率, PC2寄与率] # ローディング: どの種別が PC1/PC2 にどれくらい効くか loadings = pd.DataFrame(pca.components_.T, index=crosstab.columns, columns=["PC1","PC2"]) ''') + f"""

結果（図と読み取り）

なぜこの図か: {N_CITIES} 市町を 2 次元散布に置き、クラスタ色分けと ローディング矢印を重ねたい。 矢印は「どの種別がどの方向を指すか」を示し、軸の意味を直感的に解釈できる。

""" + figure("assets/X02_pca_scatter.png", "図4: PCA 2 次元散布 — 市町をクラスタ色で表示、矢印が各種別のローディング") + f"""

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

結果（表と読み取り）

表5: PCA ローディング（{N_KINDS} 種別 × 2 主成分）

""" + loadings_html + f"""

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

"""), ("分析6: 文化財密度 vs 観光資源数", f"""

狙い

「文化財が多い市町は観光地化も進んでいるのか？」 H3（正の相関）を単回帰で検証。さらに、クラスタごとの観光化率を比較して H2 を補強。

手法（観光施設の市町割り当て + 単回帰）

STEP A の Before/After 例（紅葉谷川庭園砂防施設の場合）:

段階処理結果
Beforeinfra_tourism.csv の 1 行（緯度=34.294, 経度=132.324, 説明=「弥山から嚴島神社の背後…」）市町列が空
STEP A-1（テキスト）16 市町名で str.contains → 「廿日市市」「広島市」がヒット？ → 不明テキスト割り当て=None
STEP A-2（最近傍）緯度経度から最寄りの埋蔵文化財を引く → 廿日市市の遺跡が最近city='廿日市市'

実装

""" + code(''' import math def haversine_km(lat1, lon1, lat2, lon2): R = 6371.0 p1, p2 = math.radians(lat1), math.radians(lat2) dp = math.radians(lat2 - lat1); dl = math.radians(lon2 - lon1) a = math.sin(dp/2)**2 + math.cos(p1)*math.cos(p2)*math.sin(dl/2)**2 return 2 * R * math.asin(min(1.0, math.sqrt(a))) # STEP A: 観光施設の市町割り当て def assign_city(row): cities = list(crosstab.index) for col in ["施設名称","施設の詳細説明","施設概要"]: s = str(row.get(col,"")) for c in cities: if c in s: return c # フォールバック: 最近傍埋蔵文化財の市町 dists = geo.apply(lambda r: haversine_km(row["lat"], row["lon"], r["lat"], r["lon"]), axis=1) return geo.loc[dists.idxmin(), "city"] tour_count = tr_geo.groupby(tr_geo.apply(assign_city, axis=1)).size() # STEP B+C: 散布 + 単回帰 joint = pd.concat([crosstab.sum(axis=1), tour_count], axis=1).fillna(0).astype(int) slope, intercept = np.polyfit(joint["burial_count"], joint["tour_count"], 1) r = np.corrcoef(joint["burial_count"], joint["tour_count"])[0,1] ''') + f"""

結果（図と読み取り）

なぜこの図か: 2 変数の関係（散布図）+ 全体傾向（単回帰直線）+ クラスタ所属（色）を 1 枚で同時に見たい。 回帰直線から外れる市町は「文化財数の割に観光化が進んでいる／いない」異常点として識別できる。

""" + figure("assets/X02_burial_vs_tourism.png", "図5: 市町別 文化財件数 vs 観光資源数 + 単回帰直線") + f"""

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

結果（表と読み取り）

表6: 市町別 文化財件数 × 観光資源数（上位10）

""" + joint_html + """

図6: クラスタ別 種別構成（左）と観光化率（右）

""" + figure("assets/X02_cluster_compare.png", "図6: クラスタ別 種別構成と観光化率の比較") + """

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

"""), ("仮説検証と考察", f"""

仮説と結果の照合

#仮説判定根拠
H1文化財種別ごとに地理分布が違う（都城・官衙跡は行政中心地、集落跡は山間部） 支持 図2 の地理散布で東部に砂留、北部に祭祀遺跡、府中市に官衙跡が集中する明確なパターン。 分析2 のクロス集計でも府中市の「官衙跡」{int(crosstab.loc['府中市','官衙跡']) if '府中市' in crosstab.index and '官衙跡' in crosstab.columns else 0} 件が突出。
H2観光対象は都城・古墳系が多く、集落跡や経塚は少ない 部分支持 図6（クラスタ別観光化率）でクラスタ間に明確な差。 古代行政型クラスタ（府中市・府中町）は観光化率 0%、複合型クラスタが {float(cluster_summary.loc[cluster_summary['cluster']==3,'tour_per_burial'].values[0]) if 3 in cluster_summary['cluster'].values else 0:.3f}。 ただし「都城・古墳系が多い」と予想した行政型クラスタが むしろ観光化されていない 結果は予想と逆。
H3市町別 文化財件数 と 観光資源数 が正の相関 反証寄り 図5 の単回帰で R² = {R2_VALUE:.3f}, Pearson r = {R_VALUE:.3f} と相関は極めて弱い。 福山市（文化財多・観光少）と呉市（文化財少・観光多）など外れ値が多い。 生件数の比例より種別構成・地域歴史が観光化を決める方が支配的。
H4市町プロファイル → 階層クラスタリングで 3 群に分かれる 支持 図3 のデンドログラムで 3 主枝が自然に出現。 各クラスタの主種別と地理特性が解釈可能（複合型 / 古代行政型 / 中世宗教型）。 ただし事前ラベル「都市型・農村型・遺跡密集型」とは違う 3 群構造になった。
H5PCA 2 軸に主題的な意味（古代行政 vs 中世宗教 など）が見える 部分支持 図4 とローディング表（表5）で PC1 = 「祭祀的小規模遺構 vs 行政遺構」、 PC2 = 「近世近代建築 vs 江戸期生活遺構」と解釈できる。 累積寄与率 {PC_CUM_PCT:.0f}% で主構造の半分弱を保持（残りは PC3 以降に分散）。

考察

"""), ("発展課題（結果から導かれる新たな問い）", f"""

結果X → 新仮説Y → 課題Z の 3 段で書く。

  1. 11 シリーズ完全網羅で再分析
  2. 観光化判定の精緻化（マルチソース）
  3. 距離的「観光化潜在性」マップ
  4. 時代軸を加えた 3 次元クラスタリング
  5. 判別分析（LDA）で「観光化される文化財」の予測
"""), ] html = render_lesson( num=2, # 表示は X02 title=f"埋蔵文化財 × 観光地化潜在性 — 階層クラスタリングと PCA で広島県 {N_CITIES} 市町を読む", tags=["X-tier", "draft", "L01水準", "PCA", "階層クラスタ"], time="120分", level="リテラシ研究", data_label=( '#1663 都城・官衙跡 13件、' '#1669 その他 190件、' '#1280 インフラツーリズム 40件' ), sections=sections, script_filename="lessons/X02_burial_tourism_potential.py", ) # X-tier の draft 属性を body に挿入 html = html.replace('
', '
') out = LESSONS / "X02_burial_tourism_potential.html" out.write_text(html, encoding="utf-8") print(f"\nsaved: {out}") print("PNGs:") for p in sorted(ASSETS.glob("X02_*.png")): print(f" - {p.name} ({p.stat().st_size//1024} KB)") print("CSVs:") for p in sorted(ASSETS.glob("X02_*.csv")): print(f" - {p.name} ({p.stat().st_size//1024} KB)")