核心算法代码.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,CASE_ROOT/"cases"/"04_modis-vegetation-phenology",CASE_ROOT/"04_modis-vegetation-phenology"]
for _parent in [CASE_ROOT,*CASE_ROOT.parents]:
    _candidates.extend([_parent/"cases"/"04_modis-vegetation-phenology",_parent/"04_modis-vegetation-phenology"])
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; run from the case root or notebook folder.")
sys.path.insert(0,str(CASE_ROOT/"src"))
from core_model import smooth_ndvi,extract_phenology,extract_phenology_by_year
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)
plt.rcParams["font.sans-serif"]=[_font_name,"DejaVu Sans"]
plt.rcParams["axes.unicode_minus"]=False
rng=np.random.default_rng(104); n_years,n_composites,ny,nx=3,23,7,9; t=np.arange(n_composites)
phase=np.linspace(-.5,.5,ny*nx); ndvi=.48+.26*np.sin(2*np.pi*(t[:,None]/23+phase[None,:]))+rng.normal(0,.035,(n_composites,ny*nx)); ndvi=np.clip(ndvi,0.05,.95)
ndvi_years=np.stack([np.clip(ndvi+rng.normal(0,.012,ndvi.shape),0.05,.95) for _ in range(n_years)])
sm=np.stack([smooth_ndvi(ndvi_years[y],9) for y in range(n_years)]); ph=extract_phenology_by_year(sm,.2)
years=np.arange(2020,2020+n_years); pixel=20

print("数据状态：固定种子合成 NDVI（不是 MODIS 实测产品）")
print(f"年份：{years.tolist()}；NDVI 数组形状（年、合成期、像元）：{ndvi_years.shape}；缺测数：{int(np.isnan(ndvi_years).sum())}")

fig,ax=plt.subplots(figsize=(9,4.5));
for yi,year in enumerate(years): ax.plot(t,ndvi_years[yi,:,pixel],".",alpha=.35); ax.plot(t,sm[yi,:,pixel],lw=2,label=f"{year} 年平滑曲线")
ax.axhline(ph["threshold"][:,pixel].mean(),ls="--",color="k",label="振幅的20%阈值"); ax.axvline(ph["sos"][0,pixel],color="green",ls=":",label="SOS（首年）"); ax.axvline(ph["eos"][0,pixel],color="brown",ls=":",label="EOS（首年）"); ax.legend(ncol=2,fontsize=8); ax.set(title="合成 NDVI 年度曲线与逐年物候提取",xlabel="年内16日合成序号",ylabel="NDVI"); 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(np.nanmean(ph["sos"],axis=0).reshape(ny,nx),cmap="YlGn"); fig.colorbar(im,ax=ax,label="平均 SOS 合成期索引（0–22）"); ax.set_title("合成网格：2020–2022 年平均 SOS 空间分布"); 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));
for yi,year in enumerate(years): ax.scatter(ph["sos"][yi],ph["eos"][yi],c=ph["amplitude"][yi],cmap="viridis",s=38,label=str(year),alpha=.7)
ax.set(title="逐年 SOS–EOS 关系（颜色为振幅）",xlabel="SOS 合成期索引（0–22）",ylabel="EOS 合成期索引（0–22）"); ax.legend(title="年份"); fig.tight_layout(); fig.savefig(out_dir/"03_spatial_diagnostic.png",dpi=160); plt.show()

fig,ax=plt.subplots(figsize=(7,4.5)); fr=[.1,.2,.3,.4]; sos=[extract_phenology_by_year(sm,f)["yearly_mean_sos"].mean() for f in fr]; ax.plot(fr,sos,"o-"); ax.set(title="逐年物候提取的阈值敏感性",xlabel="阈值比例",ylabel="平均 SOS 合成期索引（0–22）"); fig.tight_layout(); fig.savefig(out_dir/"04_sensitivity.png",dpi=160); plt.show()

from pathlib import Path
import shutil
summary_df=pd.DataFrame({"year":years,"mean_SOS_composite":ph["yearly_mean_sos"],"mean_EOS_composite":ph["yearly_mean_eos"],"mean_peak_NDVI":np.nanmean(ph["peak"],axis=1),"threshold_fraction":.2})
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)))