核心算法代码.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"/"02_chirps-drought-hotspots",CASE_ROOT/"02_chirps-drought-hotspots"]
for _parent in [CASE_ROOT,*CASE_ROOT.parents]:
    _candidates.extend([_parent/"cases"/"02_chirps-drought-hotspots",_parent/"02_chirps-drought-hotspots"])
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 spi3,drought_frequency,moran_i,lisa_quadrants,lisa_statistics
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(102); months=np.arange(120); ny,nx=8,10
season=65+25*np.sin(2*np.pi*months[:,None]/12); spatial=np.linspace(-12,12,ny*nx)[None,:]
precip=np.maximum(1,season+spatial+rng.normal(0,16,(120,ny*nx)))
Z=spi3(precip,period=3,min_samples=5); freq=drought_frequency(Z); W=np.zeros((ny*nx,ny*nx))
for r in range(ny):
 for c in range(nx):
  i=r*nx+c
  for dr,dc in ((-1,0),(1,0),(0,-1),(0,1)):
   rr,cc=r+dr,c+dc
   if 0<=rr<ny and 0<=cc<nx: W[i,rr*nx+cc]=1
W=W/np.maximum(W.sum(1,keepdims=True),1); annual=np.nanmean(Z.reshape(10,12,-1),axis=1)[-1]

print("数据状态：固定种子合成月降水（不是 CHIRPS 实测栅格）")
print(f"月份数：{months.size}；降水矩阵形状：{precip.shape}；网格：{ny}×{nx}；缺测数：{int(np.isnan(precip).sum())}")

fig,ax=plt.subplots(figsize=(9,4.5)); im=ax.imshow(Z[-1].reshape(ny,nx),cmap="BrBG",vmin=-2.2,vmax=2.2); fig.colorbar(im,ax=ax,label="SPI-3（标准化指数）"); ax.set_title("合成网格：末月 SPI-3 干旱异常"); style_synthetic_grid_axes(ax,nx,ny); fig.tight_layout(); fig.savefig(out_dir/"01_spatiotemporal_overview.png",dpi=160); plt.show()

I=moran_i(freq,W); fig,ax=plt.subplots(figsize=(7,4.5)); im=ax.imshow(freq.reshape(ny,nx),cmap="YlOrRd",vmin=0,vmax=1); fig.colorbar(im,ax=ax,label="干旱发生频率（比例）"); ax.set_title(f"合成网格：全期干旱频率（Global Moran I={I:.3f}；未检验显著性）"); style_synthetic_grid_axes(ax,nx,ny); fig.tight_layout(); fig.savefig(out_dir/"02_trend_significance.png",dpi=160); plt.show()

# Local Moran's I permutation test + Benjamini–Hochberg FDR correction.
lisa=lisa_statistics(annual,W,permutations=4999,seed=102,fdr_alpha=.05); quad=lisa_quadrants(annual,W)
lisa_table=pd.DataFrame({"cell_id":np.arange(ny*nx),"analysis_field":"last_year_12month_mean_SPI3","local_I":lisa["local_i"],"p_value":lisa["p_value"],"q_value_BH":lisa["q_value"],"quadrant":quad,"significant_FDR_0_05":lisa["significant"]}); lisa_table.to_csv(CASE_ROOT/"outputs"/"student_run"/"lisa_cell_diagnostics.csv",index=False)
colors={"HH":"#d73027","LL":"#4575b4","HL":"#fdae61","LH":"#74add1"}; quadrant_names={"HH":"高-高 HH","LL":"低-低 LL","HL":"高-低 HL","LH":"低-高 LH"}; fig,ax=plt.subplots(figsize=(9,4.5)); ax.imshow(np.zeros((ny,nx)),cmap="Greys",alpha=.1)
for q,col in colors.items():
 idx=np.where((quad==q)&lisa["significant"])[0]; ax.scatter(idx%nx,idx//nx,c=col,label=f"{quadrant_names[q]}（FDR 校正 q≤0.05）",s=45)
idx=np.where(~lisa["significant"])[0]; ax.scatter(idx%nx,idx//nx,c="#bdbdbd",label="未达显著性",s=24)
if not lisa["significant"].any(): fig.text(.5,.015,"FDR 校正 q≤0.05：本组合成样例没有显著聚集像元",ha="center",va="bottom",fontsize=9)
ax.legend(title="LISA 聚类象限",fontsize=8,loc="center left",bbox_to_anchor=(1.02,.5)); ax.set_title("合成网格：末年12个月平均 SPI-3 的局部 Moran’s I（置换检验 + BH-FDR）"); style_synthetic_grid_axes(ax,nx,ny); fig.tight_layout(rect=(0,.06,.78,1)); fig.savefig(out_dir/"03_spatial_diagnostic.png",dpi=160); plt.show()

fig,ax=plt.subplots(figsize=(7,4.5)); th=[-.5,-1,-1.5,-2]; vals=[drought_frequency(Z,t).mean() for t in th]; ax.plot(th,vals,"o-"); ax.invert_xaxis(); ax.set(title="干旱阈值敏感性",xlabel="SPI-3 阈值",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([{"global_moran_I_drought_frequency":I,"mean_drought_frequency":float(np.mean(freq)),"max_drought_frequency":float(np.max(freq)),"spi_threshold":-1.0,"lisa_fdr_significant_share":float(np.mean(lisa["significant"])),"spi_method":"calendar-month Gamma + zero-probability mixture; short samples use documented z-score fallback","data_status":"synthetic_demo","lisa_permutations":4999,"lisa_fdr_alpha":0.05,"spi_calibration_start_index":0}])
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)))