"""
案例 4：××市低空巡检遥感监测数据分析
输入：低空巡检监测数据.csv / 巡检飞行日志.csv
输出：4 张分析图 + 结论数据（打印到 stdout）
"""
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib import font_manager

BASE = "/Users/frank/WorkBuddy/2026-08-09-23-27-38/geovis-training"

# ---------- 中文字体 ----------
for fp in ["/System/Library/Fonts/PingFang.ttc",
           "/System/Library/Fonts/Hiragino Sans GB.ttc",
           "/Library/Fonts/Arial Unicode.ttf"]:
    try:
        font_manager.fontManager.addfont(fp)
        name = font_manager.FontProperties(fname=fp).get_name()
        plt.rcParams["font.family"] = name
        print(f"[字体] 使用 {name}")
        break
    except Exception:
        continue
plt.rcParams["axes.unicode_minus"] = False

df = pd.read_csv(f"{BASE}/data/低空巡检监测数据.csv")
log = pd.read_csv(f"{BASE}/data/巡检飞行日志.csv")
df["date"] = pd.to_datetime(df["date"])
log["date"] = pd.to_datetime(log["date"])

print("=" * 62)
print("一、数据概览")
print("=" * 62)
print(f"监测周期      : {df.date.min():%Y-%m-%d} ~ {df.date.max():%Y-%m-%d}（{(df.date.max()-df.date.min()).days+1} 天）")
print(f"有效巡检架次  : {df.sortie_id.nunique()}")
print(f"覆盖网格      : {df.grid_id.nunique()} 个")
print(f"识别问题总数  : {len(df)}")
print(f"停飞天数      : {int(log.grounded.sum())} 天（占比 {log.grounded.mean()*100:.1f}%）")
print(f"高置信度占比  : {(df.confidence>=0.8).mean()*100:.1f}%（confidence>=0.8）")
print(f"人工复核确认率: {df.verified.mean()*100:.1f}%")

CLR = {"违法建设":"#d94f4f","裸土未覆盖":"#d98a3c","垃圾堆放":"#4f7fd9",
       "水体异色":"#3ca39c","秸秆焚烧":"#8b5fbf"}

fig = plt.figure(figsize=(15.5, 10.5))
fig.suptitle("××市低空巡检遥感监测分析报告（2026-05-01 ~ 2026-08-08）",
             fontsize=15, y=0.975)

# ===== 图1 问题类型时序趋势（周聚合） =====
ax1 = fig.add_subplot(2, 2, 1)
wk = df.groupby([pd.Grouper(key="date", freq="W"), "problem_type"]).size().unstack(fill_value=0)
for c in wk.columns:
    ax1.plot(wk.index, wk[c], marker="o", ms=3.5, lw=1.7, label=c, color=CLR[c])
ax1.axvline(pd.Timestamp("2026-06-15"), color="#888", ls="--", lw=1.2)
ax1.text(pd.Timestamp("2026-06-16"), ax1.get_ylim()[1]*0.93, "6/15 加强执法",
         fontsize=9, color="#666")
ax1.set_title("问题类型周度趋势", fontsize=12)
ax1.set_ylabel("问题数量")
ax1.legend(fontsize=8.5, ncol=2, frameon=False)
ax1.grid(alpha=.25, ls=":")
for s in ("top","right"): ax1.spines[s].set_visible(False)

# ===== 图2 网格类型 × 问题类型热力 =====
ax2 = fig.add_subplot(2, 2, 2)
pv = df.pivot_table(index="grid_type", columns="problem_type",
                    values="sortie_id", aggfunc="count", fill_value=0)
# 归一化为每架次问题密度
sor = df.groupby("grid_type")["sortie_id"].nunique()
pvn = pv.div(sor, axis=0)
im = ax2.imshow(pvn.values, cmap="YlOrRd", aspect="auto")
ax2.set_xticks(range(len(pvn.columns)))
ax2.set_xticklabels(pvn.columns, fontsize=9, rotation=18, ha="right")
ax2.set_yticks(range(len(pvn.index)))
ax2.set_yticklabels(pvn.index, fontsize=9)
for i in range(pvn.shape[0]):
    for j in range(pvn.shape[1]):
        v = pvn.values[i, j]
        ax2.text(j, i, f"{v:.2f}", ha="center", va="center", fontsize=9,
                 color="white" if v > pvn.values.max()*0.6 else "#333")
ax2.set_title("单架次问题密度（问题数/架次）", fontsize=12)
fig.colorbar(im, ax=ax2, fraction=0.042, pad=0.03)

# ===== 图3 成像质量对识别置信度的影响（分箱） =====
ax3 = fig.add_subplot(2, 2, 3)
fly = log[log.grounded == 0].copy()
bins = [0, 6, 9, 12, 16, 40]
labels = ["<6", "6-9", "9-12", "12-16", ">16"]
df["vis_bin"] = pd.cut(df.visibility_km, bins=bins, labels=labels)
gb = df.groupby("vis_bin", observed=True).agg(
    conf=("confidence", "mean"),
    qual=("img_quality", "mean"),
    n=("confidence", "size"))
xpos = np.arange(len(gb))
b1 = ax3.bar(xpos - 0.2, gb.qual, width=0.38, label="平均成像质量",
             color="#4f7fd9", edgecolor="none")
b2 = ax3.bar(xpos + 0.2, gb.conf, width=0.38, label="平均识别置信度",
             color="#d98a3c", edgecolor="none")
for b in list(b1) + list(b2):
    ax3.text(b.get_x() + b.get_width()/2, b.get_height() + 0.012,
             f"{b.get_height():.3f}", ha="center", fontsize=8.5)
ax3.set_xticks(xpos)
ax3.set_xticklabels(labels)
ax3.set_xlabel("能见度区间 (km)")
ax3.set_ylabel("质量 / 置信度")
ax3.set_ylim(0, 1.08)
ax3.set_title("能见度对成像质量与识别置信度的影响", fontsize=12)
ax3.legend(fontsize=9, frameon=False, loc="upper center",
           bbox_to_anchor=(0.5, -0.16), ncol=2)
ax3.grid(alpha=.25, ls=":", axis="y")
for s in ("top","right"): ax3.spines[s].set_visible(False)
for i, n in enumerate(gb.n):
    ax3.text(i, 0.045, f"n={n}", ha="center", fontsize=8, color="#555")

# ===== 图4 网格问题排名 TOP12 =====
ax4 = fig.add_subplot(2, 2, 4)
gs = df.groupby(["grid_id","grid_type"]).size().reset_index(name="n")
gs = gs.sort_values("n", ascending=False).head(12).iloc[::-1]
tc = {"城中村":"#d94f4f","工业园区":"#d98a3c","城乡结合部":"#4f7fd9","水域岸线":"#3ca39c"}
ax4.barh(gs.grid_id + " · " + gs.grid_type, gs.n,
         color=[tc[t] for t in gs.grid_type], height=.66)
for y, v in enumerate(gs.n):
    ax4.text(v + max(gs.n)*0.012, y, str(v), va="center", fontsize=9)
ax4.set_title("重点网格问题数排名 TOP12", fontsize=12)
ax4.set_xlabel("问题数量")
ax4.grid(alpha=.25, ls=":", axis="x")
for s in ("top","right"): ax4.spines[s].set_visible(False)

plt.tight_layout(rect=[0, 0, 1, 0.955])
out = f"{BASE}/assets/case4-分析图.png"
plt.savefig(out, dpi=155, bbox_inches="tight", facecolor="white")
print(f"\n[输出] {out}")

# ---------- 结论数据 ----------
print("\n" + "=" * 62)
print("二、关键分析结论")
print("=" * 62)

pre = df[df.date < "2026-06-15"]
post = df[df.date >= "2026-06-15"]
pre_d = (pd.Timestamp("2026-06-15") - df.date.min()).days
post_d = (df.date.max() - pd.Timestamp("2026-06-15")).days + 1

print("\n[1] 执法干预效果（6/15 前后日均问题数对比）")
for pt in ["违法建设", "垃圾堆放", "裸土未覆盖", "水体异色"]:
    a = len(pre[pre.problem_type == pt]) / pre_d
    b = len(post[post.problem_type == pt]) / post_d
    chg = (b - a) / a * 100 if a else 0
    flag = "↓ 显著下降" if chg < -20 else ("↑ 上升" if chg > 10 else "— 基本持平")
    print(f"    {pt:6s}: {a:5.2f} → {b:5.2f} 件/天  ({chg:+6.1f}%)  {flag}")

print("\n[2] 秸秆焚烧季节性异常")
jun = df[(df.problem_type=="秸秆焚烧") & (df.date>="2026-06-01") & (df.date<="2026-06-12")]
oth = df[(df.problem_type=="秸秆焚烧") & ~((df.date>="2026-06-01") & (df.date<="2026-06-12"))]
r1 = len(jun)/12
r2 = len(oth)/((df.date.max()-df.date.min()).days+1-12)
print(f"    6/1-6/12 麦收期: {r1:.2f} 件/天")
print(f"    其余时段        : {r2:.2f} 件/天")
print(f"    倍数            : {r1/r2:.1f}×  → 明确季节性，应提前布防")

print("\n[3] 气象对巡检的量化影响")
print(f"    停飞天数        : {int(log.grounded.sum())} 天（风速>10m/s 或 能见度<3km 或 降水>8mm）")
lo = fly[fly.visibility_km < 8]; hi = fly[fly.visibility_km >= 8]
print(f"    能见度<8km 日均问题数 : {lo.problems.mean():.1f}")
print(f"    能见度>=8km 日均问题数: {hi.problems.mean():.1f}")
print(f"    差异: {(hi.problems.mean()-lo.problems.mean())/hi.problems.mean()*100:.1f}% 的识别能力损失")
q_lo = df[df.img_quality < 0.8].confidence.mean()
q_hi = df[df.img_quality >= 0.9].confidence.mean()
print(f"    成像质量<0.8 时平均置信度: {q_lo:.3f}")
print(f"    成像质量>=0.9 时平均置信度: {q_hi:.3f}")

print("\n[4] 网格风险画像")
for gt in ["城中村","工业园区","城乡结合部","水域岸线"]:
    sub = df[df.grid_type == gt]
    top = sub.problem_type.value_counts().idxmax()
    dens = len(sub) / sub.sortie_id.nunique()
    print(f"    {gt:6s}: 首要问题「{top}」  单架次问题密度 {dens:.2f}")

print("\n[5] 复核与置信度")
hc = df[df.confidence >= 0.8]
lc = df[df.confidence < 0.8]
print(f"    高置信度(>=0.8) 复核确认率: {hc.verified.mean()*100:.1f}%（{len(hc)} 件）")
print(f"    低置信度(<0.8)  复核确认率: {lc.verified.mean()*100:.1f}%（{len(lc)} 件）")
print(f"    → 建议：置信度<0.8 的告警必须人工复核，否则误报将进入执法流程")

print("\n" + "=" * 62)