核心算法代码.ipynb

# 图表字体：优先使用案例库随包提供的开源简体中文字体，保证不同电脑导出中文不变方框。
from pathlib import Path
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import font_manager

_font_relative_path = Path("assets") / "fonts" / "NotoSansCJKsc-Regular.otf"
_font_search_roots = [Path.cwd(), *Path.cwd().parents]
_font_path = next((base / _font_relative_path for base in _font_search_roots
                   if (base / _font_relative_path).is_file()), None)
if _font_path is not None:
    font_manager.fontManager.addfont(str(_font_path))
    _font_name = font_manager.FontProperties(fname=str(_font_path)).get_name()
    _font_source = "案例库随包字体"
else:
    # 单独拷贝 Notebook 时仍可使用已安装的中文字体；找不到时明确报错，不生成方框图。
    _system_font_candidates = ["Microsoft YaHei", "Noto Sans SC", "Noto Sans CJK SC",
                               "Source Han Sans SC", "SimHei", "DengXian", "SimSun",
                               "WenQuanYi Zen Hei", "Arial Unicode MS"]
    _available_fonts = {item.name for item in font_manager.fontManager.ttflist}
    _font_name = next((name for name in _system_font_candidates if name in _available_fonts), None)
    if _font_name is None:
        raise FileNotFoundError(
            "未找到中文字体。请从案例库根目录运行 Notebook，并保留 assets/fonts/NotoSansCJKsc-Regular.otf。"
        )
    _font_source = "系统已安装字体"

matplotlib.rcParams["font.family"] = "sans-serif"
matplotlib.rcParams["font.sans-serif"] = [_font_name, "DejaVu Sans"]
matplotlib.rcParams["axes.unicode_minus"] = False
print(f"中文字体：{_font_name}（{_font_source}）")

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from pathlib import Path
import sys
CASE_ROOT=Path.cwd()
_candidates=[CASE_ROOT,*CASE_ROOT.parents]
for _parent in [CASE_ROOT,*CASE_ROOT.parents]:
    _candidates.extend([_parent/"cases"/'03_era5-heat-extremes',_parent/'03_era5-heat-extremes'])
for _candidate in _candidates:
    if (_candidate/"src"/"core_model.py").is_file(): CASE_ROOT=_candidate; break
else: raise FileNotFoundError("Cannot locate case root or src/core_model.py")
sys.path.insert(0,str(CASE_ROOT/"src"))
from core_model import percentile_threshold, heatwave_mask, longest_run
out_dir=CASE_ROOT/"notebook"/"figures"; package_figures=Path('figures'); out_dir.mkdir(parents=True,exist_ok=True)
def style_synthetic_grid_axes(ax,nx,ny):
    xt=np.unique(np.linspace(0,nx-1,min(nx,5),dtype=int)); yt=np.unique(np.linspace(0,ny-1,min(ny,5),dtype=int))
    ax.set_xticks(xt,[str(v+1) for v in xt]); ax.set_yticks(yt,[str(v+1) for v in yt])
    ax.set_xlim(-.5,nx-.5); ax.set_ylim(ny-.5,-.5); ax.set_aspect("equal",adjustable="box")
    ax.set_xlabel("网格列（像元序号）"); ax.set_ylabel("网格行（自上而下）"); ax.tick_params(labelsize=8)
rng=np.random.default_rng(103); days=np.arange(2*365); ny,nx=8,10
season=27+8*np.sin(2*np.pi*(days%365)/365)[:,None]; spatial=np.linspace(-3,3,ny*nx)[None,:]
tmax=season+spatial+rng.normal(0,2.2,(len(days),ny*nx)); baseline=tmax[:365]; threshold=percentile_threshold(baseline,90,axis=0); mask=heatwave_mask(tmax,threshold,3); eval_tmax=tmax[365:]; eval_mask=mask[365:]; eval_runs=longest_run(eval_mask)

print("数据状态：固定种子合成日最高温（不是 ERA5-Land 实测数据）")
print(f"日数：{days.size}；温度矩阵形状：{tmax.shape}；网格：{ny}×{nx}；单位：°C")

fig,ax=plt.subplots(figsize=(9,4.5)); ax.plot(days[:730],tmax[:730,30],lw=.8,label="第31网格最高气温"); ax.axhline(threshold[30],color="crimson",ls="--",label="第90百分位阈值"); ax.fill_between(days[:730],threshold[30],tmax[:730,30],where=tmax[:730,30]>threshold[30],color="crimson",alpha=.25); ax.legend(); ax.set(title="合成教学数据：极端高温阈值",xlabel="日序号（合成序列）",ylabel="最高气温（℃）"); fig.tight_layout(); fig.savefig(out_dir/"01_spatiotemporal_overview.png",dpi=160); plt.show()

fig,ax=plt.subplots(figsize=(7,4.5)); im=ax.imshow(eval_mask.sum(0).reshape(ny,nx),cmap="OrRd"); fig.colorbar(im,ax=ax,label="热浪日数（天）"); ax.set_title("合成末年（365日）热浪日数（连续≥3日）"); style_synthetic_grid_axes(ax,nx,ny); fig.tight_layout(); fig.savefig(out_dir/"02_trend_significance.png",dpi=160); plt.show()

fig,ax=plt.subplots(figsize=(7,4.5)); ax.hist(eval_runs,bins=np.arange(eval_runs.max()+2)-.5,color="#e67e22"); ax.set(title="合成网格：各像元最长连续热浪日数分布",xlabel="持续天数",ylabel="网格像元数"); fig.tight_layout(); fig.savefig(out_dir/"03_spatial_diagnostic.png",dpi=160); plt.show()

fig,ax=plt.subplots(figsize=(7,4.5)); qs=[85,90,95]; vals=[]
for q in qs: vals.append(heatwave_mask(eval_tmax,percentile_threshold(baseline,q),3).mean())
ax.plot(qs,vals,"o-"); ax.set(title="百分位阈值敏感性",xlabel="百分位数",ylabel="被判定为热浪的日数比例"); fig.tight_layout(); fig.savefig(out_dir/"04_sensitivity.png",dpi=160); plt.show()

from pathlib import Path
import shutil
summary_df=pd.DataFrame([{"mean_threshold_C":float(np.nanmean(threshold)),"heatwave_fraction":float(eval_mask.mean()),"max_run_days":int(eval_runs.max()),"minimum_duration_days":3,"data_status":"synthetic_demo"}])
out=CASE_ROOT/"outputs"/"student_run"; fig_out=out/"figures"; fig_out.mkdir(parents=True,exist_ok=True)
for p in sorted(out_dir.glob("*.png")): shutil.copy2(p,fig_out/p.name)
summary_df.to_csv(out/"summary.csv",index=False); display(summary_df.round(4))

# 交付图同步：从图件追溯表读取本次真实运行图，更新 Notebook 预览和 PPT 使用的四张标准图。
from pathlib import Path
import csv
import shutil
from IPython.display import Image, display

_figure_dir = CASE_ROOT / "notebook" / "figures"
_trace_path = CASE_ROOT / "notebook" / "FIGURE_TRACEABILITY.csv"
with _trace_path.open("r", encoding="utf-8-sig", newline="") as _stream:
    _trace_rows = list(csv.DictReader(_stream))
if len(_trace_rows) != 4:
    raise ValueError(f"图件追溯表应有4条记录，实际 {len(_trace_rows)} 条")

for _row in _trace_rows:
    _source_rel = Path(_row["student_output"])
    if _source_rel.is_absolute() or ".." in _source_rel.parts:
        raise ValueError(f"图件来源路径必须位于案例目录内：{_source_rel}")
    _source = CASE_ROOT / _source_rel
    if not _source.is_file():
        raise FileNotFoundError(f"本次运行没有生成追溯表登记的图：{_source}")
    _target_name = (_row.get("delivery_output") or "").strip()
    if not _target_name or Path(_target_name).name != _target_name:
        raise ValueError(f"交付图名无效：{_target_name!r}")
    _target = _figure_dir / _target_name
    if _source.resolve() != _target.resolve():
        shutil.copy2(_source, _target)
    print(f"交付图已更新：{_target_name} ← {_source_rel.as_posix()}")
    display(Image(filename=str(_target)))