# -*- coding: utf-8 -*- """ =========================================================================== 附錄二｜完整程式碼（Appendix 2: Full Python Code） 從櫻花早開到夏季缺電：以「剪刀效應」與「死亡螺旋」視覺化 與真實數據之氣候變遷通識教學設計 =========================================================================== 本附錄整合四個跨領域氣候數據視覺化教學模組的完整 Python 程式碼。 所有資料皆自官方來源即時下載；每次執行都會反映各資料庫的最新版本。 【模組目錄】 模組一｜京都櫻花剪刀效應（The Scissors Effect） 資料：Aono & Kazui (2008); NASA GISTEMP v4 模組二｜勃艮第葡萄歷史斷層（The Breakpoint） 資料：Chuine et al. (2004), NOAA 古氣候庫; NASA GISTEMP v4 模組四｜北極海冰死亡螺旋（The Death Spiral） 資料：NSIDC Sea Ice Index v4 模組五｜台灣電力在地剪刀（The Local Scissors） 資料：台灣電力公司 / 政府資料開放平臺 【使用方式】 # Colab 首次執行 !apt-get install -y fonts-noto-cjk -q # 執行全部四個模組 python appendix_full_code.py # 執行單一模組 python appendix_full_code.py 1 # 模組一 python appendix_full_code.py 2 # 模組二 python appendix_full_code.py 4 # 模組四 python appendix_full_code.py 5 # 模組五 【運算環境】 Python 3.10+；相依套件：pandas, numpy, matplotlib, statsmodels, scipy 【授權】 本程式碼依 CC BY 4.0 條款釋出；資料源各有其授權條款，請參見各模組 說明與原始資料發布方。 =========================================================================== """ import io import sys 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.collections import LineCollection from matplotlib.colors import Normalize, ListedColormap, BoundaryNorm from matplotlib.cm import ScalarMappable from matplotlib.ticker import MultipleLocator from scipy import stats import statsmodels.api as sm warnings.filterwarnings("ignore") # ============================================================================ # PART 0｜共用基礎設施 # ============================================================================ # -- 中文字型設定（Colab 需先安裝 fonts-noto-cjk，執行後重啟 Runtime） -- 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 def load_gistemp() -> pd.DataFrame: """ 共用函式｜NASA GISTEMP v4 陸地+海洋年均溫異常值。 基期：1951–1980；單位：°C；範圍：1880 年迄今。 模組一、模組二皆呼叫此函式。 """ 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"[GISTEMP] 共讀入 {len(df):,} 筆年度紀錄；" f"範圍 {df['AD'].min()}–{df['AD'].max()}") return df.reset_index(drop=True) 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)}") def _banner(title: str, char: str = "=", width: int = 72) -> None: """共用函式｜印出視覺化分隔標題。""" print("\n" + char * width) print(f" {title}") print(char * width) # ============================================================================ # PART 1｜模組一：京都櫻花剪刀效應 # ============================================================================ # 【教學目標】以 Python 讀取、清理並合併兩個跨度不同的時序資料； # 解讀 OLS 迴歸與 Pearson/Spearman 相關；繪製雙軸剪刀圖。 # 【資料】Aono & Kazui (2008), Int. J. Climatol., 28, 905–914 # 原站點 osakafu-u.ac.jp 於 2026 年初青野靖之教授辭世後失效， # 現改用 Our World in Data 與 GitHub 鏡像雙重備援。 # ============================================================================ def m1_load_kyoto_cherry() -> pd.DataFrame: """讀取京都御所花期資料（1200+ 年）。""" errors = [] # 來源 1：OWID（含最新年份） owid_url = ("https://ourworldindata.org/grapher/" "date-of-the-peak-cherry-tree-blossom-in-kyoto.csv") try: df = pd.read_csv(owid_url) value_col = [c for c in df.columns if c not in ("Entity", "Code", "Year")][0] df = df[["Year", value_col]].rename( columns={"Year": "AD", value_col: "FullFloweringDOY"}) df = df.dropna().copy() df["AD"] = df["AD"].astype(int) df["FullFloweringDOY"] = df["FullFloweringDOY"].astype(float) print(f"[京都櫻花] 來源：OWID；共 {len(df):,} 筆；" f"{df['AD'].min()}–{df['AD'].max()}") return df.reset_index(drop=True) except Exception as e: errors.append(f"OWID: {e}") # 來源 2：GitHub 鏡像 github_url = ("https://raw.githubusercontent.com/Ryo-N7/" "sakura_bloom/master/Kyoto_Flowers.csv") try: df = pd.read_csv(github_url) doy_col = [c for c in df.columns if "DOY" in c][0] df = df[["AD", doy_col]].rename(columns={doy_col: "FullFloweringDOY"}) df["FullFloweringDOY"] = pd.to_numeric( df["FullFloweringDOY"], errors="coerce") df = df.dropna().copy() df["AD"] = df["AD"].astype(int) print(f"[京都櫻花] 來源：GitHub 鏡像；共 {len(df):,} 筆；" f"{df['AD'].min()}–{df['AD'].max()}") return df.reset_index(drop=True) except Exception as e: errors.append(f"GitHub: {e}") raise RuntimeError("兩個備援來源皆失敗：\n " + "\n ".join(errors)) def m1_correlation_analysis(df: pd.DataFrame) -> dict: """Pearson + Spearman 雙軌相關分析。""" x = df["TempAnomaly"].values y = df["FullFloweringDOY"].values pearson_r, pearson_p = stats.pearsonr(x, y) spearman_r, spearman_p = stats.spearmanr(x, y) print(f"\n[相關性] n = {len(df)}") print(f" Pearson r = {pearson_r:+.3f} (p = {pearson_p:.2e})") print(f" Spearman ρ = {spearman_r:+.3f} (p = {spearman_p:.2e})") return {"n": len(df), "pearson_r": pearson_r, "pearson_p": pearson_p, "spearman_r": spearman_r, "spearman_p": spearman_p} def m1_ols_regression(df: pd.DataFrame): """OLS 迴歸含殘差常態性檢定。""" X = sm.add_constant(df["TempAnomaly"].values) y = df["FullFloweringDOY"].values model = sm.OLS(y, X).fit() beta = model.params[1] ci = model.conf_int()[1] print(f"\n[OLS] β = {beta:+.3f} 天/°C 95% CI [{ci[0]:+.3f}, {ci[1]:+.3f}]") print(f" R² = {model.rsquared:.3f} p = {model.pvalues[1]:.2e}") resid = model.resid jb_stat, jb_p, _, _ = sm.stats.stattools.jarque_bera(resid) print(f"[殘差] Jarque-Bera p = {jb_p:.3f}") return model def m1_plot_millennium_cherry(cherry: pd.DataFrame, savepath: str = "fig1a_millennium_cherry.png"): """圖 1a：千年花期序列。""" fig, ax = plt.subplots(figsize=(11, 4.5)) ax.scatter(cherry["AD"], cherry["FullFloweringDOY"], s=6, c="#C71585", alpha=0.5, label="滿開日（DOY）") c = cherry.sort_values("AD").copy() c["Roll30"] = c["FullFloweringDOY"].rolling(30, min_periods=15).median() ax.plot(c["AD"], c["Roll30"], color="#800020", lw=2.0, label="30 年滾動中位數") ax.set_xlabel("西元年（AD）") ax.set_ylabel("滿開日年積日（Day of Year）") ax.set_title("圖 1a｜京都御所櫻花滿開日千年序列（Aono & Kazui, 2008）") ax.invert_yaxis() ax.legend(loc="lower left", frameon=False) ax.grid(True, alpha=0.25) fig.tight_layout() fig.savefig(savepath, dpi=200, bbox_inches="tight") print(f"[已輸出] {savepath}") plt.close(fig) def m1_plot_scissors_effect(merged: pd.DataFrame, model, savepath: str = "fig1b_scissors_effect.png"): """圖 1b：剪刀效應雙軸圖（主圖）。""" fig, ax1 = plt.subplots(figsize=(11, 5.5)) col_temp = "#B22222" ax1.plot(merged["AD"], merged["TempAnomaly"], color=col_temp, lw=1.1, alpha=0.55, label="年均溫異常") ax1.plot(merged["AD"], merged["TempAnomaly"].rolling(10, center=True, min_periods=5).mean(), color=col_temp, lw=2.8, label="10 年移動平均（溫度）") ax1.set_xlabel("西元年（AD）") 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_cherry = "#E75480" ax2.plot(merged["AD"], merged["FullFloweringDOY"], color=col_cherry, lw=1.1, alpha=0.55, label="京都滿開日") ax2.plot(merged["AD"], merged["FullFloweringDOY"].rolling(10, center=True, min_periods=5).mean(), color=col_cherry, lw=2.8, label="10 年移動平均（花期）") ax2.set_ylabel("京都櫻花滿開日（DOY）", color=col_cherry) ax2.tick_params(axis="y", labelcolor=col_cherry) ax2.invert_yaxis() beta = model.params[1]; r2 = model.rsquared; p = model.pvalues[1] sig = "p < .001" if p < 1e-3 else f"p = {p:.3f}" ax1.set_title( f"圖 1b｜剪刀效應：氣溫 ↑ 與花期 ↓ 之耦合\n" f"OLS：β = {beta:+.2f} 天/°C，R² = {r2:.2f}，{sig}，n = {len(merged)}", pad=12, ) 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.xaxis.set_major_locator(MultipleLocator(20)) ax1.grid(True, alpha=0.2) fig.tight_layout() fig.savefig(savepath, dpi=200, bbox_inches="tight") print(f"[已輸出] {savepath}") plt.close(fig) def run_module_1() -> None: """模組一主流程。""" _banner("模組一｜剪刀效應：京都櫻花 × 全球氣溫") cherry = m1_load_kyoto_cherry() temp = load_gistemp() merged = cherry.merge(temp, on="AD", how="inner") print(f"[合併] 重疊 {len(merged):,} 筆；" f"{merged['AD'].min()}–{merged['AD'].max()}") m1_correlation_analysis(merged) model = m1_ols_regression(merged) m1_plot_millennium_cherry(cherry) m1_plot_scissors_effect(merged, model) # ============================================================================ # PART 2｜模組二：勃艮第葡萄歷史斷層 # ============================================================================ # 【教學目標】讀取 NOAA 古氣候純文字資料；以 t 檢定與 CUSUM 偵測斷點； # 區分「漸進趨勢」與「斷層變化」兩種統計結構。 # 【資料】Chuine et al. (2004), Nature, 432, 289–290 # DOI: 10.25921/e3mc-wm08；範圍 1370–2003； # 原檔欄位：Year + 9 月 1 日起算天數（轉為 DOY） # 【延伸】Labbé et al. (2019), Clim. Past, 15, 1485–1501；指出 1988 斷層 # ============================================================================ def m2_load_burgundy_ghd() -> pd.DataFrame: """讀取 NOAA 存檔之勃艮第葡萄採收日（1370–2003）。""" url = ("https://www.ncei.noaa.gov/pub/data/paleo/historical/" "france/burgundy2004.txt") with urllib.request.urlopen(url, timeout=30) as resp: text = resp.read().decode("latin-1") # 尋找資料表頭 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 = float(parts[1]) rows.append((year, days)) except ValueError: continue df = pd.DataFrame(rows, columns=["AD", "DaysAfterSep1"]) df["HarvestDOY"] = 244 + df["DaysAfterSep1"] # 9/1 ≈ DOY 244 print(f"[勃艮第葡萄] 共 {len(df):,} 筆；" f"{df['AD'].min()}–{df['AD'].max()}") return df.reset_index(drop=True) def m2_breakpoint_ttest(df: pd.DataFrame, break_year: int = 1988, value_col: str = "HarvestDOY") -> dict: """先驗斷點 Welch's t 檢定 + Cohen's d。""" 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() 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 檢定 @ {break_year}]") print(f" 前段 n={len(before)} 均值={before.mean():.2f} SD={before.std(ddof=1):.2f}") print(f" 後段 n={len(after)} 均值={after.mean():.2f} SD={after.std(ddof=1):.2f}") print(f" 提前 {shift_days:.2f} 天；t = {t_stat:.3f}，p = {p_val:.2e}") print(f" Cohen's d = {cohens_d:.2f}") 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} def m2_cusum_scan(df: pd.DataFrame, value_col: str = "HarvestDOY", min_segment: int = 30) -> dict: """CUSUM 暴力搜索最佳斷點（不預設位置）。""" years = df["AD"].values y = df[value_col].values n = len(y) best_score, best_idx = np.inf, None for i in range(min_segment, n - min_segment): b, a = y[:i], y[i:] score = b.var(ddof=1) * (len(b) - 1) + a.var(ddof=1) * (len(a) - 1) if score < best_score: best_score, best_idx = score, i best_year = int(years[best_idx]) print(f"\n[CUSUM 掃描] 最佳斷點 {best_year}；" f"前均 {y[:best_idx].mean():.2f}，後均 {y[best_idx:].mean():.2f}") return {"best_year": best_year, "mean_before": y[:best_idx].mean(), "mean_after": y[best_idx:].mean()} def m2_temperature_coupling(ghd: pd.DataFrame, temp: pd.DataFrame): """GHD × GISTEMP 現代重疊期相關與 OLS。""" merged = ghd.merge(temp, on="AD", how="inner") print(f"\n[氣溫耦合] {merged['AD'].min()}–{merged['AD'].max()}; n = {len(merged)}") r, p = stats.pearsonr(merged["TempAnomaly"], merged["HarvestDOY"]) print(f" Pearson r = {r:+.3f} (p = {p:.2e})") X = sm.add_constant(merged["TempAnomaly"].values) model = sm.OLS(merged["HarvestDOY"].values, X).fit() print(f" OLS β = {model.params[1]:+.3f} 天/°C；R² = {model.rsquared:.3f}") return model def m2_plot_breakpoint(df: pd.DataFrame, tt: dict, savepath: str = "fig2a_breakpoint_panorama.png"): """圖 2a：六百年全景斷層圖（主圖）。""" fig, ax = plt.subplots(figsize=(11.5, 5.5)) ax.scatter(df["AD"], df["HarvestDOY"], s=7, c="#722F37", alpha=0.4, label="年採收日（DOY）") roll = df["HarvestDOY"].rolling(30, min_periods=15).mean() ax.plot(df["AD"], roll, color="#2C1810", lw=2.2, label="30 年移動平均") br = tt["break_year"] xlim = (df["AD"].min() - 5, df["AD"].max() + 5) y_pre, s_pre = tt["mean_before"], tt["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, s_post = tt["mean_after"], tt["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, label=f"{br} 斷層") ax.annotate( f"↓ {tt['shift_days']:.1f} 天\n(p = {tt['p_value']:.1e},\n" f" Cohen's d = {tt['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() 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) def m2_plot_distribution_comparison(df: pd.DataFrame, break_year: int = 1988, savepath: str = "fig2b_distribution_shift.png"): """圖 2b：斷點前後分布小提琴圖。""" 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) def run_module_2() -> None: """模組二主流程。""" _banner("模組二｜歷史斷層：勃艮第葡萄 × 全球氣溫") ghd = m2_load_burgundy_ghd() temp = load_gistemp() tt = m2_breakpoint_ttest(ghd, break_year=1988) cs = m2_cusum_scan(ghd) if abs(cs["best_year"] - 1988) <= 5: print(f"\n★ CUSUM 自動找到 {cs['best_year']}，與 1988 相距" f" {abs(cs['best_year'] - 1988)} 年；資料支持 1988 斷層。") m2_temperature_coupling(ghd, temp) m2_plot_breakpoint(ghd, tt) m2_plot_distribution_comparison(ghd) # ============================================================================ # PART 3｜模組四：北極海冰死亡螺旋 # ============================================================================ # 【教學目標】下載並合併 12 個月份 CSV；極坐標轉換；製作具時間色彩 # 映射的螺旋圖；線性外推「無冰北極」年份。 # 【資料】NSIDC Sea Ice Index v4，DOI: 10.7265/a98x-0f50 # 北半球月均海冰範圍，1979 年起，每月更新 # 【視覺靈感】Hawkins (2016), Pluck (2017) # ============================================================================ def m4_load_nsidc_monthly() -> pd.DataFrame: """下載並合併 NSIDC 12 個月份海冰範圍 CSV。""" 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) df.columns = [c.strip().replace("-", "_").lower() for c in df.columns] df = df[["year", "mo", "extent"]].copy() df["extent"] = pd.to_numeric(df["extent"], errors="coerce") df = df[df["extent"] > 0] 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 def m4_seasonal_cycle(df: pd.DataFrame) -> pd.DataFrame: """各月海冰面積的長期平均。""" s = df.groupby("mo")["extent"].agg(mean="mean", std="std", min_val="min", max_val="max") print("\n[季節循環] 各月海冰面積（百萬 km²）") print(s.round(2).to_string()) return s def m4_prepare_polar(df: pd.DataFrame) -> pd.DataFrame: """極坐標轉換：theta = 2π × (mo - 1) / 12；r = extent。""" df = df.copy() df["theta"] = 2 * np.pi * (df["mo"] - 1) / 12 df["r"] = df["extent"] return df def m4_plot_death_spiral(df: pd.DataFrame, savepath: str = "fig4b_death_spiral.png", cmap_name: str = "viridis"): """圖 4b：死亡螺旋主視覺（Hawkins–Pluck 風格）。""" df = m4_prepare_polar(df).sort_values(["year", "mo"]) year_min, year_max = df["year"].min(), df["year"].max() norm = Normalize(vmin=year_min, vmax=year_max) cmap = mpl.colormaps[cmap_name] fig = plt.figure(figsize=(10, 10)) ax = fig.add_subplot(111, projection="polar") ax.set_theta_zero_location("N") # 1 月朝上 ax.set_theta_direction(-1) # 順時針 segments, colors = [], [] arr = df[["theta", "r", "year", "mo"]].values for i in range(len(arr) - 1): t0, r0, y0, m0 = arr[i] t1, r1, _, _ = arr[i + 1] 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) ax.set_rmin(df["r"].min() * 0.9) ax.set_rmax(df["r"].max() * 1.05) ax.set_rticks([4, 8, 12, 16]) 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) cbar.set_label("年份", fontsize=10) 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) def m4_september_trend(df: pd.DataFrame) -> dict: """九月極小值 OLS 趨勢 + 無冰北極線性外推。""" sept = df[df["mo"] == 9].copy().sort_values("year").reset_index(drop=True) X = sm.add_constant(sept["year"].values) model = sm.OLS(sept["extent"].values, X).fit() slope, intercept = model.params[1], model.params[0] year_ice_free = (1.0 - intercept) / slope print(f"\n[九月趨勢] 每十年 {slope*10:+.3f} 百萬 km²；" f"R² = {model.rsquared:.3f}") print(f"[線性外推] 無冰門檻 (<1 M km²) 預期年份 ≈ {year_ice_free:.0f}") print(" 警語：純線性外推；實際系統具非線性反饋；IPCC AR6 推估 2050 前。") return {"sept_df": sept, "slope": slope, "intercept": intercept, "r_squared": model.rsquared, "per_decade": slope * 10, "year_ice_free_linear": year_ice_free, "model": model} def m4_plot_september_trend(res: dict, savepath: str = "fig4c_september_trend.png"): """圖 4c：九月極小值散點 + 迴歸線 + 外推。""" sept = res["sept_df"] slope, intercept = res["slope"], res["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 迴歸線") x_end = min(res["year_ice_free_linear"], 2060) x_ext = np.array([sept["year"].max(), x_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), ) 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"每十年 {res['per_decade']:.3f} M km²，R² = {res['r_squared']:.2f}") 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) def run_module_4() -> None: """模組四主流程。""" _banner("模組四｜死亡螺旋：北極海冰的時間軌跡") df = m4_load_nsidc_monthly() m4_seasonal_cycle(df) m4_plot_death_spiral(df) res = m4_september_trend(df) m4_plot_september_trend(res) # ============================================================================ # PART 4｜模組五：台灣電力在地剪刀 # ============================================================================ # 【教學目標】處理政府開放資料（Big5/UTF-8）；民國年轉西元；計算 # 台電燈號警戒日；以熱地圖與雙軸圖視覺化。 # 【資料】台灣電力公司 / 政府資料開放平臺： # #8307（歷年尖峰負載及備用容量率，自 1982 年起） # #24945（近三年每日尖峰備轉容量率） # 【燈號】綠 ≥10% | 黃 6–10% | 橘 3–6% | 紅 <3% # ============================================================================ def m5_download_taipower_csv(url: str, encodings: tuple = ("utf-8-sig", "utf-8", "big5", "cp950")) -> pd.DataFrame: """穩健台電下載器：多編碼備援 + 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}）") last_err = None for enc in encodings: try: return pd.read_csv(io.StringIO(raw.decode(enc))) except (UnicodeDecodeError, pd.errors.ParserError) as e: last_err = e continue raise RuntimeError(f"解碼失敗（嘗試：{encodings}）；最後錯誤：{last_err}") def m5_load_annual_reserve() -> pd.DataFrame: """讀取 1982 年起之年度尖峰負載與備用容量率（政府開放資料 #8307）。""" url = ("https://service.taipower.com.tw/data/opendata/" "apply/file/d003002/001.csv") df = m5_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"] for c in out.columns: out[c] = pd.to_numeric(out[c], 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 m5_load_daily_operating_reserve() -> pd.DataFrame: """讀取近三年每日備轉容量率（政府開放資料 #24945）。""" url = ("https://service.taipower.com.tw/data/opendata/" "apply/file/d006004/001.csv") df = m5_download_taipower_csv(url) col_date = _find_col(df, ["日期", "date"]) col_cap = _find_col(df, ["備轉容量(萬", "備轉容量 (萬"]) 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 = 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 def m5_long_term_scissors_regression(annual: pd.DataFrame) -> dict: """兩條 OLS 迴歸：負載上升、容量率下降。""" X = sm.add_constant(annual["AD_year"].values) m_load = sm.OLS(annual["peak_load_MW"].values, X).fit() m_margin = sm.OLS(annual["reserve_margin_pct"].values, X).fit() print(f"\n[長期趨勢]") print(f" 尖峰負載 每十年 +{m_load.params[1]*10:.0f} MW " f"(R² = {m_load.rsquared:.2f})") print(f" 備用容量率 每十年 {m_margin.params[1]*10:+.2f} 百分點 " f"(R² = {m_margin.rsquared:.2f})") return {"load_model": m_load, "margin_model": m_margin} def m5_summer_alert_days(daily: pd.DataFrame) -> pd.DataFrame: """夏季（6–9 月）台電燈號分布統計。""" 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[夏季警戒日分布]") print(cross.to_string()) return cross def m5_plot_long_term_scissors(annual: pd.DataFrame, models: dict, savepath: str = "fig5a_taiwan_scissors.png"): """圖 5a：台灣電力在地剪刀效應。""" 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 ax1.plot(x, models["load_model"].predict(sm.add_constant(x)) / 1000, 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="備用容量率") ax2.plot(x, models["margin_model"].predict(sm.add_constant(x)), 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") ax2.axvline(2022, color="blue", lw=1.0, ls=":", alpha=0.5) b_load = models["load_model"].params[1] * 10 b_margin = models["margin_model"].params[1] * 10 ax1.set_title( f"圖 5a｜台灣電力剪刀效應（{annual['AD_year'].min()}–" f"{annual['AD_year'].max()}）\n" f"尖峰負載每十年 +{b_load:.0f} MW｜備用容量率每十年 {b_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) def m5_plot_summer_heatmap(daily: pd.DataFrame, savepath: str = "fig5b_summer_heatmap.png"): """圖 5b：夏季每日備轉容量率熱地圖（以台電燈號離散色階）。""" 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) 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) def run_module_5() -> None: """模組五主流程。""" _banner("模組五｜台灣夏季備轉容量率惡化趨勢") annual = m5_load_annual_reserve() daily = m5_load_daily_operating_reserve() models = m5_long_term_scissors_regression(annual) m5_summer_alert_days(daily) m5_plot_long_term_scissors(annual, models) m5_plot_summer_heatmap(daily) # ============================================================================ # PART 5｜主流程 # ============================================================================ MODULES = { "1": ("模組一｜京都櫻花剪刀效應", run_module_1), "2": ("模組二｜勃艮第葡萄歷史斷層", run_module_2), "4": ("模組四｜北極海冰死亡螺旋", run_module_4), "5": ("模組五｜台灣電力在地剪刀", run_module_5), } def main() -> None: """命令列進入點；支援執行單一模組或全部模組。""" args = sys.argv[1:] if not args: selected = list(MODULES.keys()) else: selected = [a for a in args if a in MODULES] invalid = [a for a in args if a not in MODULES] if invalid: print(f"[警告] 無效模組代號：{invalid}；可用：{list(MODULES.keys())}") if not selected: print("無可執行之模組，結束。") return print("=" * 72) print(f"執行模組：{', '.join(selected)}") print("=" * 72) failures = [] for code in selected: label, fn = MODULES[code] try: fn() except Exception as e: failures.append((code, label, str(e))) print(f"\n[❌ 模組 {code} 失敗] {type(e).__name__}: {e}") continue # 收尾摘要 print("\n" + "=" * 72) if not failures: print(f"✅ 全部 {len(selected)} 個模組執行完成") else: ok = len(selected) - len(failures) print(f"⚠ 完成 {ok}/{len(selected)}；失敗 {len(failures)} 個：") for code, label, msg in failures: print(f" - 模組 {code}（{label}）：{msg}") print("=" * 72) if __name__ == "__main__": main()