# -*- coding: utf-8 -*-
"""
===========================================================================
模組五｜台灣夏季備轉容量率惡化趨勢
當氣候變遷撞上電力系統：剪刀效應的在地版
===========================================================================

【教學目標｜ILO】
  1. 學習者能自動從台電與政府資料開放平臺下載最新電力供需 CSV。
  2. 學習者能處理 Big5 / UTF-8 雙編碼情境與民國年 / 西元年轉換。
  3. 學習者能計算「警戒日」（yellow/orange/red alert）並以熱地圖視覺化。
  4. 學習者能以 OLS 迴歸量化「尖峰負載攀升 vs 備用容量率下降」的
     在地剪刀效應，並連結至模組一（櫻花）所呈現的全球剪刀效應。

【資料來源】
  本模組直接從台電官方 API 抓取資料，每次執行都會得到最新數字。

  (1) 歷年尖峰負載及備用容量率（自民國 71 年＝西元 1982 年起）
      資料集代號：政府資料開放平臺 #8307
      URL: https://service.taipower.com.tw/data/opendata/
            apply/file/d003002/001.csv
      欄位：年度、尖峰負載(MW)、備用容量率(%)
      ※ 民國 111 年（2022）起改以「夜間備用容量率」呈現

  (2) 近三年每日尖峰備轉容量率
      資料集代號：政府資料開放平臺 #24945
      URL: https://service.taipower.com.tw/data/opendata/
            apply/file/d006004/001.csv
      欄位：日期、備轉容量(萬瓩)、備轉容量率(%)

【重要名詞釐清】
  - 備用容量率（Reserve Margin）：年度指標，＝（裝置容量 − 尖峰負載）/ 尖峰負載
  - 備轉容量率（Operating Reserve Margin）：每日指標，＝（運轉淨尖峰能力 −
    尖峰負載）/ 尖峰負載。本模組主要使用此指標。

【台電燈號制度】
  - 綠燈 Green   : 備轉容量率 ≥ 10%（供電充裕）
  - 黃燈 Yellow  : 6% ≤ 備轉容量率 < 10%（供電吃緊）
  - 橘燈 Orange  : 備轉容量率 < 6%（供電警戒）
  - 紅燈 Red     : 備轉容量 < 50 萬瓩（限電警戒）
  - 黑燈 Black   : 實施限電

【運算環境】
  Python 3.10+（Colab / 本地皆可）
  相依套件：pandas, numpy, matplotlib, statsmodels, requests

【Colab 首次執行】
  !apt-get install -y fonts-noto-cjk -q

===========================================================================
"""

# ============================================================================
# 第 0 節｜環境設定
# ============================================================================
import io
import warnings
import urllib.request
import urllib.error
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib import rcParams
from matplotlib.colors import ListedColormap, BoundaryNorm
import statsmodels.api as sm

warnings.filterwarnings("ignore")

rcParams["font.sans-serif"] = [
    "Noto Sans CJK TC", "Microsoft JhengHei",
    "PingFang TC", "Arial Unicode MS", "DejaVu Sans",
]
rcParams["axes.unicode_minus"] = False
rcParams["figure.dpi"] = 110


# ============================================================================
# 第 1 節｜穩健下載器（處理 Big5 / UTF-8 / User-Agent 問題）
# ============================================================================
def download_taipower_csv(url: str,
                          encodings: tuple = ("utf-8-sig", "utf-8",
                                              "big5", "cp950")) -> pd.DataFrame:
    """
    台電 CSV 下載器，依序嘗試多種編碼。

    台電檔案的編碼曾在 Big5 與 UTF-8 之間變動過，為使程式長期可用，
    本函式以「依序嘗試，成功即返回」的方式處理編碼不確定性。
    若遇 403 Forbidden，會加上 User-Agent 重試。
    """
    req = urllib.request.Request(
        url, headers={"User-Agent":
                      "Mozilla/5.0 (educational research; Python/pandas)"})
    try:
        with urllib.request.urlopen(req, timeout=30) as resp:
            raw = resp.read()
    except urllib.error.HTTPError as e:
        raise RuntimeError(f"無法下載：HTTP {e.code} {e.reason}（{url}）")
    except urllib.error.URLError as e:
        raise RuntimeError(f"無法連線：{e.reason}（{url}）")

    last_err = None
    for enc in encodings:
        try:
            text = raw.decode(enc)
            df = pd.read_csv(io.StringIO(text))
            return df
        except UnicodeDecodeError as e:
            last_err = e
            continue
        except pd.errors.ParserError as e:
            last_err = e
            continue
    raise RuntimeError(f"所有編碼皆解碼失敗（嘗試：{encodings}）；"
                       f"最後錯誤：{last_err}")


# ============================================================================
# 第 2 節｜讀取歷年尖峰負載及備用容量率（Dataset 8307）
# ============================================================================
def load_annual_reserve() -> pd.DataFrame:
    """
    讀取 1982 年起之年度尖峰負載與備用容量率。

    民國年轉西元：西元年 = 民國年 + 1911
    例：民國 71 年 = 1982 年；民國 114 年 = 2025 年

    欄位名稱可能隨年份變動，本函式以關鍵字模糊匹配找出正確欄位。
    """
    url = ("https://service.taipower.com.tw/data/opendata/"
           "apply/file/d003002/001.csv")
    df = download_taipower_csv(url)

    # 欄位模糊匹配
    col_year = _find_col(df, ["年度", "年", "year"])
    col_load = _find_col(df, ["尖峰負載", "負載", "peak"])
    col_margin = _find_col(df, ["備用容量率", "備用", "margin"])

    out = df[[col_year, col_load, col_margin]].copy()
    out.columns = ["ROC_year", "peak_load_MW", "reserve_margin_pct"]

    # 清理
    out["ROC_year"] = pd.to_numeric(out["ROC_year"], errors="coerce")
    out["peak_load_MW"] = pd.to_numeric(out["peak_load_MW"], errors="coerce")
    out["reserve_margin_pct"] = pd.to_numeric(
        out["reserve_margin_pct"], errors="coerce")
    out = out.dropna(subset=["ROC_year"]).copy()

    # 民國年轉西元
    out["AD_year"] = out["ROC_year"].astype(int) + 1911
    out = out.sort_values("AD_year").reset_index(drop=True)

    print(f"[年度資料] 共讀入 {len(out)} 年；"
          f"範圍 {out['AD_year'].min()}–{out['AD_year'].max()}")
    return out


def _find_col(df: pd.DataFrame, keywords: list) -> str:
    """在 DataFrame 欄位中找含任一關鍵字的欄；找不到則 raise。"""
    for col in df.columns:
        for kw in keywords:
            if kw in str(col):
                return col
    raise KeyError(f"找不到欄位；關鍵字 {keywords}；"
                   f"現有欄位 {list(df.columns)}")


# ============================================================================
# 第 3 節｜讀取近三年每日備轉容量率（Dataset 24945）
# ============================================================================
def load_daily_operating_reserve() -> pd.DataFrame:
    """
    讀取近三年每日備轉容量率。

    日期欄位格式可能為 'YYYY/MM/DD' 或 'YYYY-MM-DD'，pandas 可自動解析。
    """
    url = ("https://service.taipower.com.tw/data/opendata/"
           "apply/file/d006004/001.csv")
    df = download_taipower_csv(url)

    col_date = _find_col(df, ["日期", "date"])
    col_cap = _find_col(df, ["備轉容量(萬", "備轉容量 (萬",
                             "operating reserve"])
    col_rate = _find_col(df, ["備轉容量率", "%", "rate"])

    out = df[[col_date, col_cap, col_rate]].copy()
    out.columns = ["date", "operating_reserve_万瓩", "operating_reserve_pct"]

    out["date"] = pd.to_datetime(out["date"], errors="coerce")
    out["operating_reserve_万瓩"] = pd.to_numeric(
        out["operating_reserve_万瓩"], errors="coerce")
    out["operating_reserve_pct"] = pd.to_numeric(
        out["operating_reserve_pct"], errors="coerce")
    out = out.dropna(subset=["date", "operating_reserve_pct"]).copy()

    out["year"] = out["date"].dt.year
    out["month"] = out["date"].dt.month
    out["day_of_year"] = out["date"].dt.dayofyear

    out = out.sort_values("date").reset_index(drop=True)

    print(f"[日資料] 共讀入 {len(out):,} 日；"
          f"範圍 {out['date'].min().strftime('%Y-%m-%d')} → "
          f"{out['date'].max().strftime('%Y-%m-%d')}")
    return out


# ============================================================================
# 第 4 節｜長期剪刀效應分析：負載上升 × 容量率下降
# ============================================================================
def long_term_scissors_regression(annual: pd.DataFrame) -> dict:
    """
    兩條 OLS 迴歸並列：
      (1) 尖峰負載 ~ 年份：上升斜率
      (2) 備用容量率 ~ 年份：下降斜率
    兩者相反方向延伸，即為「在地剪刀效應」。
    """
    print("\n─── 長期趨勢 OLS 迴歸（1982 年起）───")

    # 負載
    X = sm.add_constant(annual["AD_year"].values)
    y = annual["peak_load_MW"].values
    m_load = sm.OLS(y, X).fit()
    print(f"尖峰負載：β = {m_load.params[1]:+.1f} MW/年，"
          f"每十年上升 {m_load.params[1]*10:.0f} MW ({m_load.rsquared*100:.1f}% 變異)")

    # 容量率
    y2 = annual["reserve_margin_pct"].values
    m_margin = sm.OLS(y2, X).fit()
    print(f"備用容量率：β = {m_margin.params[1]:+.3f} %/年，"
          f"每十年下降 {-m_margin.params[1]*10:.2f} 百分點 "
          f"({m_margin.rsquared*100:.1f}% 變異)")

    # 近年最低
    last_5 = annual.tail(5)
    print(f"\n最近 5 年備用容量率均值：{last_5['reserve_margin_pct'].mean():.1f} %")
    print(f"歷史最低年份：{int(annual.loc[annual['reserve_margin_pct'].idxmin(), 'AD_year'])}"
          f"，為 {annual['reserve_margin_pct'].min():.1f} %")

    return {"load_model": m_load, "margin_model": m_margin}


# ============================================================================
# 第 5 節｜夏季（6–9 月）警戒日分析
# ============================================================================
def summer_alert_days(daily: pd.DataFrame) -> pd.DataFrame:
    """
    計算每年夏季（6–9 月）的燈號分布。

    台電燈號定義（簡化版，以備轉容量率為準）：
      綠燈 green  : ≥ 10%
      黃燈 yellow : 6% – 10%
      橘燈 orange : 3% – 6%
      紅燈 red    : < 3%
    """
    summer = daily[daily["month"].between(6, 9)].copy()

    def classify(pct):
        if pct >= 10:
            return "綠燈"
        elif pct >= 6:
            return "黃燈"
        elif pct >= 3:
            return "橘燈"
        else:
            return "紅燈"

    summer["alert"] = summer["operating_reserve_pct"].apply(classify)

    # 依年份交叉統計
    cross = (summer.groupby(["year", "alert"])
             .size().unstack(fill_value=0))
    # 確保四種燈號都存在
    for c in ["綠燈", "黃燈", "橘燈", "紅燈"]:
        if c not in cross.columns:
            cross[c] = 0
    cross = cross[["綠燈", "黃燈", "橘燈", "紅燈"]]
    cross["夏季總天數"] = cross.sum(axis=1)
    cross["吃緊比例%"] = ((cross["黃燈"] + cross["橘燈"] + cross["紅燈"])
                      / cross["夏季總天數"] * 100).round(1)

    print("\n─── 夏季（6–9 月）警戒日分布 ───")
    print(cross.to_string())
    return cross


# ============================================================================
# 第 6 節｜視覺化（一）：長期剪刀效應主圖
# ============================================================================
def plot_long_term_scissors(annual: pd.DataFrame, models: dict,
                            savepath: str = "fig5a_taiwan_scissors.png") -> None:
    """雙軸：尖峰負載（上升紅）× 備用容量率（下降綠）。"""
    fig, ax1 = plt.subplots(figsize=(11.5, 5.5))

    col_load = "#B22222"
    ax1.plot(annual["AD_year"], annual["peak_load_MW"] / 1000,
             color=col_load, lw=2.0, marker="o", ms=3.5,
             label="尖峰負載")
    # 迴歸線
    x = annual["AD_year"].values
    y_pred_load = models["load_model"].predict(sm.add_constant(x)) / 1000
    ax1.plot(x, y_pred_load, color=col_load, lw=1.2, ls="--", alpha=0.6,
             label="OLS 趨勢線")
    ax1.set_xlabel("西元年")
    ax1.set_ylabel("尖峰負載（GW ／ 十億瓦）", color=col_load)
    ax1.tick_params(axis="y", labelcolor=col_load)

    ax2 = ax1.twinx()
    col_margin = "#006400"
    ax2.plot(annual["AD_year"], annual["reserve_margin_pct"],
             color=col_margin, lw=2.0, marker="s", ms=3.5,
             label="備用容量率")
    y_pred_margin = models["margin_model"].predict(sm.add_constant(x))
    ax2.plot(x, y_pred_margin, color=col_margin, lw=1.2, ls="--", alpha=0.6,
             label="OLS 趨勢線")
    ax2.set_ylabel("備用容量率（%）", color=col_margin)
    ax2.tick_params(axis="y", labelcolor=col_margin)
    ax2.axhline(15, color="gray", lw=0.7, ls=":", alpha=0.6)
    ax2.text(annual["AD_year"].max(), 15.5, "15% 規劃基準",
             fontsize=8, color="gray", ha="right")

    # 標註民國 111 年（2022）資料定義變動
    ax2.axvline(2022, color="blue", lw=1.0, ls=":", alpha=0.5)
    ax2.annotate(
        "2022 起\n改計夜間\n備用容量率",
        xy=(2022, annual[annual["AD_year"] == 2022]["reserve_margin_pct"].iloc[0]
            if 2022 in annual["AD_year"].values else 7),
        xytext=(2014, 25), fontsize=8, color="blue",
        arrowprops=dict(arrowstyle="->", color="blue", lw=0.8),
    )

    # 標題與圖例
    beta_load = models["load_model"].params[1] * 10
    beta_margin = models["margin_model"].params[1] * 10
    ax1.set_title(
        f"圖 5a｜台灣電力剪刀效應（{annual['AD_year'].min()}–"
        f"{annual['AD_year'].max()}）\n"
        f"尖峰負載每十年 +{beta_load:.0f} MW｜"
        f"備用容量率每十年 {beta_margin:+.1f} 百分點"
    )

    lines1, labels1 = ax1.get_legend_handles_labels()
    lines2, labels2 = ax2.get_legend_handles_labels()
    ax1.legend(lines1 + lines2, labels1 + labels2,
               loc="upper left", frameon=False, fontsize=9)
    ax1.grid(True, alpha=0.25)
    fig.tight_layout()
    fig.savefig(savepath, dpi=200, bbox_inches="tight")
    print(f"[已輸出] {savepath}")
    plt.close(fig)


# ============================================================================
# 第 7 節｜視覺化（二）：夏季每日備轉容量熱地圖
# ============================================================================
def plot_summer_heatmap(daily: pd.DataFrame,
                        savepath: str = "fig5b_summer_heatmap.png") -> None:
    """
    熱地圖：y 軸年份、x 軸 6–9 月每日，色階表示備轉容量率。
    使用台電燈號制作 4 色離散色階，能直觀看出哪些日子「吃緊」。
    """
    summer = daily[daily["month"].between(6, 9)].copy()
    summer["md"] = summer["date"].dt.strftime("%m-%d")

    # 構造年 × 日 的矩陣
    pivot = summer.pivot_table(index="year", columns="md",
                                values="operating_reserve_pct",
                                aggfunc="mean")
    pivot = pivot.sort_index(ascending=False)  # 較新年在上

    # 燈號離散色階
    cmap = ListedColormap(["#CC0000", "#FF8800", "#FFD700", "#228B22"])
    # 紅、橘、黃、綠（由低到高）
    bounds = [0, 3, 6, 10, 30]
    norm = BoundaryNorm(bounds, cmap.N)

    fig, ax = plt.subplots(figsize=(13, 3.5 + 0.4 * len(pivot.index)))
    im = ax.imshow(pivot.values, aspect="auto", cmap=cmap, norm=norm,
                   interpolation="nearest")

    # 標籤
    ax.set_yticks(range(len(pivot.index)))
    ax.set_yticklabels(pivot.index)
    # x 軸僅顯示每月 1 日 + 15 日
    xticks_pos, xticks_lab = [], []
    for i, md in enumerate(pivot.columns):
        if md.endswith("-01") or md.endswith("-15"):
            xticks_pos.append(i)
            xticks_lab.append(md)
    ax.set_xticks(xticks_pos)
    ax.set_xticklabels(xticks_lab, rotation=45, ha="right", fontsize=8)

    ax.set_xlabel("日期（月-日）")
    ax.set_ylabel("年份")
    ax.set_title("圖 5b｜夏季（6–9 月）每日備轉容量率熱地圖\n"
                 "紅：<3% 限電警戒｜橘：3–6% 供電警戒｜黃：6–10% 供電吃緊｜綠：≥10% 充裕")

    # 色條
    cbar = fig.colorbar(im, ax=ax, shrink=0.7, ticks=[1.5, 4.5, 8, 20])
    cbar.ax.set_yticklabels(["紅燈", "橘燈", "黃燈", "綠燈"])
    cbar.set_label("台電燈號", fontsize=9)

    fig.tight_layout()
    fig.savefig(savepath, dpi=200, bbox_inches="tight")
    print(f"[已輸出] {savepath}")
    plt.close(fig)


# ============================================================================
# 第 8 節｜視覺化（三）：夏季分布小提琴圖
# ============================================================================
def plot_summer_distribution(daily: pd.DataFrame,
                             savepath: str = "fig5c_summer_dist.png") -> None:
    """分年度比較夏季備轉容量率分布。"""
    summer = daily[daily["month"].between(6, 9)].copy()
    years = sorted(summer["year"].unique())

    groups = [summer[summer["year"] == y]["operating_reserve_pct"].values
              for y in years]

    fig, ax = plt.subplots(figsize=(9, 5))
    parts = ax.violinplot(groups, positions=range(len(years)),
                          widths=0.7, showmeans=True, showextrema=False)
    cmap = mpl.colormaps["YlOrRd"]
    for i, pc in enumerate(parts["bodies"]):
        pc.set_facecolor(cmap(0.3 + i * 0.2))
        pc.set_alpha(0.7)
    parts["cmeans"].set_color("black")

    # 警戒線
    ax.axhline(10, color="#228B22", lw=1.0, ls="--", alpha=0.7,
               label="綠燈門檻 10%")
    ax.axhline(6, color="#FFD700", lw=1.0, ls="--", alpha=0.7,
               label="黃燈門檻 6%")
    ax.axhline(3, color="#FF8800", lw=1.0, ls="--", alpha=0.7,
               label="橘燈門檻 3%")

    ax.set_xticks(range(len(years)))
    ax.set_xticklabels([f"{y}\n夏季" for y in years])
    ax.set_ylabel("備轉容量率（%）")
    ax.set_title(f"圖 5c｜近 {len(years)} 年夏季備轉容量率分布比較")
    ax.legend(loc="upper right", frameon=False, fontsize=8)
    ax.grid(True, axis="y", alpha=0.3)
    fig.tight_layout()
    fig.savefig(savepath, dpi=200, bbox_inches="tight")
    print(f"[已輸出] {savepath}")
    plt.close(fig)


# ============================================================================
# 第 9 節｜主流程
# ============================================================================
def main() -> None:
    print("=" * 72)
    print("模組五｜台灣夏季備轉容量率惡化趨勢")
    print("=" * 72)

    # 1. 下載資料
    print("\n[步驟 1/5] 下載年度資料（自 1982 年起）...")
    annual = load_annual_reserve()

    print("\n[步驟 2/5] 下載每日資料（近三年）...")
    daily = load_daily_operating_reserve()

    # 2. 長期趨勢分析
    print("\n[步驟 3/5] 長期剪刀效應分析...")
    models = long_term_scissors_regression(annual)

    # 3. 夏季警戒日分析
    print("\n[步驟 4/5] 夏季警戒日統計...")
    alert_summary = summer_alert_days(daily)

    # 4. 視覺化
    print("\n[步驟 5/5] 產出視覺化...")
    plot_long_term_scissors(annual, models)
    plot_summer_heatmap(daily)
    plot_summer_distribution(daily)

    print("\n[完成] 主流程結束。共產出三張圖：")
    print("  - fig5a_taiwan_scissors.png  （長期剪刀效應）")
    print("  - fig5b_summer_heatmap.png   （夏季熱地圖）")
    print("  - fig5c_summer_dist.png      （年度分布比較）")
    print("\n※ 本程式每次執行都會下載最新資料。若台電更新，"
          "\n  下次執行即自動反映。")


# ============================================================================
# 第 10 節｜課堂討論題
# ============================================================================
"""
【討論題一｜跨模組整合】
  Q1. 模組一（京都櫻花）的剪刀效應是「氣溫↑ × 花期↓」，
      模組五（台灣電力）的剪刀效應是「負載↑ × 容量率↓」。
      兩個剪刀效應之間有什麼因果關聯？（提示：冷氣用電 → 夏季負載
      → 發電廠冷卻水效率 → 容量率下降 → 更熱的日子更缺電）
      這個「氣候—基礎設施—社會」的正回饋鏈，如何讓單純的氣候議題
      變成能源治理議題？

【討論題二｜資料定義變動】
  Q2. 台電自民國 111 年（2022）起，「備用容量率」改以「夜間備用
      容量率」呈現，原因是白天太陽光電供給充足、夜間才是真正的供電
      壓力點。這個變動對本模組的時序圖會造成什麼解讀挑戰？如果你是
      研究者，應該如何在論文中處理這段「定義不連續」的資料？

【討論題三｜指標選擇的政治性】
  Q3. 備用容量率（年度、規劃用）與備轉容量率（每日、運轉用）在
      公共討論中常被混用。「夏季用電吃緊」究竟該用哪個指標較精確？
      試比較執政黨與在野黨在同一年度會如何各自選擇有利自己論述的
      指標，並批判性地評估兩種指標的適用場合。

【討論題四｜因果 vs 相關】
  Q4. 尖峰負載年年上升，但上升的原因不只是氣溫升高：人口增加、
      工業用電成長（特別是半導體產業）、電動車普及都是驅動因素。
      若要嚴謹分離「氣候變遷對尖峰負載的貢獻」，你需要哪些額外
      資料？（提示：每日最高溫、工業用電、人口、電動車數量、
      冷氣普及率）

【討論題五｜行動與政策】
  Q5. 看完圖 5b 熱地圖後，你最想對政府提的一個電力政策建議是什麼？
      這個建議會如何影響（a）你家電費、（b）空氣品質、（c）氣候
      減緩目標、（d）科技業競爭力？政策選擇永遠在多個價值之間取捨，
      寫下你的價值優先順序並說明理由。
"""


if __name__ == "__main__":
    main()
