# -*- coding: utf-8 -*-
"""
CCESD 通識模組一：從京都櫻花紀錄看氣候變遷
=============================================

本程式下載四個公開資料集，並在同一個時間軸上比較：

1. 京都櫻花盛開日 (Aono & Kazui, 2008; Aono & Saito, 2010)
   西元 812 年起，目前已知最長的連續植物物候紀錄。
2. NASA GISTEMP v4 全球地表溫度異常（1880 年起）
3. NOAA Mauna Loa 年均 CO₂ 濃度（1959 年起，直接大氣觀測）
4. Law Dome 冰芯 CO₂（MacFarling Meure et al. 2006；~AD 13 起，古氣候重建）

第 3 與第 4 資料串接後，可得到一條從西元初年到今日的完整 CO₂ 序列，
時間尺度與櫻花紀錄對齊，方便在同一張圖上比對前工業期的穩定基線與
工業化後的陡升。

執行：
    python kyoto_sakura_climate.py

相依套件：pandas, numpy, matplotlib, scipy, statsmodels
    pip install pandas numpy matplotlib scipy statsmodels

輸出：
    data/     原始下載檔（首次執行會下載，之後從本地讀取）
    output/   整併後的 CSV、四宮格圖、相關矩陣、迴歸結果

授權：程式碼 CC BY 4.0；資料遵循各來源原授權條款。
"""

from __future__ import annotations

import sys
import urllib.request
from io import StringIO
from pathlib import Path

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy import stats


# ----------------------------------------------------------------------------
# 0. 設定：路徑、中文字型、下載 header
# ----------------------------------------------------------------------------

DATA_DIR = Path("data")
OUT_DIR = Path("output")
DATA_DIR.mkdir(exist_ok=True)
OUT_DIR.mkdir(exist_ok=True)

# 有些政府與學術伺服器會拒絕沒有 User-Agent 的請求
HEADERS = {"User-Agent": "Mozilla/5.0 (educational research)"}


def setup_chinese_font() -> None:
    """在 Windows / macOS / Linux 自動挑選可用的中文字型。"""
    candidates = [
        "Microsoft JhengHei",       # Windows 正體
        "Microsoft YaHei",          # Windows 簡體
        "PingFang TC", "PingFang SC",  # macOS
        "Heiti TC", "Heiti SC",
        "Noto Sans CJK TC", "Noto Sans CJK SC",  # Linux
        "WenQuanYi Zen Hei",
        "SimHei",
    ]
    available = {f.name for f in mpl.font_manager.fontManager.ttflist}
    for name in candidates:
        if name in available:
            mpl.rcParams["font.sans-serif"] = [name]
            break
    else:
        print("[警告] 找不到中文字型，圖表中文可能顯示為方塊。")
    mpl.rcParams["axes.unicode_minus"] = False


def download(url: str, filename: str) -> Path:
    """下載到本地 data/ 目錄並快取；檔案已存在則直接回傳路徑。"""
    path = DATA_DIR / filename
    if path.exists() and path.stat().st_size > 0:
        return path
    print(f"下載中：{url}")
    req = urllib.request.Request(url, headers=HEADERS)
    with urllib.request.urlopen(req, timeout=60) as resp, open(path, "wb") as f:
        f.write(resp.read())
    print(f"  → 已儲存至 {path} ({path.stat().st_size:,} bytes)")
    return path


# ----------------------------------------------------------------------------
# 1. 資料來源一：京都櫻花盛開日
# ----------------------------------------------------------------------------
# Aono, Y. & Kazui, K. (2008). Phenological data series of cherry tree
#   flowering in Kyoto, Japan. Int. J. Climatol., 28, 905-914.
# Aono, Y. & Saito, S. (2010). Clarifying springtime temperature
#   reconstructions of the medieval period. Int. J. Biometeorol., 54, 211-219.
# 原始資料：http://atmenv.envi.osakafu-u.ac.jp/aono/kyophenotemp4/
# 以下使用 Ryo Nakagawara 整理後的 CSV 鏡像（欄位較易解析）

def load_sakura() -> pd.DataFrame:
    url = "https://raw.githubusercontent.com/Ryo-N7/sakura_bloom/master/Kyoto_Flowers.csv"
    path = download(url, "kyoto_flowers.csv")

    df = pd.read_csv(path)
    # 原始欄位：AD, Full-flowering date (DOY), Full-flowering date,
    #           Source code, Data type code, Reference Name
    df = df.rename(columns={
        "AD": "year",
        "Full-flowering date (DOY)": "doy",   # day-of-year，一年中第幾天
    })
    df = df[["year", "doy"]].copy()
    df["doy"] = pd.to_numeric(df["doy"], errors="coerce")
    df = df.dropna(subset=["doy"])
    df["doy"] = df["doy"].astype(int)
    df["year"] = df["year"].astype(int)
    return df.sort_values("year").reset_index(drop=True)


# ----------------------------------------------------------------------------
# 2. 資料來源二：NASA GISTEMP v4 全球地表溫度異常
# ----------------------------------------------------------------------------
# GISTEMP Team (2026). GISS Surface Temperature Analysis, version 4.
# https://data.giss.nasa.gov/gistemp/
# 基期：1951-1980；單位：°C

def load_giss() -> pd.DataFrame:
    url = "https://data.giss.nasa.gov/gistemp/tabledata_v4/GLB.Ts+dSST.csv"
    path = download(url, "giss_glb_ts_dsst.csv")

    # 檔頭第一行為說明文字，實際欄位從第二行起
    df = pd.read_csv(path, skiprows=1, na_values=["***", "****", "*****"])
    df = df.rename(columns={"Year": "year", "J-D": "temp_anomaly"})
    df["temp_anomaly"] = pd.to_numeric(df["temp_anomaly"], errors="coerce")
    df = df[["year", "temp_anomaly"]].dropna()

    # v4 CSV 版本已是 °C，但若偵測到整數型（0.01°C）則自動換算
    if df["temp_anomaly"].abs().max() > 10:
        df["temp_anomaly"] = df["temp_anomaly"] / 100.0
    df["year"] = df["year"].astype(int)
    return df.sort_values("year").reset_index(drop=True)


# ----------------------------------------------------------------------------
# 3. 資料來源三：NOAA Mauna Loa 年均 CO₂ 濃度
# ----------------------------------------------------------------------------
# Lan, X., Tans, P. & Thoning, K.W. (2026). Trends in globally-averaged CO2.
# NOAA Global Monitoring Laboratory. https://gml.noaa.gov/ccgg/trends/
# 單位：ppm（體積比百萬分之一）

def load_co2() -> pd.DataFrame:
    url = "https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_annmean_mlo.txt"
    path = download(url, "co2_annmean_mlo.txt")

    # 格式：以 # 為註解，欄位為 year, mean, unc（空白分隔）
    df = pd.read_csv(
        path, comment="#", sep=r"\s+", header=None,
        names=["year", "co2_ppm", "uncertainty"],
        engine="python",
    )
    df["year"] = df["year"].astype(int)
    df["co2_ppm"] = pd.to_numeric(df["co2_ppm"], errors="coerce")
    return df[["year", "co2_ppm"]].dropna().sort_values("year").reset_index(drop=True)


# ----------------------------------------------------------------------------
# 3b. 資料來源四：Law Dome 冰芯 CO₂（古氣候重建，~AD 13 起）
# ----------------------------------------------------------------------------
# MacFarling Meure, C. et al. (2006). Law Dome CO2, CH4 and N2O Ice Core
#   Records Extended to 2000 years BP. Geophys. Res. Lett., 33, L14810.
# DOI: 10.1029/2006GL026152
# 原始資料：NOAA NCEI Paleoclimatology
#   https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/law/
# 欄位：sample, age_CE, CO2, citation（tab 分隔）
# 同一年可能有多筆（DSS、DE08、DE08-2 三個冰芯 + DSSW20K 雪坑空氣 + Cape Grim 大氣），
# 本函式將同年多個樣品取算術平均作為年均值。

def load_law_dome() -> pd.DataFrame:
    url = ("https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/"
           "law/law2006-co2-noaa.txt")
    path = download(url, "law2006_co2_noaa.txt")

    # 檔頭有大量 '#' 註解行，之後才是資料
    # 資料前面有一行 tab 分隔的欄位標頭 "sample age_CE CO2 citation"
    with open(path, "r", encoding="utf-8", errors="replace") as f:
        lines = f.readlines()

    header_idx = None
    for i, line in enumerate(lines):
        stripped = line.lstrip()
        if not stripped.startswith("#") and stripped.startswith("sample"):
            header_idx = i
            break
    if header_idx is None:
        raise RuntimeError("找不到 Law Dome 資料欄位標頭 'sample'；檔案格式可能已變更")

    data_text = "".join(lines[header_idx:])
    df = pd.read_csv(StringIO(data_text), sep="\t", engine="python")

    df = df.rename(columns={"age_CE": "year_frac", "CO2": "co2_ppm"})
    df["year_frac"] = pd.to_numeric(df["year_frac"], errors="coerce")
    df["co2_ppm"] = pd.to_numeric(df["co2_ppm"], errors="coerce")
    df = df.dropna(subset=["year_frac", "co2_ppm"])

    # 年代取整；同年多樣品取算術平均
    df["year"] = df["year_frac"].astype(int)
    yearly = (df.groupby("year", as_index=False)["co2_ppm"]
                .mean()
                .sort_values("year")
                .reset_index(drop=True))
    return yearly


def load_co2_combined() -> pd.DataFrame:
    """合併 Law Dome 與 Mauna Loa，形成單一的長時段 CO₂ 序列。

    策略：以 1958 年為切點
    - 1958 年（含）之前：Law Dome 冰芯／雪坑空氣（古氣候重建）
    - 1959 年起：Mauna Loa 直接大氣觀測（精度更高）
    """
    law = load_law_dome()
    mlo = load_co2()
    law_pre = law[law["year"] <= 1958].copy()
    law_pre["source"] = "Law Dome"
    mlo_post = mlo.copy()
    mlo_post["source"] = "Mauna Loa"
    combined = pd.concat([law_pre, mlo_post], ignore_index=True)
    combined = combined.sort_values("year").reset_index(drop=True)
    return combined


# ----------------------------------------------------------------------------
# 4. 整併三個資料集
# ----------------------------------------------------------------------------

def build_master() -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame,
                             pd.DataFrame, pd.DataFrame]:
    sakura = load_sakura()
    giss = load_giss()
    co2_combined = load_co2_combined()  # Law Dome + Mauna Loa 接合
    mlo = load_co2()                     # 單獨保留 Mauna Loa 供現代時段用

    # 以櫻花年份為主鍵（時間跨度最長），左併溫度與合併後的 CO₂
    master = (sakura
              .merge(giss, on="year", how="left")
              .merge(co2_combined[["year", "co2_ppm"]], on="year", how="left"))

    out_path = OUT_DIR / "master_merged.csv"
    master.to_csv(out_path, index=False, encoding="utf-8-sig")
    co2_combined.to_csv(OUT_DIR / "co2_combined.csv",
                        index=False, encoding="utf-8-sig")
    print(f"\n[整併結果] 共 {len(master)} 筆（西元 {master['year'].min()}–{master['year'].max()}）")
    print(f"          已輸出 → {out_path}")
    print(f"          合併後 CO₂ 序列 → {OUT_DIR/'co2_combined.csv'}  "
          f"({co2_combined['year'].min()}–{co2_combined['year'].max()})")
    return master, sakura, giss, co2_combined, mlo


# ----------------------------------------------------------------------------
# 5. 視覺化：四宮格圖
# ----------------------------------------------------------------------------

def plot_four_panel(master: pd.DataFrame, sakura: pd.DataFrame,
                    giss: pd.DataFrame, co2_combined: pd.DataFrame,
                    mlo: pd.DataFrame) -> None:
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))

    # (a) 櫻花盛開日全序列（812–present）
    ax = axes[0, 0]
    ax.scatter(sakura["year"], sakura["doy"], s=6, alpha=0.45, color="#d48fa6",
               label=f"年度觀測（n={len(sakura)}）")
    roll = (sakura.set_index("year")["doy"]
            .rolling(30, min_periods=10).mean())
    ax.plot(roll.index, roll.values, color="#8b2a4a", lw=1.8,
            label="30 年移動平均")
    ax.invert_yaxis()  # 反轉 Y 軸：上方 = 開得早
    ax.set_xlabel("西元年")
    ax.set_ylabel("盛開日（一年中第幾天）\n↑ 越上越早")
    ax.set_title("(a) 京都櫻花盛開日 812–present\n"
                 "資料：Aono & Kazui (2008); Aono & Saito (2010)", fontsize=10)
    ax.legend(loc="lower left", fontsize=8)
    ax.grid(alpha=0.3)

    # (b) 全球溫度異常
    ax = axes[0, 1]
    ax.plot(giss["year"], giss["temp_anomaly"], color="#7f1d1d", lw=1)
    ax.fill_between(giss["year"], 0, giss["temp_anomaly"],
                    where=giss["temp_anomaly"] > 0, alpha=0.35, color="#c0392b")
    ax.fill_between(giss["year"], 0, giss["temp_anomaly"],
                    where=giss["temp_anomaly"] <= 0, alpha=0.35, color="#1f618d")
    ax.axhline(0, color="k", lw=0.6)
    ax.set_xlabel("西元年")
    ax.set_ylabel("溫度異常（°C，相對 1951–1980）")
    ax.set_title("(b) 全球地表溫度異常\n資料：NASA GISTEMP v4", fontsize=10)
    ax.grid(alpha=0.3)

    # (c) CO₂ 濃度 ─ Law Dome（冰芯）+ Mauna Loa（大氣直接觀測）
    ax = axes[1, 0]
    law = co2_combined[co2_combined["source"] == "Law Dome"]
    ax.scatter(law["year"], law["co2_ppm"], s=10, color="#5b7891",
               alpha=0.8, label="Law Dome 冰芯 (AD 13–1958)")
    ax.plot(mlo["year"], mlo["co2_ppm"], color="#c0392b", lw=2.0,
            label="Mauna Loa 大氣觀測 (1959–present)")
    ax.axhline(280, color="gray", ls="--", lw=0.8,
               label="工業化前 ~280 ppm 基線")
    ax.set_xlabel("西元年")
    ax.set_ylabel("CO₂ 濃度（ppm）")
    ax.set_title("(c) 大氣 CO₂ 濃度 ─ 古氣候重建接合現代觀測\n"
                 "資料：MacFarling Meure et al. (2006); NOAA Mauna Loa", fontsize=10)
    ax.legend(loc="upper left", fontsize=8)
    ax.grid(alpha=0.3)

    # (d) 共同觀測時段：三變數標準化疊圖
    ax = axes[1, 1]
    common = master.dropna(subset=["doy", "temp_anomaly", "co2_ppm"])
    if len(common) >= 10:
        def z(x):
            return (x - x.mean()) / x.std()

        # 為了直覺比較：盛開日取負值，使三條線皆「上升 = 暖化訊號強」
        ax.plot(common["year"], -z(common["doy"]),
                color="#d48fa6", lw=1.5, label="盛開日（反轉後，越上越早）")
        ax.plot(common["year"], z(common["temp_anomaly"]),
                color="#c0392b", lw=1.5, label="全球溫度異常")
        ax.plot(common["year"], z(common["co2_ppm"]),
                color="#2c3e50", lw=1.5, label="CO₂ 濃度")
        ax.axhline(0, color="k", lw=0.5)
        ax.set_xlabel("西元年")
        ax.set_ylabel("標準化 z-score")
        ax.set_title(f"(d) 共同時段三變數比較\n"
                     f"n = {len(common)} 年（{common['year'].min()}–{common['year'].max()}）",
                     fontsize=10)
        ax.legend(loc="upper left", fontsize=8)
        ax.grid(alpha=0.3)

    plt.tight_layout()
    path = OUT_DIR / "figure_four_panel.png"
    plt.savefig(path, dpi=150, bbox_inches="tight")
    plt.close()
    print(f"[視覺化] 四宮格圖 → {path}")


def plot_sakura_vs_co2_long(sakura: pd.DataFrame,
                             co2_combined: pd.DataFrame) -> None:
    """長時段版本：櫻花 + 合併 CO₂ 同軸疊圖，讓學生直觀看到
    「前工業期 CO₂ 穩定、櫻花盛開日也穩定；兩者在 20 世紀同步轉折」。
    """
    fig, ax1 = plt.subplots(figsize=(12, 6))

    # 櫻花（30 年平均）
    sakura_roll = (sakura.set_index("year")["doy"]
                   .rolling(30, min_periods=10).mean())
    ax1.plot(sakura_roll.index, sakura_roll.values,
             color="#8b2a4a", lw=1.8, label="櫻花盛開日 30 年平均")
    ax1.set_xlabel("西元年")
    ax1.set_ylabel("櫻花盛開日（一年中第幾天）", color="#8b2a4a")
    ax1.tick_params(axis="y", labelcolor="#8b2a4a")
    ax1.invert_yaxis()
    ax1.grid(alpha=0.3)

    # CO₂（第二 Y 軸）
    ax2 = ax1.twinx()
    law = co2_combined[co2_combined["source"] == "Law Dome"]
    mlo_part = co2_combined[co2_combined["source"] == "Mauna Loa"]
    ax2.scatter(law["year"], law["co2_ppm"], s=12, color="#5b7891",
                alpha=0.75, label="Law Dome 冰芯 CO₂")
    ax2.plot(mlo_part["year"], mlo_part["co2_ppm"], color="#c0392b", lw=2.0,
             label="Mauna Loa CO₂")
    ax2.set_ylabel("CO₂ 濃度（ppm）", color="#2c3e50")
    ax2.tick_params(axis="y", labelcolor="#2c3e50")

    # 合併圖例
    lines1, labels1 = ax1.get_legend_handles_labels()
    lines2, labels2 = ax2.get_legend_handles_labels()
    ax1.legend(lines1 + lines2, labels1 + labels2,
               loc="lower left", fontsize=9)

    ax1.set_title("京都櫻花盛開日 vs 大氣 CO₂ 濃度（西元 812 年至今）",
                  fontsize=12)
    plt.tight_layout()
    path = OUT_DIR / "figure_sakura_co2_1200yr.png"
    plt.savefig(path, dpi=150, bbox_inches="tight")
    plt.close()
    print(f"[視覺化] 千年尺度疊圖 → {path}")


def plot_scatter_doy_vs_temp(master: pd.DataFrame) -> None:
    """散佈圖：盛開日 vs 溫度異常，附迴歸線。"""
    sub = master.dropna(subset=["doy", "temp_anomaly"])
    if len(sub) < 10:
        return
    fig, ax = plt.subplots(figsize=(8, 6))
    ax.scatter(sub["temp_anomaly"], sub["doy"], alpha=0.55, s=25,
               c=sub["year"], cmap="viridis")
    cb = plt.colorbar(ax.collections[0], ax=ax)
    cb.set_label("西元年")

    # 加迴歸線
    x = sub["temp_anomaly"].values
    y = sub["doy"].values
    slope, intercept, r, p, se = stats.linregress(x, y)
    xs = np.linspace(x.min(), x.max(), 100)
    ax.plot(xs, intercept + slope * xs, "r-", lw=2,
            label=f"y = {intercept:.2f} + {slope:.2f} x\n"
                  f"r = {r:+.3f}  p = {p:.2e}  n = {len(sub)}")
    ax.invert_yaxis()
    ax.set_xlabel("全球溫度異常（°C）")
    ax.set_ylabel("京都櫻花盛開日（DOY）\n↑ 越上越早")
    ax.set_title("盛開日與全球溫度異常的關係")
    ax.legend(loc="upper right")
    ax.grid(alpha=0.3)

    path = OUT_DIR / "figure_scatter_doy_temp.png"
    plt.tight_layout()
    plt.savefig(path, dpi=150, bbox_inches="tight")
    plt.close()
    print(f"[視覺化] 散佈圖 → {path}")


# ----------------------------------------------------------------------------
# 6. 進階分析：相關 + 線性迴歸 + 多元迴歸
# ----------------------------------------------------------------------------

def correlation_analysis(master: pd.DataFrame) -> pd.DataFrame:
    """Pearson 與 Spearman 相關，以兩者交叉驗證線性/單調關係。"""
    sub = master.dropna(subset=["doy", "temp_anomaly", "co2_ppm"])
    print("\n" + "=" * 64)
    print("相關分析（共同觀測時段）")
    print("=" * 64)
    print(f"樣本數 n = {len(sub)}  "
          f"（{sub['year'].min()}–{sub['year'].max()}）\n")

    rows = []
    pairs = [
        ("doy", "temp_anomaly", "盛開日 vs 全球溫度異常"),
        ("doy", "co2_ppm",       "盛開日 vs CO₂ 濃度"),
        ("temp_anomaly", "co2_ppm", "溫度異常 vs CO₂ 濃度"),
    ]
    for a, b, label in pairs:
        pear = stats.pearsonr(sub[a], sub[b])
        spear = stats.spearmanr(sub[a], sub[b])
        rows.append({
            "pair": label,
            "pearson_r": round(pear.statistic, 4),
            "pearson_p": f"{pear.pvalue:.2e}",
            "spearman_rho": round(spear.statistic, 4),
            "spearman_p": f"{spear.pvalue:.2e}",
            "n": len(sub),
        })
        print(f"{label:<25s}  Pearson r = {pear.statistic:+.3f}  (p = {pear.pvalue:.2e})")
        print(f"{'':25s}  Spearman ρ = {spear.statistic:+.3f}  (p = {spear.pvalue:.2e})")

    corr_df = pd.DataFrame(rows)
    corr_df.to_csv(OUT_DIR / "correlation_table.csv", index=False, encoding="utf-8-sig")
    return sub


def regression_analysis(sub: pd.DataFrame) -> None:
    """簡單與多元 OLS 迴歸；係數、標準誤、R²、Durbin-Watson 一併輸出。"""
    # (A) 簡單線性：盛開日 ~ 溫度異常
    print("\n" + "=" * 64)
    print("模型 1：盛開日 = β₀ + β₁ × 溫度異常")
    print("=" * 64)
    X1 = sm.add_constant(sub["temp_anomaly"])
    m1 = sm.OLS(sub["doy"], X1).fit()
    print(m1.summary())

    # (B) 多元：盛開日 ~ 溫度異常 + CO₂
    print("\n" + "=" * 64)
    print("模型 2：盛開日 = β₀ + β₁ × 溫度異常 + β₂ × CO₂")
    print("=" * 64)
    X2 = sm.add_constant(sub[["temp_anomaly", "co2_ppm"]])
    m2 = sm.OLS(sub["doy"], X2).fit()
    print(m2.summary())

    # 將兩個模型的摘要輸出到文字檔
    with open(OUT_DIR / "regression_summary.txt", "w", encoding="utf-8") as f:
        f.write("模型 1：盛開日 ~ 溫度異常\n")
        f.write("=" * 64 + "\n")
        f.write(str(m1.summary()) + "\n\n")
        f.write("模型 2：盛開日 ~ 溫度異常 + CO₂\n")
        f.write("=" * 64 + "\n")
        f.write(str(m2.summary()) + "\n")
    print(f"\n[迴歸結果] 已輸出 → {OUT_DIR / 'regression_summary.txt'}")

    # 教學重點提示
    print("\n" + "-" * 64)
    print("【課堂討論要點】")
    print("-" * 64)
    b1 = m1.params["temp_anomaly"]
    print(f"1. 模型 1 係數：溫度每上升 1°C，盛開日提前約 {abs(b1):.1f} 天")
    print(f"   （β₁ = {b1:+.2f}，p = {m1.pvalues['temp_anomaly']:.2e}，R² = {m1.rsquared:.3f}）")
    print(f"2. 模型 2 加入 CO₂ 後 R² 從 {m1.rsquared:.3f} 變為 {m2.rsquared:.3f}")
    if m2.rsquared_adj < m1.rsquared_adj + 0.01:
        print("   → 調整 R² 幾乎沒變化：溫度與 CO₂ 高度共線，CO₂ 的邊際貢獻有限")
    print("3. 請學生思考：")
    print("   (a) 盛開日是氣候『結果』，CO₂ 是氣候『驅動力』，中間的機制是什麼？")
    print("   (b) Pearson r 與 Spearman ρ 差距大不大？代表關係是否線性？")
    print("   (c) 1940–1970 年全球溫度曾短暫停滯，櫻花資料有沒有跟著反映？")


# ----------------------------------------------------------------------------
# 7. 主流程
# ----------------------------------------------------------------------------

def main() -> int:
    setup_chinese_font()

    try:
        master, sakura, giss, co2_combined, mlo = build_master()
    except Exception as e:
        print(f"[錯誤] 資料下載或整併失敗：{e}", file=sys.stderr)
        print("請檢查網路連線或來源網址是否變動。", file=sys.stderr)
        return 1

    plot_four_panel(master, sakura, giss, co2_combined, mlo)
    plot_sakura_vs_co2_long(sakura, co2_combined)
    plot_scatter_doy_vs_temp(master)

    sub = correlation_analysis(master)
    if len(sub) >= 30:
        regression_analysis(sub)
    else:
        print("[警告] 共同時段樣本少於 30 筆，跳過迴歸分析")

    print("\n完成。所有輸出在 output/ 目錄下。")
    return 0


if __name__ == "__main__":
    sys.exit(main())