# -*- coding: utf-8 -*-
"""
===========================================================================
模組二｜歷史斷層（The Breakpoint）
勃艮第葡萄採收日 × 全球氣溫異常值
===========================================================================

【教學目標｜ILO】
  1. 學習者能讀取 NOAA 古氣候資料庫的純文字古物候資料。
  2. 學習者能以 t 檢定與 CUSUM 法偵測時間序列中的「結構突變」。
  3. 學習者能區分「漸進趨勢」（模組一）與「斷層變化」（本模組）
     兩種統計結構，並理解 1988 年作為全球氣候斷層的意義。
  4. 學習者能批判性檢視歷史資料的來源偏誤（本模組以 1540 年
     Chuine 2004 與 Labbé-Gaveau 2011 的資料分歧為例）。

【資料來源】
  (1) 勃艮第葡萄採收日：Chuine, I. et al. (2004), Nature, 432, 289–290.
      DOI: 10.25921/e3mc-wm08
      URL: https://www.ncei.noaa.gov/pub/data/paleo/historical/france/
            burgundy2004.txt
      範圍：1370–2003 CE（634 年）；以 9 月 1 日為 0 起算之天數

  (2) 全球氣溫異常值：NASA GISTEMP v4（陸地+海洋表面溫度）
      URL：https://data.giss.nasa.gov/gistemp/tabledata_v4/GLB.Ts+dSST.csv
      基期：1951–1980；範圍：1880 年迄今

【延伸閱讀（教師備課用）】
  Labbé, T. et al. (2019). The longest homogeneous series of grape harvest
  dates, Beaune 1354–2018. Climate of the Past, 15, 1485–1501.
  此研究發現：1354–1987 年平均採收日為 9/28（DOY 271），
              1988–2018 年平均採收日為 9/15（DOY 258），
              整整提前 13 天，構成清晰的「1988 歷史斷層」。

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

【Colab 首次執行】
  !apt-get install -y fonts-noto-cjk -q   # 中文字型（安裝後需 Runtime → Restart）

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

# ============================================================================
# 第 0 節｜環境設定
# ============================================================================
import io
import warnings
import urllib.request
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
from matplotlib.patches import Rectangle
from scipy import stats
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 節｜讀取勃艮第葡萄採收日資料（Chuine et al. 2004，NOAA 古氣候庫）
# ============================================================================
def load_burgundy_ghd() -> pd.DataFrame:
    """
    讀取 NOAA 存檔之勃艮第葡萄採收日資料（1370–2003）。

    原始檔欄位：
        Column 1: Year AD
        Column 2: Harvest date（以 9 月 1 日為第 0 天起算之天數）

    本函式將原始「9/1 起算天數」轉換為標準年積日（DOY）：
        DOY = 244 + Harvest_date_days   （9 月 1 日 ≈ DOY 244，平年）
    """
    url = ("https://www.ncei.noaa.gov/pub/data/paleo/historical/"
           "france/burgundy2004.txt")

    # 原始檔有大量 header 說明，資料從 "Year     Harvest date" 下一列開始
    with urllib.request.urlopen(url, timeout=30) as resp:
        text = resp.read().decode("latin-1")  # NOAA 檔案常見 latin-1 編碼

    # 定位資料起點：找到 "Year" 與 "Harvest" 皆出現的表頭列
    lines = text.splitlines()
    start = None
    for i, ln in enumerate(lines):
        if "Year" in ln and "Harvest" in ln:
            start = i + 1
            break
    if start is None:
        raise RuntimeError("無法在 NOAA 檔案中定位資料起點。")

    # 從表頭下一列開始解析；容許空白或多重分隔
    rows = []
    for ln in lines[start:]:
        parts = ln.strip().split()
        if len(parts) >= 2:
            try:
                year = int(parts[0])
                days_after_sep1 = float(parts[1])
                rows.append((year, days_after_sep1))
            except ValueError:
                continue  # 非資料列

    df = pd.DataFrame(rows, columns=["AD", "DaysAfterSep1"])
    df["HarvestDOY"] = 244 + df["DaysAfterSep1"]  # 轉為年積日

    print(f"[勃艮第葡萄] 共讀入 {len(df):,} 筆紀錄；"
          f"年代範圍 {df['AD'].min()}–{df['AD'].max()}")
    print(f"[勃艮第葡萄] 採收日 DOY 範圍："
          f"{df['HarvestDOY'].min():.0f}–{df['HarvestDOY'].max():.0f}"
          f"（約 8 月底至 10 月中）")
    return df.reset_index(drop=True)


# ============================================================================
# 第 2 節｜讀取全球氣溫異常值（NASA GISTEMP v4）
# ============================================================================
def load_gistemp() -> pd.DataFrame:
    """讀取 NASA GISTEMP 陸地+海洋年均溫異常值。基期：1951–1980；單位：°C。"""
    url = ("https://data.giss.nasa.gov/gistemp/tabledata_v4/"
           "GLB.Ts+dSST.csv")
    df = pd.read_csv(url, skiprows=1)
    df = df[["Year", "J-D"]].rename(
        columns={"Year": "AD", "J-D": "TempAnomaly"})
    df["TempAnomaly"] = pd.to_numeric(df["TempAnomaly"], errors="coerce")
    df = df.dropna().copy()
    df["AD"] = df["AD"].astype(int)
    print(f"[全球氣溫] 共讀入 {len(df):,} 筆年度紀錄；"
          f"範圍 {df['AD'].min()}–{df['AD'].max()}")
    return df.reset_index(drop=True)


# ============================================================================
# 第 3 節｜結構突變偵測（一）：先驗斷點 t 檢定
# ============================================================================
def breakpoint_ttest(df: pd.DataFrame,
                     break_year: int = 1988,
                     value_col: str = "HarvestDOY") -> dict:
    """
    以指定斷點年份切割序列，並對斷點前後兩段做 Welch's t 檢定
    （不假設變異數齊性）。

    教學意義：此法測試「以 1988 為界的兩段均值差異是否顯著」，
    對應 Labbé et al. (2019) 的核心主張。
    """
    before = df[df["AD"] < break_year][value_col].values
    after = df[df["AD"] >= break_year][value_col].values

    t_stat, p_val = stats.ttest_ind(before, after, equal_var=False)
    shift_days = before.mean() - after.mean()

    # Cohen's d 效應量
    pooled_sd = np.sqrt((before.std(ddof=1) ** 2 + after.std(ddof=1) ** 2) / 2)
    cohens_d = shift_days / pooled_sd

    print(f"\n─── 斷點 t 檢定（H0：{break_year} 年前後均值相等）───")
    print(f"前段（{df['AD'].min()}–{break_year-1}）： "
          f"n = {len(before)},  均值 = {before.mean():.2f},  "
          f"SD = {before.std(ddof=1):.2f}")
    print(f"後段（{break_year}–{df['AD'].max()}）： "
          f"n = {len(after)},  均值 = {after.mean():.2f},  "
          f"SD = {after.std(ddof=1):.2f}")
    print(f"平均提前日數 ≈ {shift_days:.2f} 天")
    print(f"t = {t_stat:.3f},  p = {p_val:.2e}")
    print(f"Cohen's d = {cohens_d:.2f}  "
          f"({'微小' if abs(cohens_d) < 0.2 else '小' if abs(cohens_d) < 0.5 else '中' if abs(cohens_d) < 0.8 else '大'}效應)")

    return {
        "break_year": break_year,
        "mean_before": before.mean(), "mean_after": after.mean(),
        "sd_before": before.std(ddof=1), "sd_after": after.std(ddof=1),
        "shift_days": shift_days,
        "t_stat": t_stat, "p_value": p_val, "cohens_d": cohens_d,
    }


# ============================================================================
# 第 4 節｜結構突變偵測（二）：CUSUM 無先驗掃描
# ============================================================================
def cusum_scan(df: pd.DataFrame,
               value_col: str = "HarvestDOY",
               min_segment: int = 30) -> dict:
    """
    CUSUM（cumulative sum）法：**不預設斷點位置**，以暴力搜索找到
    使得前後段變異數總和最小的斷點年份。

    教學意義：這是「資料自己說話」的作法——若我們不先告訴程式 1988 年
    重要，演算法是否仍會自動找到 1988 附近？結果可讓學生體會「信號
    偵測」與「先驗假設」的不同。
    """
    years = df["AD"].values
    y = df[value_col].values
    n = len(y)

    best_score = np.inf
    best_idx = None

    for i in range(min_segment, n - min_segment):
        before, after = y[:i], y[i:]
        # 計算前後段合併的組內變異（越小代表切得越好）
        score = (before.var(ddof=1) * (len(before) - 1)
                 + after.var(ddof=1) * (len(after) - 1))
        if score < best_score:
            best_score = score
            best_idx = i

    best_year = int(years[best_idx])
    mean_before = y[:best_idx].mean()
    mean_after = y[best_idx:].mean()

    print(f"\n─── CUSUM 斷點掃描（無先驗）───")
    print(f"最佳斷點年份：{best_year}")
    print(f"斷點前平均 = {mean_before:.2f} DOY")
    print(f"斷點後平均 = {mean_after:.2f} DOY")
    print(f"平均提前 = {mean_before - mean_after:.2f} 天")

    return {"best_year": best_year,
            "mean_before": mean_before, "mean_after": mean_after,
            "shift_days": mean_before - mean_after}


# ============================================================================
# 第 5 節｜極端年份分析（前 5% 最早採收）
# ============================================================================
def extreme_years_analysis(df: pd.DataFrame,
                           value_col: str = "HarvestDOY",
                           pct: float = 5.0) -> pd.DataFrame:
    """
    列出前 5% 最早採收年份，依世紀分布。
    Labbé et al. (2019) 指出：33 個極端年份中，21 個落在 1393–1719 期間；
    而自 2003 起，過去 16 年中就有 8 年屬於此類。本函式重現此分析。
    """
    threshold = np.percentile(df[value_col], pct)
    extremes = df[df[value_col] <= threshold].copy()
    extremes["Century"] = (extremes["AD"] // 100 + 1).astype(int)

    dist = extremes.groupby("Century").size().reset_index(name="n_extreme")
    print(f"\n─── 極端早採年份（前 {pct}% = DOY ≤ {threshold:.1f}）───")
    print(f"極端年份總數：{len(extremes)}")
    print("\n各世紀分布：")
    print(dist.to_string(index=False))

    # 最近 30 年（1974–2003）的極端比例
    recent_30 = extremes[extremes["AD"] >= 1974]
    print(f"\n※ 最近 30 年（1974–2003）極端年份數：{len(recent_30)}")
    print(f"※ 佔整體極端年份的 {len(recent_30) / len(extremes) * 100:.1f}%")

    return extremes


# ============================================================================
# 第 6 節｜氣溫耦合分析（現代重疊期 1880–2003）
# ============================================================================
def temperature_coupling(ghd: pd.DataFrame,
                         temp: pd.DataFrame) -> sm.regression.linear_model.RegressionResultsWrapper:
    """
    針對 1880 年後 GHD 與 GISTEMP 的重疊期做相關與 OLS 迴歸。
    """
    merged = ghd.merge(temp, on="AD", how="inner")
    print(f"\n─── 氣溫耦合（重疊期 {merged['AD'].min()}–{merged['AD'].max()}）───")
    print(f"樣本數 n = {len(merged)}")

    pearson_r, pearson_p = stats.pearsonr(
        merged["TempAnomaly"], merged["HarvestDOY"])
    print(f"Pearson  r = {pearson_r:+.3f}  (p = {pearson_p:.2e})")

    X = sm.add_constant(merged["TempAnomaly"].values)
    y = merged["HarvestDOY"].values
    model = sm.OLS(y, X).fit()
    print(f"OLS 斜率 β = {model.params[1]:+.3f} 天/°C")
    print(f"R² = {model.rsquared:.3f}")
    print(f"※ 對照模組一（京都櫻花）β ≈ −6.5 天/°C。")
    print(f"   葡萄斜率較大係因（a）葡萄對春夏溫度敏感度更高，")
    print(f"   （b）現代重疊期較短，近期升溫效應放大了斜率。")
    return model


# ============================================================================
# 第 7 節｜視覺化（一）：全景斷層圖（主圖）
# ============================================================================
def plot_breakpoint(df: pd.DataFrame,
                    ttest_result: dict,
                    savepath: str = "fig2a_breakpoint_panorama.png") -> None:
    """
    六百年全景圖，明確標出 1988 斷層：
      - 散點：每年採收日 DOY
      - 30 年移動平均：凸顯世紀尺度波動
      - 垂直紅線：1988 斷點
      - 水平帶：斷點前後兩段的平均 ± 1 SD
    """
    fig, ax = plt.subplots(figsize=(11.5, 5.5))

    # 散點
    ax.scatter(df["AD"], df["HarvestDOY"],
               s=7, c="#722F37", alpha=0.4, label="年採收日（DOY）")

    # 30 年移動平均
    roll = df["HarvestDOY"].rolling(30, min_periods=15).mean()
    ax.plot(df["AD"], roll, color="#2C1810", lw=2.2,
            label="30 年移動平均")

    # 斷點前後平均帶
    br = ttest_result["break_year"]
    xlim = (df["AD"].min() - 5, df["AD"].max() + 5)

    # 前段陰影
    y_pre = ttest_result["mean_before"]
    s_pre = ttest_result["sd_before"]
    ax.fill_between([xlim[0], br], y_pre - s_pre, y_pre + s_pre,
                    color="#8B7355", alpha=0.18,
                    label=f"1370–{br-1} 平均 ± 1 SD")
    ax.hlines(y_pre, xlim[0], br, colors="#8B7355", lw=1.8, ls="--")

    # 後段陰影
    y_post = ttest_result["mean_after"]
    s_post = ttest_result["sd_after"]
    ax.fill_between([br, xlim[1]], y_post - s_post, y_post + s_post,
                    color="#C71585", alpha=0.25,
                    label=f"{br}– 平均 ± 1 SD")
    ax.hlines(y_post, br, xlim[1], colors="#C71585", lw=1.8, ls="--")

    # 斷點垂直線
    ax.axvline(br, color="#B22222", lw=2.2, ls="-",
               label=f"{br} 斷層")

    # 斷層標註
    shift = ttest_result["shift_days"]
    ax.annotate(
        f"↓ {shift:.1f} 天\n(p = {ttest_result['p_value']:.1e},\n"
        f" Cohen's d = {ttest_result['cohens_d']:.2f})",
        xy=(br, y_post),
        xytext=(br + 25, y_post - 15),
        fontsize=10, color="#B22222",
        arrowprops=dict(arrowstyle="->", color="#B22222", lw=1.5),
    )

    ax.set_xlabel("西元年（AD）")
    ax.set_ylabel("葡萄採收日年積日（DOY）")
    ax.set_title(
        f"圖 2a｜勃艮第葡萄採收日之歷史斷層（1370–2003）\n"
        f"Chuine et al. (2004)；Labbé et al. (2019) 指出 1988 為關鍵突變點",
        pad=10,
    )
    ax.invert_yaxis()  # 採收日愈早（DOY 愈小）放在上方
    ax.set_xlim(xlim)
    ax.legend(loc="lower left", frameon=False, fontsize=9)
    ax.grid(True, alpha=0.25)
    fig.tight_layout()
    fig.savefig(savepath, dpi=200, bbox_inches="tight")
    print(f"[已輸出] {savepath}")
    plt.close(fig)


# ============================================================================
# 第 8 節｜視覺化（二）：分布比較小提琴圖
# ============================================================================
def plot_distribution_comparison(df: pd.DataFrame,
                                 break_year: int = 1988,
                                 savepath: str = "fig2b_distribution_shift.png") -> None:
    """以箱形 + 小提琴圖比較斷點前後的採收日分布，凸顯「分布平移」。"""
    before = df[df["AD"] < break_year]["HarvestDOY"].values
    after = df[df["AD"] >= break_year]["HarvestDOY"].values

    fig, ax = plt.subplots(figsize=(8, 5))
    parts = ax.violinplot([before, after], positions=[1, 2],
                          widths=0.6, showmeans=True, showmedians=False)
    for i, pc in enumerate(parts["bodies"]):
        pc.set_facecolor(["#8B7355", "#C71585"][i])
        pc.set_alpha(0.55)
    parts["cmeans"].set_color("black")

    # 箱形圖疊加
    ax.boxplot([before, after], positions=[1, 2], widths=0.18,
               patch_artist=True,
               boxprops=dict(facecolor="white", edgecolor="black"),
               medianprops=dict(color="red", lw=2))

    ax.set_xticks([1, 2])
    ax.set_xticklabels(
        [f"斷層前\n({df['AD'].min()}–{break_year-1})\nn = {len(before)}",
         f"斷層後\n({break_year}–{df['AD'].max()})\nn = {len(after)}"])
    ax.set_ylabel("葡萄採收日年積日（DOY）")
    ax.invert_yaxis()
    ax.set_title(f"圖 2b｜{break_year} 斷點前後之採收日分布比較")
    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 plot_coupled_timeseries(ghd: pd.DataFrame,
                            temp: pd.DataFrame,
                            model: sm.regression.linear_model.RegressionResultsWrapper,
                            savepath: str = "fig2c_grape_vs_temp.png") -> None:
    """現代重疊期（1880 後）的雙軸圖，對應模組一的剪刀效應。"""
    merged = ghd.merge(temp, on="AD", how="inner")
    fig, ax1 = plt.subplots(figsize=(11, 5))

    col_temp = "#B22222"
    ax1.plot(merged["AD"], merged["TempAnomaly"],
             color=col_temp, lw=1.0, alpha=0.55, label="年均溫異常")
    ax1.plot(merged["AD"],
             merged["TempAnomaly"].rolling(10, center=True, min_periods=5).mean(),
             color=col_temp, lw=2.5, label="10 年平均（溫度）")
    ax1.set_xlabel("西元年")
    ax1.set_ylabel("全球氣溫異常（°C，基期 1951–1980）", color=col_temp)
    ax1.tick_params(axis="y", labelcolor=col_temp)
    ax1.axhline(0, color=col_temp, lw=0.6, ls="--", alpha=0.4)

    ax2 = ax1.twinx()
    col_grape = "#722F37"
    ax2.plot(merged["AD"], merged["HarvestDOY"],
             color=col_grape, lw=1.0, alpha=0.55, label="葡萄採收日")
    ax2.plot(merged["AD"],
             merged["HarvestDOY"].rolling(10, center=True, min_periods=5).mean(),
             color=col_grape, lw=2.5, label="10 年平均（採收日）")
    ax2.set_ylabel("勃艮第葡萄採收日（DOY）", color=col_grape)
    ax2.tick_params(axis="y", labelcolor=col_grape)
    ax2.invert_yaxis()

    ax1.axvline(1988, color="black", lw=1.5, ls=":", alpha=0.6)
    ax1.annotate("1988", xy=(1988, ax1.get_ylim()[1]),
                 xytext=(1990, ax1.get_ylim()[1] * 0.92),
                 fontsize=9, color="black")

    ax1.set_title(
        f"圖 2c｜氣溫 × 採收日（1880–2003 重疊期）\n"
        f"OLS：β = {model.params[1]:+.2f} 天/°C，"
        f"R² = {model.rsquared:.2f}，n = {int(model.nobs)}"
    )

    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.2)
    fig.tight_layout()
    fig.savefig(savepath, dpi=200, bbox_inches="tight")
    print(f"[已輸出] {savepath}")
    plt.close(fig)


# ============================================================================
# 第 10 節｜主流程
# ============================================================================
def main() -> None:
    print("=" * 72)
    print("模組二｜歷史斷層：勃艮第葡萄 × 全球氣溫")
    print("=" * 72)

    ghd = load_burgundy_ghd()
    temp = load_gistemp()

    print("\n─── 敘述統計（全序列）───")
    print(ghd[["HarvestDOY"]].describe().round(2))

    # 斷點偵測（先驗 1988 vs CUSUM 自動掃描）
    tt = breakpoint_ttest(ghd, break_year=1988)
    cs = cusum_scan(ghd)
    if abs(cs["best_year"] - 1988) <= 5:
        print(f"\n★ CUSUM 自動掃描找到 {cs['best_year']} 年，與先驗假設的"
              f" 1988 相距僅 {abs(cs['best_year'] - 1988)} 年——資料支持"
              f"「1988 斷層」的假設。")
    else:
        print(f"\n⚠ CUSUM 找到 {cs['best_year']} 年，與 1988 差距較大，"
              f"值得進一步探討。")

    # 極端年份分析
    extreme_years_analysis(ghd)

    # 氣溫耦合
    model = temperature_coupling(ghd, temp)

    # 視覺化三張
    plot_breakpoint(ghd, tt)
    plot_distribution_comparison(ghd)
    plot_coupled_timeseries(ghd, temp, model)

    print("\n[完成] 主流程結束。")


# ============================================================================
# 第 11 節｜課堂討論題
# ============================================================================
"""
【討論題一｜統計結構】
  Q1. 模組一（京都櫻花）呈現「漸進式剪刀」，模組二（勃艮第葡萄）呈現
      「結構性斷層」。這兩種統計特徵在公共溝通上各有什麼優勢與風險？
      （提示：斷層敘事更戲劇化，但也可能被質疑為「挑選時間窗」；漸進
      敘事較穩健，但容易被「氣溫一直都在波動」的質疑稀釋。）

【討論題二｜方法論反思】
  Q2. 本模組用了兩種斷點偵測法：先驗 t 檢定（我們先告訴程式 1988 年
      重要）與 CUSUM 暴力掃描（讓程式自己找）。若 CUSUM 找到的年份
      並非 1988 而是 1975 或 1995，應如何解釋？這暴露了「資料支持
      假設」與「資料產生假設」的哪些差異？

【討論題三｜史料批判】
  Q3. Chuine et al. (2004) 對 1540 年的採收日紀錄為 10 月 4 日
      (DOY 278)；Labbé & Gaveau (2011) 重查原始檔案，發現正確日期
      應為 9 月 3 日 (DOY 247)，差距 31 天。這一錯誤在 Nature 發表
      的原始研究中流傳 15 年，直到檔案學者重新查核才被發現。
      這個案例告訴我們什麼關於「資料永遠是人類的產物」？

【討論題四｜文化感】
  Q4. 勃艮第葡萄採收日的提前，意味著葡萄在更熱的夏季成熟。釀酒學上，
      過熱的夏季會導致葡萄糖分過高、酸度不足，最終葡萄酒風味趨於
      「過度成熟」（jammy, flabby）。請討論：當「風土（terroir）」
      這個文化符碼因氣候變遷而改變時，人類無形文化遺產（UNESCO 已
      將勃艮第葡萄園列入世界遺產）是否應該被視為氣候衝擊的
      「另一種受害者」？

【討論題五｜跨模組整合】
  Q5. 請將模組一（櫻花）與模組二（葡萄）並列：兩地一東一西、
      兩種作物一觀賞一食用、兩種文化一和歌一詩酒。但兩者的氣候信號
      卻殊途同歸。這說明了氣候變遷作為「超物件」（Morton, 2013）的
      什麼特性？你如何向非科學背景的親友講述這兩個故事？
"""


if __name__ == "__main__":
    main()