# -*- coding: utf-8 -*-
"""
===========================================================================
模組四｜死亡螺旋（The Death Spiral）
北極海冰面積的時間軌跡（1979–至今）
===========================================================================

【教學目標｜ILO】
  1. 學習者能自動下載 NSIDC 十二個月份 CSV 並合併為長格式時序。
  2. 學習者能將直角坐標數據轉換為極坐標，並理解坐標系選擇的敘事差異。
  3. 學習者能使用 matplotlib 繪製具有時間色彩映射的極坐標螺旋圖。
  4. 學習者能以 OLS 迴歸量化九月海冰極小值的長期趨勢，並推估
     「無冰北極」（ice-free Arctic）的可能年份。

【資料來源】
  NSIDC Sea Ice Index v4（北半球月均海冰面積）
  DOI: 10.7265/a98x-0f50
  URL: https://noaadata.apps.nsidc.org/NOAA/G02135/north/monthly/data/
       檔案命名：N_{01..12}_extent_v4.0.csv
  起始年：1978 年 11 月（1979 年 1 月始為完整年份）

【視覺化靈感來源】
  - Hawkins, E. (2016). Climate spirals. Climate Lab Book.
  - Pluck, K. (2017). Arctic sea ice death spiral.

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

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

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

# ============================================================================
# 第 0 節｜環境設定
# ============================================================================
import io
import warnings
import urllib.request
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib import rcParams
from matplotlib.collections import LineCollection
from matplotlib.colors import Normalize
from matplotlib.cm import ScalarMappable
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 節｜下載並合併 12 個月份的海冰面積 CSV
# ============================================================================
def load_nsidc_monthly() -> pd.DataFrame:
    """
    從 NSIDC 下載北半球月均海冰面積（extent）資料。
    每個月份一個檔案，共 12 檔。

    欄位格式（NSIDC v4）：year, mo, data-type, region, extent, area
      extent: 海冰範圍（包括濃度 ≥15% 的所有海洋格）
      area:   海冰面積（不含北極點觀測盲區）
      單位：百萬平方公里（10⁶ km²）

    missing 值為 -9999，本函式會過濾掉。
    """
    base = ("https://noaadata.apps.nsidc.org/NOAA/G02135/north/"
            "monthly/data/N_{mo:02d}_extent_v4.0.csv")

    all_dfs = []
    print("[下載 NSIDC 月資料] ", end="", flush=True)
    for mo in range(1, 13):
        url = base.format(mo=mo)
        try:
            df = pd.read_csv(url, skipinitialspace=True)
            # 清理欄位名（NSIDC 檔案欄名有空白）
            df.columns = [c.strip().replace("-", "_").lower() for c in df.columns]
            # 欄位 ['year', 'mo', 'data_type', 'region', 'extent', 'area']
            df = df[["year", "mo", "extent"]].copy()
            df["extent"] = pd.to_numeric(df["extent"], errors="coerce")
            df = df[df["extent"] > 0]  # 過濾 -9999 與 NaN
            all_dfs.append(df)
            print(f"{mo:02d}", end=" ", flush=True)
        except Exception as e:
            print(f"\n  [警告] {mo:02d} 月讀取失敗：{e}")
    print("✓")

    combined = pd.concat(all_dfs, ignore_index=True)
    combined = combined.sort_values(["year", "mo"]).reset_index(drop=True)
    combined["year"] = combined["year"].astype(int)
    combined["mo"] = combined["mo"].astype(int)

    print(f"[NSIDC] 共讀入 {len(combined):,} 筆月資料；"
          f"年度範圍 {combined['year'].min()}–{combined['year'].max()}")
    return combined


# ============================================================================
# 第 2 節｜敘述統計與季節循環
# ============================================================================
def seasonal_cycle(df: pd.DataFrame) -> pd.DataFrame:
    """計算各月份的長期平均與標準差。"""
    seasonal = df.groupby("mo")["extent"].agg(
        mean="mean", std="std",
        min_year=lambda s: df.loc[s.idxmin(), "year"] if len(s) else np.nan,
        min_val="min",
        max_year=lambda s: df.loc[s.idxmax(), "year"] if len(s) else np.nan,
        max_val="max",
    )
    print("\n─── 各月海冰面積（百萬 km²）───")
    print(seasonal.round(2).to_string())
    print("\n※ 最大值通常出現在 3 月（冬末），最小值在 9 月（夏末）。")
    return seasonal


# ============================================================================
# 第 3 節｜直角坐標先行：標準時序圖（作為對比）
# ============================================================================
def plot_cartesian_baseline(df: pd.DataFrame,
                            savepath: str = "fig4a_cartesian.png") -> None:
    """
    傳統直角坐標時序圖。作為「為什麼我們需要極坐標」的對比。
    """
    df = df.copy()
    df["date"] = pd.to_datetime(
        df["year"].astype(str) + "-" + df["mo"].astype(str) + "-15")

    fig, ax = plt.subplots(figsize=(12, 4.5))
    ax.plot(df["date"], df["extent"], color="#1f77b4", lw=0.8, alpha=0.8)
    ax.set_xlabel("年份")
    ax.set_ylabel("海冰面積（百萬 km²）")
    ax.set_title("圖 4a｜直角坐標版本：可見鋸齒狀季節循環 + 長期緩降趨勢")
    ax.grid(True, alpha=0.3)
    fig.tight_layout()
    fig.savefig(savepath, dpi=200, bbox_inches="tight")
    print(f"[已輸出] {savepath}")
    plt.close(fig)


# ============================================================================
# 第 4 節｜極坐標轉換：數學原理
# ============================================================================
def prepare_polar_coordinates(df: pd.DataFrame) -> pd.DataFrame:
    """
    將月份 (mo) 轉換為角度 (theta)：
        theta = 2π × (mo - 1) / 12
    使 1 月 = 0 rad（在 matplotlib polar 中預設為東方，即 3 點鐘方向），
    但習慣上我們希望 1 月朝上（12 點鐘），因此後續在繪圖時會旋轉。

    半徑 (r) = 海冰面積，直接映射。
    """
    df = df.copy()
    df["theta"] = 2 * np.pi * (df["mo"] - 1) / 12
    df["r"] = df["extent"]
    return df


# ============================================================================
# 第 5 節｜死亡螺旋主視覺（Hawkins–Pluck 風格）
# ============================================================================
def plot_death_spiral(df: pd.DataFrame,
                      savepath: str = "fig4b_death_spiral.png",
                      cmap_name: str = "viridis") -> None:
    """
    核心視覺化：
      - 極坐標（polar）
      - 1 月朝上；順時針（模擬日曆走法）
      - 線段顏色依「所屬年份」漸變（viridis：深藍 → 亮黃）
      - 徑向軸刻度表示海冰面積
    """
    df = prepare_polar_coordinates(df).sort_values(["year", "mo"])

    # 色彩映射：1979（深藍）→ 最新年（亮黃）
    year_min, year_max = df["year"].min(), df["year"].max()
    norm = Normalize(vmin=year_min, vmax=year_max)
    cmap = mpl.colormaps[cmap_name]

    # 建立 polar axes
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection="polar")
    ax.set_theta_zero_location("N")  # 0° 朝上（1 月）
    ax.set_theta_direction(-1)        # 順時針

    # 把「年份」連續化：用 year + (mo-1)/12 做漸變顏色起點
    # 每段線使用起點年份的色
    segments = []
    colors = []
    arr = df[["theta", "r", "year", "mo"]].values
    for i in range(len(arr) - 1):
        t0, r0, y0, m0 = arr[i]
        t1, r1, y1, m1 = arr[i + 1]
        # 只連接相鄰的月份（避免跨年切換角度時產生大跨角連線問題）
        # 雖然 theta 從 12 月(330°) → 次年 1 月(0°) 會「回捲」，
        # matplotlib polar 會自動處理最短路徑，因此連線仍會穿過 0°。
        segments.append([(t0, r0), (t1, r1)])
        colors.append(cmap(norm(y0 + (m0 - 1) / 12)))

    lc = LineCollection(segments, colors=colors, linewidths=1.3, alpha=0.85)
    ax.add_collection(lc)

    # 徑向軸設定
    r_min, r_max = df["r"].min() * 0.9, df["r"].max() * 1.05
    ax.set_rmin(r_min)
    ax.set_rmax(r_max)
    ax.set_rticks([4, 8, 12, 16])  # 徑向刻度（百萬 km²）
    ax.set_rlabel_position(135)     # 將徑向標籤放在西南方向
    ax.tick_params(axis="y", labelsize=8, colors="#666666")

    # 月份標籤
    month_labels = ["一月", "二月", "三月", "四月", "五月", "六月",
                    "七月", "八月", "九月", "十月", "十一月", "十二月"]
    ax.set_xticks(np.linspace(0, 2 * np.pi, 12, endpoint=False))
    ax.set_xticklabels(month_labels, fontsize=10)
    ax.tick_params(axis="x", pad=10)

    # 色條（年份圖例）
    sm_ = ScalarMappable(norm=norm, cmap=cmap)
    sm_.set_array([])
    cbar = fig.colorbar(sm_, ax=ax, pad=0.12, shrink=0.7,
                        orientation="vertical")
    cbar.set_label("年份", fontsize=10)
    cbar.ax.tick_params(labelsize=8)

    # 標題
    ax.set_title(
        f"圖 4b｜北極海冰死亡螺旋（{year_min}–{year_max}）\n"
        f"內旋＝海冰退縮；顏色由深至亮＝時間推進\n"
        f"資料：NSIDC Sea Ice Index v4；視覺設計啟發自 Hawkins (2016)、Pluck (2017)",
        fontsize=11, pad=22,
    )

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


# ============================================================================
# 第 6 節｜九月極小值趨勢分析（核心量化證據）
# ============================================================================
def september_trend(df: pd.DataFrame) -> dict:
    """
    九月是北極海冰的年度最小值月份，最能代表「夏季極端狀態」。
    以 OLS 迴歸估算下降速率與 R²。
    """
    sept = df[df["mo"] == 9].copy().sort_values("year").reset_index(drop=True)

    X = sm.add_constant(sept["year"].values)
    y = sept["extent"].values
    model = sm.OLS(y, X).fit()
    slope = model.params[1]
    intercept = model.params[0]
    r2 = model.rsquared
    p_val = model.pvalues[1]

    # 每十年下降量
    per_decade = slope * 10

    print("\n─── 九月海冰極小值 OLS 迴歸 ───")
    print(f"樣本數 n = {len(sept)}")
    print(f"斜率 = {slope:+.4f} 百萬 km² / 年")
    print(f"每十年下降 = {per_decade:+.3f} 百萬 km²")
    print(f"R² = {r2:.3f}  (p = {p_val:.2e})")

    # 線性外推：預測何年海冰 < 1 百萬 km²（「無冰北極」慣用門檻）
    # 海冰 = intercept + slope × year → year = (1.0 - intercept) / slope
    year_ice_free = (1.0 - intercept) / slope
    print(f"\n※ 線性外推「無冰北極」（<1 百萬 km²）預期年份：{year_ice_free:.0f}")
    print("  （警語：此為純線性外推；實際氣候系統具非線性反饋，")
    print("   較保守的學術推估為 2030–2050 間；IPCC AR6 指出 2050 前至少")
    print("   會出現一次九月無冰事件）")

    return {
        "sept_df": sept,
        "slope": slope, "intercept": intercept,
        "r_squared": r2, "p_value": p_val,
        "per_decade": per_decade,
        "year_ice_free_linear": year_ice_free,
        "model": model,
    }


def plot_september_trend(result: dict,
                         savepath: str = "fig4c_september_trend.png") -> None:
    """九月極小值散點 + 迴歸線 + 外推虛線。"""
    sept = result["sept_df"]
    slope, intercept = result["slope"], result["intercept"]

    fig, ax = plt.subplots(figsize=(11, 5))

    # 觀測散點
    ax.scatter(sept["year"], sept["extent"],
               s=40, c="#1f4e79", alpha=0.8, edgecolor="white",
               label="NSIDC 觀測")

    # 迴歸線（觀測期內）
    x_obs = np.array([sept["year"].min(), sept["year"].max()])
    ax.plot(x_obs, intercept + slope * x_obs,
            color="#C71585", lw=2.5, label="OLS 迴歸線")

    # 外推虛線（延伸到無冰年份或 2060，取近者）
    x_ext_end = min(result["year_ice_free_linear"], 2060)
    x_ext = np.array([sept["year"].max(), x_ext_end])
    ax.plot(x_ext, intercept + slope * x_ext,
            color="#C71585", lw=2.0, ls="--", alpha=0.6,
            label="線性外推（僅作教學示意）")

    # 標註關鍵年份
    min_idx = sept["extent"].idxmin()
    min_year = int(sept.loc[min_idx, "year"])
    min_val = sept.loc[min_idx, "extent"]
    ax.annotate(
        f"{min_year} 年\n歷史最低\n{min_val:.2f} M km²",
        xy=(min_year, min_val),
        xytext=(min_year - 10, min_val - 1.5),
        fontsize=9, color="#B22222",
        arrowprops=dict(arrowstyle="->", color="#B22222", lw=1.2),
    )

    # 1 百萬 km² 無冰門檻
    ax.axhline(1.0, color="red", lw=1.0, ls=":", alpha=0.7,
               label="「無冰北極」門檻（<1 M km²）")

    ax.set_xlabel("年份")
    ax.set_ylabel("九月海冰範圍（百萬 km²）")
    ax.set_title(
        f"圖 4c｜九月北極海冰極小值之線性下降\n"
        f"每十年下降 {result['per_decade']:.3f} M km²，"
        f"R² = {result['r_squared']:.2f}，p < .001"
    )
    ax.legend(loc="upper right", frameon=False, fontsize=9)
    ax.grid(True, alpha=0.3)
    ax.set_ylim(bottom=0)
    fig.tight_layout()
    fig.savefig(savepath, dpi=200, bbox_inches="tight")
    print(f"[已輸出] {savepath}")
    plt.close(fig)


# ============================================================================
# 第 7 節｜十年平均比較（視覺化「退縮的速度」）
# ============================================================================
def decade_comparison(df: pd.DataFrame,
                      savepath: str = "fig4d_decade_comparison.png") -> None:
    """
    依十年分組，繪製各十年的「平均年循環」以顯示整體退縮。
    每條線為一個十年（1980s, 1990s, 2000s, 2010s, 2020s）。
    """
    df = df.copy()
    df["decade"] = (df["year"] // 10) * 10

    fig, ax = plt.subplots(figsize=(10, 5.5))
    cmap = mpl.colormaps["viridis"]
    decades = sorted(df["decade"].unique())
    n = len(decades)

    for i, dec in enumerate(decades):
        sub = df[df["decade"] == dec].groupby("mo")["extent"].mean()
        # 跳過資料太少的年代（如 2020s 可能不滿十年）
        label = f"{dec}s"
        if dec == decades[-1]:
            end_year = df[df["decade"] == dec]["year"].max()
            label = f"{dec}–{end_year}"
        ax.plot(sub.index, sub.values,
                color=cmap(i / max(n - 1, 1)), lw=2.2, marker="o",
                label=label)

    ax.set_xticks(range(1, 13))
    ax.set_xticklabels(["1", "2", "3", "4", "5", "6",
                        "7", "8", "9", "10", "11", "12"])
    ax.set_xlabel("月份")
    ax.set_ylabel("平均海冰範圍（百萬 km²）")
    ax.set_title("圖 4d｜各十年的平均年循環：\n"
                 "每條線向下位移 = 整個十年的海冰都比上一個十年少")
    ax.legend(loc="lower center", ncol=len(decades), frameon=False, fontsize=9)
    ax.grid(True, alpha=0.3)
    fig.tight_layout()
    fig.savefig(savepath, dpi=200, bbox_inches="tight")
    print(f"[已輸出] {savepath}")
    plt.close(fig)


# ============================================================================
# 第 8 節｜主流程
# ============================================================================
def main() -> None:
    print("=" * 72)
    print("模組四｜死亡螺旋：北極海冰的時間軌跡")
    print("=" * 72)

    df = load_nsidc_monthly()
    seasonal_cycle(df)

    # 五張圖：直角、螺旋、九月趨勢、十年比較
    plot_cartesian_baseline(df)
    plot_death_spiral(df)

    sept_result = september_trend(df)
    plot_september_trend(sept_result)

    decade_comparison(df)

    print("\n[完成] 主流程結束。共產出四張圖：")
    print("  - fig4a_cartesian.png         （對比基線）")
    print("  - fig4b_death_spiral.png      （主視覺）")
    print("  - fig4c_september_trend.png   （量化證據）")
    print("  - fig4d_decade_comparison.png （退縮速度）")


# ============================================================================
# 第 9 節｜課堂討論題
# ============================================================================
"""
【討論題一｜坐標系的敘事差異】
  Q1. 同一份資料，直角坐標版本（圖 4a）與極坐標版本（圖 4b）講的是
      不同的故事。用兩到三句話描述：直角版最突出什麼？極坐標版最突出
      什麼？為什麼 Hawkins 與 Pluck 選擇用極坐標來做氣候視覺化？

【討論題二｜色彩映射的倫理】
  Q2. 本模組使用 viridis 色盤（深藍→亮黃）。若改用紅色系（深紅→亮紅）
      或溫度系（藍→紅），整個圖的情緒會變成什麼？Knaflic（2015）強調
      「色彩是資料敘事的情感引擎」——你認為哪種配色最能忠實呈現事實
      而不過度煽動？

【討論題三｜線性外推的限度】
  Q3. 程式計算出「線性外推無冰年份」通常落在 2060–2080 之間。但 IPCC
      AR6 報告指出 2050 前就可能出現第一次九月無冰事件。為什麼實際
      科學推估比純線性外推「更悲觀」？（提示：反照率正回饋、北極放大
      效應、冰層厚度非線性崩潰）

【討論題四｜超物件的視覺化】
  Q4. Morton（2013）所稱的「超物件」在時間上超越個人生命、在空間上
      超越個人感知。死亡螺旋如何「把超物件縮尺到一張圖」？這樣的縮尺
      有什麼風險？（提示：當超物件被「馴化」為圖表，會不會讓人誤以為
      它變得可控？）

【討論題五｜動手的學習意義】
  Q5. 本模組要求學生親手跑程式碼、調整參數（例如改變 cmap_name）。
      請寫下你在執行程式時遇到的至少一個錯誤（環境問題、資料問題、
      圖形看起來不對），並描述你是如何解決的。這段「除錯經驗」與
      Ridsdale 等人（2015）所說「資料操作是日常流程」有何呼應？
"""


if __name__ == "__main__":
    main()