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

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())