"""
L02: 県内カメラ351点 — folium 多層マップと点配置の空間統計

DoBoX #1279 の県内カメラ 351 台を、

  1. folium MarkerCluster + HeatMap 2 レイヤで地図化
  2. k-NN 距離分布 (k=1,3,5) を haversine で計算しヒストグラム化
  3. Voronoi 担当面積分布で「カメラの守備範囲のばらつき」を測る
  4. 同領域に Poisson 配置を生成し、最近傍距離分布を実データと比較 (KS検定)
  5. 管理区分 (国/県/市/その他) 別に空間集中度 (NN平均距離) を比較

の 5 ステップで深掘りする v2 教材。

実行:
    cd "2026 DoBoX 教材"
    py -X utf8 lessons/L02_camera_map.py

データは存在しなければ DoBoX #1279 から自動取得する (ensure_dataset)。
"""
from pathlib import Path
import numpy as np
import pandas as pd
import folium
from folium.plugins import MarkerCluster, HeatMap
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon as MplPolygon
from matplotlib.collections import PatchCollection
from scipy.spatial import Voronoi, cKDTree
from scipy.stats import ks_2samp
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. データ自動取得 =======================================================
CSV_PATH = ROOT / "data" / "camera_list.csv"
print("=== 0. データ自動取得 ===")
ensure_dataset(CSV_PATH, dataset_id=1279, label="県内カメラ #1279")

# === 1. 読込 & 前処理 =========================================================
print("\n=== 1. 読込 & 前処理 ===")
df = pd.read_csv(CSV_PATH)
df = df.dropna(subset=["緯度", "経度"]).reset_index(drop=True)
print(f"件数: {len(df)}台")
print(df["管理区分"].value_counts().to_string())

LAT = df["緯度"].to_numpy()
LON = df["経度"].to_numpy()

# 管理者カテゴリの粗分類 (国/県/市/その他)
def classify_admin(s):
    s = str(s)
    if "国土交通省" in s:
        return "国"
    if s.endswith("県") or s == "広島県":
        return "県"
    if s.endswith("市"):
        return "市"
    return "その他"


df["admin_cat"] = df["所管"].map(classify_admin)
print("\n管理者粗分類:")
print(df["admin_cat"].value_counts().to_string())

# === 2. 図 1 + folium HTML: MarkerCluster + HeatMap 2 レイヤ ==================
print("\n=== 2. folium 多層マップ ===")
center = [LAT.mean(), LON.mean()]
m = folium.Map(location=center, zoom_start=9, tiles="OpenStreetMap",
               control_scale=True)

# レイヤ A: MarkerCluster (管理区分で色分け)
cluster_layer = folium.FeatureGroup(name="MarkerCluster (区分別)", show=True)
mc = MarkerCluster().add_to(cluster_layer)
COLOR = {"道路": "blue", "河川": "green", "ため池": "red",
         "海岸": "orange", "港湾": "purple", "その他": "gray"}
for _, r in df.iterrows():
    folium.CircleMarker(
        [r["緯度"], r["経度"]], radius=4,
        color=COLOR.get(r["管理区分"], "gray"),
        fill=True, fill_opacity=0.85, weight=1,
        popup=folium.Popup(
            f"<b>{r['カメラ名']}</b><br>"
            f"住所: {r['住所']}<br>"
            f"区分: {r['管理区分']} / 所管: {r['所管']}<br>"
            f"<a href='{r['公開URL']}' target='_blank'>ライブ映像</a>",
            max_width=320),
    ).add_to(mc)
cluster_layer.add_to(m)

# レイヤ B: HeatMap (密度)
heat_layer = folium.FeatureGroup(name="HeatMap (密度)", show=False)
HeatMap(list(zip(LAT, LON)), radius=18, blur=22,
        min_opacity=0.3).add_to(heat_layer)
heat_layer.add_to(m)

folium.LayerControl(collapsed=False).add_to(m)

legend = """
<div style='position:fixed; top:60px; right:20px; z-index:9999; background:white;
padding:8px 12px; border:1px solid #888; font-size:12px; line-height:1.6;
border-radius:4px; box-shadow:0 1px 4px rgba(0,0,0,0.2);'>
<b>管理区分</b><br>
<span style='color:blue'>●</span> 道路 <span style='color:green'>●</span> 河川<br>
<span style='color:red'>●</span> ため池 <span style='color:orange'>●</span> 海岸<br>
<span style='color:purple'>●</span> 港湾 <span style='color:gray'>●</span> その他<br>
<span class="tnote">右上のレイヤ切替で HeatMap も表示可</span>
</div>"""
m.get_root().html.add_child(folium.Element(legend))
MAP_PATH = ASSETS / "L02_map.html"
m.save(str(MAP_PATH))
print(f"  saved: {MAP_PATH}")

# === 3. 図 2: 管理区分別カメラ数 + 緯度経度ジョイント散布 ====================
print("\n=== 3. 図 2: 管理区分別 + 散布 ===")
fig, axes = plt.subplots(1, 2, figsize=(11.5, 4.6),
                         gridspec_kw={"width_ratios": [1, 1.3]})
vc = df["管理区分"].value_counts()
bars = axes[0].bar(range(len(vc)), vc.values,
                   color=[COLOR.get(k, "gray") for k in vc.index],
                   edgecolor="black", linewidth=0.7)
axes[0].set_xticks(range(len(vc)))
axes[0].set_xticklabels(vc.index, rotation=0)
axes[0].set_ylabel("カメラ数")
axes[0].set_title(f"管理区分別カメラ数 (全 {len(df)} 台)")
for i, v in enumerate(vc.values):
    axes[0].text(i, v + 1.5, str(int(v)), ha="center", fontsize=10)
axes[0].grid(alpha=0.3, axis="y")

# 散布: 区分別色, 海岸線がうっすら見えるように
for cat, c in COLOR.items():
    sub = df[df["管理区分"] == cat]
    axes[1].scatter(sub["経度"], sub["緯度"], s=20, c=c,
                    alpha=0.75, edgecolors="black", linewidth=0.3, label=cat)
axes[1].set_xlabel("経度 (東経)")
axes[1].set_ylabel("緯度 (北緯)")
axes[1].set_title("カメラ配置 (経度×緯度)")
axes[1].legend(fontsize=9, loc="lower right", ncol=2)
axes[1].set_aspect(1.2)  # 緯度1度 ≈ 1.2 × 経度1度 (北緯34.5)
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.savefig(ASSETS / "L02_overview.png", dpi=140)
plt.close()


# === 4. 図 3: 最近傍距離 (k-NN) ===============================================
def haversine_km(lat1, lon1, lat2, lon2):
    """ベクトル化 haversine 距離 (km)。入力は度。"""
    R = 6371.0
    p1 = np.radians(lat1)
    p2 = np.radians(lat2)
    dphi = np.radians(lat2 - lat1)
    dlam = np.radians(lon2 - lon1)
    a = np.sin(dphi / 2)**2 + np.cos(p1) * np.cos(p2) * np.sin(dlam / 2)**2
    return 2 * R * np.arcsin(np.sqrt(a))


def knn_distances_km(lats, lons, k=1):
    """各点から k 番目の近傍までの距離 (km) を返す。

    高速化のため平面近似 (lat,lon -> 等距離スケール) で kd-tree を組み、
    k+1 近傍 index を取得した上で正確な haversine で距離計算する。
    """
    lat_mean = float(np.mean(lats))
    # 経度1度 ≈ cos(lat) × 緯度1度 のスケーリング
    x = (lons - np.mean(lons)) * np.cos(np.radians(lat_mean))
    y = (lats - np.mean(lats))
    tree = cKDTree(np.column_stack([x, y]))
    _, idx = tree.query(np.column_stack([x, y]), k=k + 1)
    # idx[:, 0] は自分自身、idx[:, k] が k 番目の隣
    nb = idx[:, k]
    return haversine_km(lats, lons, lats[nb], lons[nb])


print("\n=== 4. 図 3: k-NN 距離分布 ===")
KS = [1, 3, 5]
knn_results = {k: knn_distances_km(LAT, LON, k=k) for k in KS}
for k in KS:
    d = knn_results[k]
    print(f"  k={k}: 平均={d.mean():.2f}km 中央値={np.median(d):.2f}km "
          f"95%点={np.percentile(d, 95):.2f}km")

fig, axes = plt.subplots(1, 3, figsize=(13, 4), sharey=True)
COLORS_K = ["#0969da", "#1f883d", "#cf222e"]
for ax, k, c in zip(axes, KS, COLORS_K):
    d = knn_results[k]
    ax.hist(d, bins=40, color=c, edgecolor="black", linewidth=0.5, alpha=0.85)
    ax.axvline(d.mean(), color="black", linestyle="--", linewidth=1.2,
               label=f"平均 {d.mean():.2f} km")
    ax.axvline(np.median(d), color="orange", linestyle=":", linewidth=1.5,
               label=f"中央値 {np.median(d):.2f} km")
    ax.set_xlabel(f"{k} 番目の近傍までの距離 (km)")
    ax.set_title(f"k = {k}")
    ax.legend(fontsize=9)
    ax.grid(alpha=0.3)
axes[0].set_ylabel("カメラ数")
fig.suptitle("k-NN 距離分布 (haversine, n=351)", fontsize=13)
plt.tight_layout()
plt.savefig(ASSETS / "L02_knn_hist.png", dpi=140)
plt.close()


# === 5. 図 4: Voronoi 担当面積分布 ============================================
def lonlat_to_xy_km(lons, lats, lat_ref=None):
    """lon,lat を局所平面 (km) に投影 (等距離近似)。"""
    if lat_ref is None:
        lat_ref = float(np.mean(lats))
    R = 6371.0
    x = R * np.radians(lons - np.mean(lons)) * np.cos(np.radians(lat_ref))
    y = R * np.radians(lats - np.mean(lats))
    return x, y


def polygon_area(xs, ys):
    """Shoelace 公式で多角形面積 (km^2)。"""
    return 0.5 * abs(np.dot(xs, np.roll(ys, -1)) - np.dot(ys, np.roll(xs, -1)))


print("\n=== 5. 図 4: Voronoi 担当面積分布 ===")
xs, ys = lonlat_to_xy_km(LON, LAT)
pts = np.column_stack([xs, ys])
vor = Voronoi(pts)

# 凸包外の領域 (-1 を含む) を除外し、有限領域だけ面積を求める
areas = np.full(len(pts), np.nan)
finite_polys = []
for i, reg_idx in enumerate(vor.point_region):
    region = vor.regions[reg_idx]
    if not region or -1 in region:
        continue
    poly = vor.vertices[region]
    a = polygon_area(poly[:, 0], poly[:, 1])
    # 極端な外周セル (>5000 km^2) は除外 (Voronoi の発散セル対策)
    if 0 < a < 5000:
        areas[i] = a
        finite_polys.append((i, poly, a))

n_valid = int(np.sum(~np.isnan(areas)))
print(f"  有限 Voronoi セル数: {n_valid} / {len(pts)}")
print(f"  面積 平均={np.nanmean(areas):.1f}km² 中央値={np.nanmedian(areas):.1f}km² "
      f"95%点={np.nanpercentile(areas, 95):.1f}km²")

fig, axes = plt.subplots(1, 2, figsize=(12.5, 5.0))

# 左: Voronoi 図 (面積で色)
ax = axes[0]
cmap = plt.get_cmap("viridis")
log_a = np.log10(np.array([a for _, _, a in finite_polys]))
norm = plt.Normalize(vmin=np.percentile(log_a, 5),
                     vmax=np.percentile(log_a, 95))
patches = []
colors_p = []
for i, poly, a in finite_polys:
    patches.append(MplPolygon(poly, closed=True))
    colors_p.append(cmap(norm(np.log10(a))))
pc = PatchCollection(patches, facecolor=colors_p, edgecolor="white",
                     linewidth=0.3, alpha=0.85)
ax.add_collection(pc)
ax.scatter(xs, ys, s=4, c="black", alpha=0.6)
ax.set_xlabel("東西 (km, 重心起点)")
ax.set_ylabel("南北 (km, 重心起点)")
ax.set_title("Voronoi 図 (色 = 担当面積 log10 km²)")
ax.set_aspect("equal")
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([])
cb = plt.colorbar(sm, ax=ax, shrink=0.85, label="log10 面積 (km²)")
ax.set_xlim(xs.min() - 5, xs.max() + 5)
ax.set_ylim(ys.min() - 5, ys.max() + 5)
ax.grid(alpha=0.2)

# 右: 面積ヒストグラム (両対数気味)
ax = axes[1]
valid_a = areas[~np.isnan(areas)]
bins = np.logspace(np.log10(max(valid_a.min(), 0.05)),
                   np.log10(valid_a.max()), 35)
ax.hist(valid_a, bins=bins, color="#8250df", edgecolor="black",
        linewidth=0.4, alpha=0.85)
ax.set_xscale("log")
ax.axvline(np.nanmedian(valid_a), color="orange", linestyle=":",
           linewidth=1.5, label=f"中央値 {np.nanmedian(valid_a):.1f} km²")
ax.axvline(np.nanmean(valid_a), color="black", linestyle="--",
           linewidth=1.2, label=f"平均 {np.nanmean(valid_a):.1f} km²")
ax.set_xlabel("Voronoi 担当面積 (km², 対数)")
ax.set_ylabel("カメラ数")
ax.set_title(f"担当面積分布 (n={n_valid})")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)

plt.tight_layout()
plt.savefig(ASSETS / "L02_voronoi.png", dpi=140)
plt.close()

# === 6. 図 5: ランダム (Poisson) 配置との比較 ================================
print("\n=== 6. 図 5: Poisson 配置との NN 距離比較 ===")
rng = np.random.default_rng(42)

# 配置領域: 実データの bounding box
LAT_MIN, LAT_MAX = LAT.min(), LAT.max()
LON_MIN, LON_MAX = LON.min(), LON.max()

REPS = 50
poisson_nn_all = []
for _ in range(REPS):
    rlat = rng.uniform(LAT_MIN, LAT_MAX, len(LAT))
    rlon = rng.uniform(LON_MIN, LON_MAX, len(LON))
    poisson_nn_all.append(knn_distances_km(rlat, rlon, k=1))
poisson_nn = np.concatenate(poisson_nn_all)
poisson_nn_one = poisson_nn_all[0]

real_nn = knn_results[1]

ks_stat, ks_p = ks_2samp(real_nn, poisson_nn)
print(f"  実データ NN 平均 = {real_nn.mean():.3f} km")
print(f"  Poisson  NN 平均 = {poisson_nn.mean():.3f} km")
print(f"  KS 統計量 D = {ks_stat:.3f}, p = {ks_p:.2e}")

# Clark-Evans R = mean(NN_obs) / mean(NN_expected)
# CSR 期待値 = 1 / (2 * sqrt(密度)). 密度はラフに矩形領域で近似.
xs_full, ys_full = lonlat_to_xy_km(np.array([LON_MIN, LON_MAX]),
                                   np.array([LAT_MIN, LAT_MAX]),
                                   lat_ref=float(np.mean(LAT)))
area_km2 = abs((xs_full[1] - xs_full[0]) * (ys_full[1] - ys_full[0]))
density = len(LAT) / area_km2
expected_nn = 1.0 / (2.0 * np.sqrt(density))
clark_evans_R = real_nn.mean() / expected_nn
print(f"  領域面積 ≈ {area_km2:.0f} km², 密度 = {density:.4f} 点/km²")
print(f"  期待 NN 平均 (CSR) = {expected_nn:.3f} km")
print(f"  Clark-Evans R = {clark_evans_R:.3f}  "
      f"(R<1 ⇒ クラスタ的, R≈1 ⇒ ランダム, R>1 ⇒ 規則的)")

fig, axes = plt.subplots(1, 2, figsize=(12, 4.6))
ax = axes[0]
bins = np.linspace(0, max(real_nn.max(), poisson_nn_one.max()) * 1.05, 40)
ax.hist(real_nn, bins=bins, color="#cf222e", alpha=0.65,
        edgecolor="black", linewidth=0.4, label=f"実データ (n={len(real_nn)})")
ax.hist(poisson_nn_one, bins=bins, color="#0969da", alpha=0.55,
        edgecolor="black", linewidth=0.4,
        label=f"Poisson 1 サンプル (n={len(poisson_nn_one)})")
ax.axvline(real_nn.mean(), color="#cf222e", linestyle="--", linewidth=1.5)
ax.axvline(poisson_nn_one.mean(), color="#0969da", linestyle="--", linewidth=1.5)
ax.set_xlabel("最近傍距離 (km)")
ax.set_ylabel("カメラ数")
ax.set_title("NN 距離: 実データ vs Poisson")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)

# 右: 累積分布 (KS のためのCDF)
ax = axes[1]
xr = np.sort(real_nn)
yr = np.arange(1, len(xr) + 1) / len(xr)
xp = np.sort(poisson_nn)
yp = np.arange(1, len(xp) + 1) / len(xp)
ax.plot(xr, yr, color="#cf222e", linewidth=2, label="実データ CDF")
ax.plot(xp, yp, color="#0969da", linewidth=2, alpha=0.7,
        label=f"Poisson CDF ({REPS}回プール)")
ax.set_xlabel("最近傍距離 (km)")
ax.set_ylabel("累積確率")
ax.set_title(f"CDF 比較  KS D = {ks_stat:.3f}, p = {ks_p:.1e}\n"
             f"Clark-Evans R = {clark_evans_R:.3f} "
             f"({'クラスタ的' if clark_evans_R < 1 else ('規則的' if clark_evans_R > 1 else 'ランダム')})",
             fontsize=11)
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(ASSETS / "L02_poisson_compare.png", dpi=140)
plt.close()

# === 7. 図 6: 管理区分別 NN 距離 (空間集中度比較) =============================
print("\n=== 7. 図 6: 管理区分別 NN 距離 ===")
admin_results = {}
for cat in ["国", "県", "市", "その他"]:
    sub = df[df["admin_cat"] == cat]
    if len(sub) < 5:
        continue
    nn_self = knn_distances_km(sub["緯度"].to_numpy(),
                               sub["経度"].to_numpy(), k=1)
    admin_results[cat] = {
        "n": len(sub),
        "nn": nn_self,
        "mean": float(nn_self.mean()),
        "median": float(np.median(nn_self)),
    }
    print(f"  {cat}: n={len(sub)} NN平均={nn_self.mean():.2f}km "
          f"中央値={np.median(nn_self):.2f}km")

fig, ax = plt.subplots(figsize=(8.5, 4.6))
labels = list(admin_results.keys())
data = [admin_results[c]["nn"] for c in labels]
bp = ax.boxplot(data, tick_labels=labels, patch_artist=True,
                showfliers=True, widths=0.55,
                medianprops={"color": "black", "linewidth": 2})
PALETTE = ["#0969da", "#cf222e", "#1f883d", "#fb8500"]
for patch, c in zip(bp["boxes"], PALETTE):
    patch.set_facecolor(c)
    patch.set_alpha(0.65)
for i, c in enumerate(labels):
    n = admin_results[c]["n"]
    mn = admin_results[c]["mean"]
    ax.text(i + 1, ax.get_ylim()[1] * 0.95,
            f"n={n}\n平均{mn:.2f}km", ha="center", fontsize=9)
ax.set_ylabel("最近傍距離 (km, 同区分内)")
ax.set_title("管理者区分別 同内 NN 距離分布\n(短いほど同種カメラが近接 = 集中度高い)")
ax.grid(alpha=0.3, axis="y")
plt.tight_layout()
plt.savefig(ASSETS / "L02_admin_nn.png", dpi=140)
plt.close()

# === 8. サマリ表 =============================================================
print("\n=== 8. サマリ表 ===")
summary_rows = [
    {"指標": "総カメラ数 (有効緯度経度)", "値": f"{len(df)} 台"},
    {"指標": "領域面積 (bounding box)", "値": f"{area_km2:.0f} km²"},
    {"指標": "点密度", "値": f"{density:.4f} 点/km²"},
    {"指標": "実 NN 平均距離", "値": f"{real_nn.mean():.2f} km"},
    {"指標": "CSR 期待 NN (1/(2√λ))", "値": f"{expected_nn:.2f} km"},
    {"指標": "Clark-Evans R", "値": f"{clark_evans_R:.3f}"},
    {"指標": "KS 検定 (実 vs Poisson)",
     "値": f"D={ks_stat:.3f}, p={ks_p:.2e}"},
    {"指標": "Voronoi 担当面積 中央値",
     "値": f"{np.nanmedian(valid_a):.1f} km²"},
    {"指標": "Voronoi 担当面積 平均",
     "値": f"{np.nanmean(valid_a):.1f} km²"},
]
summary_df = pd.DataFrame(summary_rows)
summary_html = summary_df.to_html(index=False)
print(summary_df.to_string(index=False))

admin_rows = []
for c, R in admin_results.items():
    admin_rows.append({
        "管理者区分": c, "n": R["n"],
        "NN 平均 (km)": round(R["mean"], 3),
        "NN 中央値 (km)": round(R["median"], 3),
    })
admin_df = pd.DataFrame(admin_rows)
admin_html = admin_df.to_html(index=False)

# === 9. HTML レンダリング ====================================================
print("\n=== 9. HTML レンダリング ===")
sections = [
    ("学習目標", """
<ul>
<li>位置情報付きCSVを folium で <b>多層 (MarkerCluster + HeatMap) マップ</b> として描き、用途別に切り替えて見せられる</li>
<li><b>haversine 距離</b>で球面上の距離を扱い、<b>k-NN (k=1,3,5) 距離分布</b>を計算できる</li>
<li><code>scipy.spatial.Voronoi</code> で各点の <b>担当面積分布</b>を出し、点配置の偏りを定量化できる</li>
<li><b>同じ領域に Poisson 配置を生成</b>し、最近傍距離分布を実データと比較 (KS 検定 / Clark-Evans R) して
「ランダムか／クラスタ的か」を <b>統計的に判定</b>できる</li>
<li>管理者区分 (国/県/市) で <b>空間集中度</b>がどう違うかを箱ひげ図で比較できる</li>
</ul>"""),

    ("使用データ", """
<ul class="kv">
<li><b>名称</b> 県内のカメラ情報 (道路・河川・ため池 ほか)</li>
<li><b>出典</b> <a href="https://hiroshima-dobox.jp/datasets/1279" target="_blank">DoBoX dataset 1279</a></li>
<li><b>件数</b> 351 台 (本実装で全行使用)</li>
<li><b>列</b> No., カメラ名, 住所, 路河川名等, <b>緯度, 経度</b>, 公開URL, 所管, 管理区分</li>
<li><b>空間範囲</b> 北緯 34.09 〜 35.08, 東経 132.09 〜 133.41 (広島県全域)</li>
</ul>"""),

    ("データ取得手順", data_recipe([
        {"label": "県内カメラ情報", "dataset_id": 1279,
         "file": "data/camera_list.csv",
         "format": "CSV (UTF-8, ヘッダ1行, 351 行)",
         "size": "約 70 KB"},
    ])),

    ("方法", """
<ol>
<li><b>多層 folium マップ</b>: <code>MarkerCluster</code> (区分別色) と <code>HeatMap</code> (密度) を 2 レイヤとして重ね、<code>LayerControl</code> で切替可能にする → <code>L02_map.html</code></li>
<li><b>k-NN 距離</b>: <code>cKDTree</code> で k+1 近傍 index を高速取得 → 正確な haversine で距離 (km) を確定。k=1, 3, 5 でヒストグラム化</li>
<li><b>Voronoi 担当面積</b>: 緯度経度 → 局所平面 (km, 等距離近似) に投影 → <code>scipy.spatial.Voronoi</code> → 有限セルだけ Shoelace で面積を出し対数ヒストグラム</li>
<li><b>Poisson 比較</b>: 同 bounding box 内に同数 (351点) の一様乱数を 50 回生成 → 実 NN 分布と KS 検定 + Clark-Evans R</li>
<li><b>管理者別比較</b>: 所管文字列を 国/県/市/その他 に粗分類 → 同区分内の NN 距離分布を箱ひげで比較。短い ⇒ 集中、長い ⇒ 散在</li>
</ol>"""),

    ("コード解説", code('''
import numpy as np, pandas as pd, folium
from folium.plugins import MarkerCluster, HeatMap
from scipy.spatial import Voronoi, cKDTree
from scipy.stats import ks_2samp
from _common import ensure_dataset

# === (0) データ自動取得 ===
ensure_dataset("data/camera_list.csv", dataset_id=1279, label="県内カメラ #1279")
df = pd.read_csv("data/camera_list.csv").dropna(subset=["緯度", "経度"])
LAT, LON = df["緯度"].to_numpy(), df["経度"].to_numpy()

# === (1) folium 多層マップ: MarkerCluster + HeatMap ===
m = folium.Map(location=[LAT.mean(), LON.mean()], zoom_start=9)
cluster_layer = folium.FeatureGroup(name="MarkerCluster", show=True)
mc = MarkerCluster().add_to(cluster_layer)
COLOR = {"道路":"blue","河川":"green","ため池":"red","海岸":"orange","港湾":"purple","その他":"gray"}
for _, r in df.iterrows():
    folium.CircleMarker([r["緯度"], r["経度"]], radius=4,
        color=COLOR.get(r["管理区分"], "gray"), fill=True).add_to(mc)
cluster_layer.add_to(m)
heat_layer = folium.FeatureGroup(name="HeatMap", show=False)
HeatMap(list(zip(LAT, LON)), radius=18).add_to(heat_layer)
heat_layer.add_to(m)
folium.LayerControl().add_to(m)
m.save("assets/L02_map.html")

# === (2) k-NN 距離 (haversine) ===
def haversine_km(lat1, lon1, lat2, lon2):
    R = 6371.0
    p1, p2 = np.radians(lat1), np.radians(lat2)
    dphi = np.radians(lat2 - lat1); dlam = np.radians(lon2 - lon1)
    a = np.sin(dphi/2)**2 + np.cos(p1)*np.cos(p2)*np.sin(dlam/2)**2
    return 2*R*np.arcsin(np.sqrt(a))

def knn_km(lats, lons, k=1):
    # 緯度経度 → 局所平面で kd-tree を組み、index 取得後に haversine で測る
    lat0 = lats.mean()
    x = (lons - lons.mean()) * np.cos(np.radians(lat0))
    y = (lats - lats.mean())
    tree = cKDTree(np.column_stack([x, y]))
    _, idx = tree.query(np.column_stack([x, y]), k=k+1)
    nb = idx[:, k]   # 自分自身は idx[:, 0]
    return haversine_km(lats, lons, lats[nb], lons[nb])

knn1 = knn_km(LAT, LON, k=1)   # 各カメラから一番近い隣カメラまで (km)

# === (3) Voronoi 担当面積 ===
R = 6371.0
lat0 = LAT.mean()
x = R * np.radians(LON - LON.mean()) * np.cos(np.radians(lat0))
y = R * np.radians(LAT - LAT.mean())
vor = Voronoi(np.column_stack([x, y]))
areas = []
for i, ri in enumerate(vor.point_region):
    region = vor.regions[ri]
    if not region or -1 in region:   # 凸包外は除外
        continue
    poly = vor.vertices[region]
    a = 0.5*abs(np.dot(poly[:,0], np.roll(poly[:,1], -1))
              - np.dot(poly[:,1], np.roll(poly[:,0], -1)))
    if 0 < a < 5000:
        areas.append(a)

# === (4) Poisson 配置との比較 (KS + Clark-Evans R) ===
rng = np.random.default_rng(42)
poisson_nn = np.concatenate([
    knn_km(rng.uniform(LAT.min(), LAT.max(), len(LAT)),
           rng.uniform(LON.min(), LON.max(), len(LON)), k=1)
    for _ in range(50)])
ks_stat, ks_p = ks_2samp(knn1, poisson_nn)
density = len(LAT) / 12000.0   # 矩形領域面積 (km^2)
expected_nn = 1 / (2 * np.sqrt(density))
clark_evans = knn1.mean() / expected_nn
print(f"D={ks_stat:.3f} p={ks_p:.2e}  CE-R={clark_evans:.3f}")
# R<1 ⇒ クラスタ的, R≈1 ⇒ ランダム, R>1 ⇒ 規則的

# === (5) 管理者区分別 NN ===
def admin_cat(s):
    s = str(s)
    if "国土交通省" in s: return "国"
    if s.endswith("県"):  return "県"
    if s.endswith("市"):  return "市"
    return "その他"
df["admin_cat"] = df["所管"].map(admin_cat)
for cat in ["国","県","市","その他"]:
    sub = df[df["admin_cat"]==cat]
    nn = knn_km(sub["緯度"].values, sub["経度"].values, k=1)
    print(cat, "n=", len(sub), "NN平均", nn.mean())
''')),

    ("結果",
     "<h3>1. 多層インタラクティブ地図 (MarkerCluster + HeatMap)</h3>"
     '<iframe class="map" src="assets/L02_map.html"></iframe>'
     '<p class="tnote">※ 右上のレイヤ切替で <b>MarkerCluster</b> (区分別色付き、ズーム連動でクラスタ展開) と '
     '<b>HeatMap</b> (密度) を切り替え可能。各点クリックでカメラ名・所管・ライブ映像URLを表示。</p>' +
     figure("assets/L02_overview.png",
            "左: 管理区分別カメラ数 (道路131・河川120・ため池70・海岸18・港湾6・その他6) ／ 右: 経度×緯度散布。沿岸 (南部) と内陸の県中央部に集中") +
     figure("assets/L02_knn_hist.png",
            "k-NN 距離分布 (k=1,3,5)。中央値が k に対し概ね 1.5〜2 倍程度に伸びるのは局所的にカメラが密集している証拠 (一様分布なら √k スケール)") +
     figure("assets/L02_voronoi.png",
            "左: Voronoi 図 (色 = log10 担当面積)。沿岸の道路・河川カメラは小さなセル、北部山地は大きなセル ／ 右: 担当面積の対数ヒストグラム。3 桁にまたがり強い右裾分布") +
     figure("assets/L02_poisson_compare.png",
            "左: 実データ NN 距離 vs 同領域 Poisson 配置 NN 距離 ／ 右: CDF 比較と KS 検定。実データが Poisson より小さい距離側に偏る = クラスタ的。Clark-Evans R が判定指標") +
     figure("assets/L02_admin_nn.png",
            "管理者区分別 同内 NN 距離分布。国 (国交省) は限られた幹線沿いに集中するため箱が低め、県・市は広域分散") +
     "<h3>主要指標サマリ</h3>" + summary_html +
     "<h3>管理者区分別 同内 NN 距離</h3>" + admin_html),

    ("考察", """
<ul>
<li><b>「ランダム」より明らかにクラスタ的</b>: Clark-Evans R が 1 を大きく下回り、KS 検定でも実データと Poisson 配置の差が極めて有意 (p ≪ 0.001)。
カメラは <b>監視対象 (道路・川筋・ため池)</b>に沿って配置されるため、空間的に均一にはならない — これは <b>定量的に</b>裏取りできた。</li>
<li><b>Voronoi 面積の対数ヒストグラム</b>が右裾を引くのは、北部山地に「面積数百 km² の担当セル」が少数できる一方、
広島市〜呉のベルト地帯に「面積 1 km² 未満のセル」が多数できているため。<b>「平均1台あたり何 km²」は誤解を招く</b>:
中央値とパーセンタイルで語るのが正しい。</li>
<li><b>管理者区分による空間構造の違い</b>: 国 (国交省) のカメラは <b>主要国道・1級河川</b>沿いの線状配置で同内 NN が短い。
県・市は <b>広域に分散</b>するため同内 NN が長い。「データの集中度はそのまま行政の役割分担を映す」。</li>
<li><b>k-NN の k に対する伸び方</b>は <b>密度が局所的に違う</b>こと (= 不均質な点過程) を示唆。
均一 Poisson なら NN_k の中央値は √k に比例して伸びるが、実データはそれより緩やかに伸びる箇所もある — 高密度クラスタの内側に居る点は k=5 でも近傍が尽きないためで、
<b>「クラスタの大きさ」を測る教材</b> (Ripley K) への自然な接続点になる。</li>
<li><b>folium 2 レイヤの教育的意義</b>: <b>MarkerCluster は「個々を見る」</b>、<b>HeatMap は「全体傾向を見る」</b>。
切り替えで「ミクロ ⇄ マクロ」を体験できる。可視化は<b>1枚の地図に詰め込まない</b>のが鉄則。</li>
</ul>"""),

    ("発展課題", """
<ol>
<li><b>Ripley K 関数</b>を実装して、距離 r を変えながら「どのスケールでクラスタ的か」を出す。今回の Clark-Evans R は r=NN距離スケールでの 1 点指標にすぎない。</li>
<li><b>道路ネットワーク距離</b>での NN 計算: ユークリッド距離ではなく <b>OSMnx</b> で取った道路グラフ上の距離で再計算すると、海越し・山越しの不自然な近接が消える。</li>
<li><b>属性条件付き</b>: 「ため池カメラから一番近い河川カメラまで」のような <b>異種 NN</b>。L09 (避難所×カメラ) と同じ枠組みで多変量化。</li>
<li><b>動的 Poisson</b>: 矩形ではなく広島県境界 GeoJSON で領域を取り直し、より厳密な CSR を定義。</li>
<li><b>稼働率と組合せ</b>: <code>requests.head</code> で各 <code>公開URL</code> を叩き、応答 200 のカメラだけで NN を再計算。「公開されているカメラのカバレッジ」が見える。</li>
<li><b>市区町村別 1人あたりカメラ数</b>: 住所列から市町村抽出 → 国勢調査人口で割って <b>カメラ密度の格差</b>を地図化 (L09 接続)。</li>
</ol>"""),
]

html = render_lesson(
    num=2,
    title="県内カメラ351点 — folium 多層マップと点配置の空間統計",
    tags=["導入", "基礎", "可視化", "地理", "空間統計", "v2-rewrite"],
    time="90分",
    level="基礎〜応用基礎",
    data_label='<a href="https://hiroshima-dobox.jp/datasets/1279" target="_blank">県内のカメラ情報 (DoBoX #1279)</a>',
    sections=sections,
    script_filename="lessons/L02_camera_map.py",
)
out = LESSONS / "L02_camera_map.html"
out.write_text(html, encoding="utf-8")
print(f"\nsaved: {out}")
print(f"  PNGs: 5 + folium HTML 1 generated in {ASSETS}")