初始化协作仓库：赛题材料、初步试解与协作说明
- 大文件（附件原始包、data/）由 files.ddwsgood.top 承载，不入库 - 新增 .gitignore/.gitattributes/COLLAB.md
This commit is contained in:
mathorcup committed
1 parent
d87c3b54eb
commit
ac3e481e11
75 files changed
+2389
No files matched your search
@@ -0,0 +1,267 @@
|
||||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
全球冰川灾害编目（TPDC 3ecd6ef8，DOI 10.11888/Cryos.tpdc.272812）审计与试算 —— 修正版 v2。
|
||||
AI 声明：DeepSeek-V4.1-Flash（DeepSeek），2026-10-23。需人工复核。
|
||||
|
||||
v2 修正（对应主线程审查意见）：
|
||||
1 去重键含完整规范化日期/区间；仅到年的记录只标“疑似重复”，不删除；统一称“候选去重”。
|
||||
2 时间解析顺序：区间 → 完整日期 → 年月 → 年份；非完整日期不记 day；保留原文。
|
||||
3 有效事件行 = 有 ID 或 Site 的真实数据行；缺日期不丢，单独标 missing；合并单元格只在真实合并范围内还原，不做无条件 ffill。
|
||||
4 audit 分 raw / candidate_unique 两套口径输出。
|
||||
5 同时保留 Trigger 与 Triggers 两个字段。
|
||||
|
||||
输出（覆盖前已在 backups/YYYYMMDD 备份）：
|
||||
结果/catalog_slope_failures_clean.csv
|
||||
结果/catalog_glof_clean.csv
|
||||
结果/catalog_audit.json
|
||||
结果/runout_angle_stats.json
|
||||
"""
|
||||
import os, io, re, json, zipfile, datetime as dt
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import openpyxl
|
||||
|
||||
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
|
||||
ZIP = os.path.abspath(os.path.join(ROOT, "..", "TPDC_补充数据", "下载",
|
||||
"3ecd6ef8-fb24-44f0-b9f1-c8839fe83760.zip"))
|
||||
RES = os.path.join(ROOT, "结果")
|
||||
os.makedirs(RES, exist_ok=True)
|
||||
RNG = np.random.default_rng(20261023) # 固定留出，保证可复现
|
||||
|
||||
z = zipfile.ZipFile(ZIP)
|
||||
names = [n for n in z.namelist() if n.endswith(".xlsx")]
|
||||
F_SF = [n for n in names if "slope failures" in n][0]
|
||||
F_GL = [n for n in names if "Moraine" in n][0]
|
||||
|
||||
def load_sheet(path, sheet):
|
||||
"""用 openpyxl 读；仅对真实 merged range 用左上角值回填，其余不 ffill。"""
|
||||
wb = openpyxl.load_workbook(io.BytesIO(z.read(path)), data_only=True)
|
||||
ws = wb[sheet]
|
||||
merged = {}
|
||||
for rng in ws.merged_cells.ranges:
|
||||
v = ws.cell(rng.min_row, rng.min_col).value
|
||||
for r in range(rng.min_row, rng.max_row + 1):
|
||||
for c in range(rng.min_col, rng.max_col + 1):
|
||||
merged[(r, c)] = v
|
||||
rows = list(ws.iter_rows(values_only=False))
|
||||
header = [c.value for c in rows[0]]
|
||||
recs, n_merged_filled = [], 0
|
||||
for r_i, row in enumerate(rows[1:], start=2):
|
||||
vals = []
|
||||
for c_i, cell in enumerate(row, start=1):
|
||||
v = cell.value
|
||||
if (v is None or v == "") and (r_i, c_i) in merged:
|
||||
v = merged[(r_i, c_i)]; n_merged_filled += 1
|
||||
vals.append(v)
|
||||
recs.append(vals)
|
||||
df = pd.DataFrame(recs, columns=[str(h).strip() if h is not None else f"col{i}" for i, h in enumerate(header)])
|
||||
df.attrs["merged_filled_cells"] = n_merged_filled
|
||||
df.attrs["merged_ranges"] = len(ws.merged_cells.ranges)
|
||||
return df
|
||||
|
||||
DATE_RE = re.compile(r"^\d{4}[-/]\d{1,2}[-/]\d{1,2}([ T]\d{1,2}:\d{2}(:\d{2})?)?$")
|
||||
YM_RE = re.compile(r"^\d{4}[-/]\d{1,2}$")
|
||||
Y_RE = re.compile(r"^\d{4}$")
|
||||
|
||||
def parse_time(v):
|
||||
"""返回 (原文, 精度, 起始年, 结束年, 规范化日期键)。
|
||||
精度：hour / day / month / year / range / other / missing。"""
|
||||
if v is None or (isinstance(v, float) and np.isnan(v)):
|
||||
return None, "missing", None, None, None
|
||||
if isinstance(v, (dt.datetime, pd.Timestamp)):
|
||||
s = v.strftime("%Y-%m-%d %H:%M:%S")
|
||||
prec = "hour" if (v.hour or v.minute or v.second) else "day"
|
||||
return s, prec, v.year, v.year, s[:10]
|
||||
s = str(v).strip()
|
||||
if s in ("", "None", "nan", "NaN"):
|
||||
return None, "missing", None, None, None
|
||||
# 1) 区间：出现两个 4 位年份（2014/2-2015/4、2000-2001、2014/2—2015/4）
|
||||
yrs = [int(t) for t in re.findall(r"(1[5-9]\d{2}|20\d{2})", s)]
|
||||
if len(yrs) >= 2:
|
||||
return s, "range", min(yrs), max(yrs), None
|
||||
# 2) 完整日期（可带时间）
|
||||
if DATE_RE.match(s):
|
||||
d = s.replace("/", "-")
|
||||
prec = "hour" if (" " in d or "T" in d) else "day"
|
||||
try:
|
||||
y = int(d[:4])
|
||||
except Exception:
|
||||
return s, "other", None, None, None
|
||||
return s, prec, y, y, d[:10]
|
||||
# 3) 年月
|
||||
if YM_RE.match(s):
|
||||
y = int(s[:4])
|
||||
return s, "month", y, y, None
|
||||
# 4) 年份
|
||||
if Y_RE.match(s):
|
||||
return s, "year", int(s), int(s), None
|
||||
# 5) 其他（含中文、英文月份、约等）
|
||||
return s, "other", (yrs[0] if yrs else None), (yrs[0] if yrs else None), None
|
||||
|
||||
SHEETS = ["Glacier detachments", "Ice avalanches", "Ice-rock avalanches",
|
||||
"Rock avalanches", "GLOF triggered by avalanches", "GLOF triggered by calving"]
|
||||
|
||||
rows, sheet_diag = [], {}
|
||||
for sh in SHEETS:
|
||||
df = load_sheet(F_SF, sh)
|
||||
n_raw = len(df)
|
||||
n_evt = n_time_missing = 0
|
||||
for _, r in df.iterrows():
|
||||
eid = r.get("ID"); site = r.get("Site")
|
||||
# 真实事件行：ID 或 Site 至少一个非空（整行空行不计）
|
||||
if (pd.isna(eid) and pd.isna(site)):
|
||||
continue
|
||||
n_evt += 1
|
||||
t_raw = r.get("Time")
|
||||
ts, prec, y0, y1, dkey = parse_time(t_raw)
|
||||
if prec == "missing":
|
||||
n_time_missing += 1
|
||||
try:
|
||||
lon = float(r.get("Longitude")); lat = float(r.get("Latitude"))
|
||||
except Exception:
|
||||
lon = lat = np.nan
|
||||
rows.append(dict(
|
||||
sheet=sh, src_ID=eid, Site=str(site).strip() if not pd.isna(site) else None,
|
||||
time_raw=ts, time_prec=prec, year_start=y0, year_end=y1, date_key=dkey,
|
||||
mechanism=r.get("Type"), lon=lon, lat=lat,
|
||||
angle_reach_deg=r.get("Angle of reach °"),
|
||||
source_area_km2=r.get("Source area (km2)"),
|
||||
volume_Mm3=r.get("Volume (106m3)"),
|
||||
source_slope_deg=r.get("Source slope (°)") if "Source slope (°)" in df.columns else r.get("Source slope °"),
|
||||
horizontal_km=r.get("Horizontal distance (km)"),
|
||||
affected_km2=r.get("Affected area (km2)"),
|
||||
cascade=r.get("Cascade hazards"),
|
||||
trigger_plural=r.get("Triggers"), # 多数表用 Triggers
|
||||
trigger_singular=r.get("Trigger"), # Ice-rock 表用 Trigger
|
||||
refs=r.get("References")))
|
||||
sheet_diag[sh] = dict(raw_rows=n_raw, event_rows=n_evt, time_missing=n_time_missing,
|
||||
merged_ranges=df.attrs["merged_ranges"],
|
||||
merged_cells_filled=df.attrs["merged_filled_cells"])
|
||||
|
||||
sf = pd.DataFrame(rows)
|
||||
for c in ["angle_reach_deg", "source_area_km2", "volume_Mm3", "source_slope_deg",
|
||||
"horizontal_km", "affected_km2"]:
|
||||
sf[c] = pd.to_numeric(sf[c], errors="coerce")
|
||||
sf["trigger"] = sf.trigger_plural.where(sf.trigger_plural.notna(), sf.trigger_singular)
|
||||
|
||||
def norm(s):
|
||||
return re.sub(r"[\s\-_]+", "", str(s).lower())
|
||||
|
||||
# ---- 去重：日期精确的才作为候选重复；年份精度的只标疑似
|
||||
sf["dedupe_key_exact"] = np.where(
|
||||
sf.date_key.notna(),
|
||||
sf.Site.map(norm) + "|" + sf.lat.round(3).astype(str) + "|" + sf.lon.round(3).astype(str) + "|" + sf.date_key.fillna(""),
|
||||
None)
|
||||
sf["dedupe_key_year"] = np.where(
|
||||
sf.date_key.isna() & sf.year_start.notna(),
|
||||
sf.Site.map(norm) + "|" + sf.lat.round(3).astype(str) + "|" + sf.lon.round(3).astype(str) + "|Y" + sf.year_start.fillna(0).astype(int).astype(str),
|
||||
None)
|
||||
sf["is_dupe_exact"] = sf.dedupe_key_exact.notna() & sf.duplicated("dedupe_key_exact", keep="first")
|
||||
sf["is_dupe_suspect"] = (~sf.is_dupe_exact) & sf.dedupe_key_year.notna() & sf.duplicated("dedupe_key_year", keep="first")
|
||||
sf["candidate_unique"] = ~(sf.is_dupe_exact | sf.is_dupe_suspect)
|
||||
|
||||
sf["is_hma"] = sf.lat.between(25, 45) & sf.lon.between(60, 105)
|
||||
sf["is_himalaya"] = sf.lat.between(26, 31) & sf.lon.between(75, 95)
|
||||
|
||||
# ---- GLOF 表
|
||||
gl_df = load_sheet(F_GL, "Sheet1")
|
||||
gl_rows = []
|
||||
for _, r in gl_df.iterrows():
|
||||
lake, date = r.get("Glacial lake"), r.get("Date")
|
||||
if pd.isna(lake) and pd.isna(date):
|
||||
continue
|
||||
ts, prec, y0, y1, dkey = parse_time(date)
|
||||
try:
|
||||
lon = float(r.get("Longitude")); lat = float(r.get("Latitude"))
|
||||
except Exception:
|
||||
lon = lat = np.nan
|
||||
gl_rows.append(dict(src_ID=r.get("ID"), country=r.get("Country"),
|
||||
lake=str(lake).strip() if not pd.isna(lake) else None,
|
||||
time_raw=ts, time_prec=prec, year_start=y0, year_end=y1, date_key=dkey,
|
||||
lake_type=r.get("Lake type"), lon=lon, lat=lat, triggers=r.get("Triggers"),
|
||||
volume_Mm3=pd.to_numeric(r.get("Release volume/106m3"), errors="coerce"),
|
||||
peak_discharge=r.get("Peak discharge/m3·s-1"), deaths=r.get("Deaths"),
|
||||
refs=r.get("References")))
|
||||
gl = pd.DataFrame(gl_rows)
|
||||
gl["is_hma"] = gl.lat.between(25, 45) & gl.lon.between(60, 105)
|
||||
gl["is_himalaya"] = gl.lat.between(26, 31) & gl.lon.between(75, 95)
|
||||
|
||||
def summarize(df, tag):
|
||||
raw = df
|
||||
uniq = df[df.candidate_unique] if "candidate_unique" in df else df
|
||||
return dict(
|
||||
raw=dict(n=int(len(raw)), by_sheet=raw.sheet.value_counts().to_dict() if "sheet" in raw else {},
|
||||
by_mechanism=raw.mechanism.astype(str).value_counts().to_dict() if "mechanism" in raw else {},
|
||||
by_precision=raw.time_prec.value_counts().to_dict(),
|
||||
hma=int(raw.is_hma.sum()), himalaya=int(raw.is_himalaya.sum())),
|
||||
candidate_unique=dict(n=int(len(uniq)),
|
||||
by_sheet=uniq.sheet.value_counts().to_dict() if "sheet" in uniq else {},
|
||||
by_mechanism=uniq.mechanism.astype(str).value_counts().to_dict() if "mechanism" in uniq else {},
|
||||
by_precision=uniq.time_prec.value_counts().to_dict(),
|
||||
hma=int(uniq.is_hma.sum()), himalaya=int(uniq.is_himalaya.sum())),
|
||||
duplicates=dict(exact_removed=int(df.is_dupe_exact.sum()) if "is_dupe_exact" in df else 0,
|
||||
suspect_year_only=int(df.is_dupe_suspect.sum()) if "is_dupe_suspect" in df else 0))
|
||||
|
||||
audit = dict(
|
||||
script_version="v2 (2026-10-23 correction)",
|
||||
source_zip=os.path.basename(ZIP), doi="10.11888/Cryos.tpdc.272812",
|
||||
metadata_claim=dict(slope_failures=737, glof=380),
|
||||
claim_note="737/380 是元数据描述；本表只报可复现的实际计数，不为对齐宣传数修改数据。",
|
||||
sheet_diagnostics=sheet_diag,
|
||||
slope_failures=summarize(sf, "slope"),
|
||||
slope_fields=dict(volume=int(sf.volume_Mm3.notna().sum()),
|
||||
volume_and_H=int((sf.volume_Mm3.notna() & sf.horizontal_km.notna()).sum()),
|
||||
angle=int(sf.angle_reach_deg.notna().sum())),
|
||||
glof=summarize(gl, "glof"),
|
||||
glof_fields=dict(volume=int(gl.volume_Mm3.notna().sum())),
|
||||
)
|
||||
|
||||
sf.to_csv(os.path.join(RES, "catalog_slope_failures_clean.csv"), index=False, encoding="utf-8-sig")
|
||||
gl.to_csv(os.path.join(RES, "catalog_glof_clean.csv"), index=False, encoding="utf-8-sig")
|
||||
with open(os.path.join(RES, "catalog_audit.json"), "w", encoding="utf-8") as f:
|
||||
json.dump(audit, f, ensure_ascii=False, indent=2)
|
||||
|
||||
print("== sheet diagnostics")
|
||||
for k, v in sheet_diag.items():
|
||||
print(f" {k:<34} raw={v['raw_rows']:>4} events={v['event_rows']:>4} no_time={v['time_missing']:>4} "
|
||||
f"merged_ranges={v['merged_ranges']} merged_cells={v['merged_cells_filled']}")
|
||||
print("\n== slope failures")
|
||||
print(" raw:", audit["slope_failures"]["raw"])
|
||||
print(" cand-unique:", audit["slope_failures"]["candidate_unique"])
|
||||
print(" dup:", audit["slope_failures"]["duplicates"])
|
||||
print("\n== GLOF")
|
||||
print(" raw:", audit["glof"]["raw"])
|
||||
print(" cand-unique:", audit["glof"]["candidate_unique"])
|
||||
|
||||
# ---------------- 传播/摩擦（样本因去重口径变化，重算一次）
|
||||
print("\n" + "=" * 60)
|
||||
sub = sf[sf.candidate_unique & sf.volume_Mm3.gt(0) & sf.horizontal_km.gt(0)]
|
||||
print("体积+距离齐全（候选唯一）：")
|
||||
print(sub.groupby("sheet").size().to_string())
|
||||
ice = sub[sub.sheet.isin(["Ice-rock avalanches", "Glacier detachments"])].copy()
|
||||
out = dict(script_version="v2", n_ice=len(ice), n_ice_hma=int(ice.is_hma.sum()))
|
||||
if len(ice) >= 12:
|
||||
x = np.log10(ice.volume_Mm3.values); y = np.log10(ice.horizontal_km.values)
|
||||
idx = RNG.permutation(len(ice)); cut = int(0.7*len(ice))
|
||||
tr, te = idx[:cut], idx[cut:]
|
||||
b, a = np.polyfit(x[tr], y[tr], 1)
|
||||
pred = a + b*x[te]
|
||||
rmse = float(np.sqrt(np.mean((pred - y[te])**2)))
|
||||
r2 = float(1 - np.sum((y[te]-pred)**2)/np.sum((y[te]-y[te].mean())**2))
|
||||
out.update(dict(fit=(float(a), float(b)), holdout_rmse_log10=rmse, holdout_r2=r2, n_test=int(len(te))))
|
||||
print(f"留出：log10H={a:.3f}+{b:.3f}log10V，RMSE={rmse:.2f}，R²={r2:.2f}（n_test={len(te)}）")
|
||||
else:
|
||||
out["note"] = "n<12, no fit"
|
||||
|
||||
ang = sf[sf.candidate_unique & sf.angle_reach_deg.between(0.5, 60)]["angle_reach_deg"]
|
||||
mu = np.tan(np.deg2rad(ang))
|
||||
out["angle_reach_deg"] = dict(n=int(len(ang)), p10=float(ang.quantile(.1)),
|
||||
median=float(ang.median()), p90=float(ang.quantile(.9)))
|
||||
out["mu_eq_tan"] = dict(p10=float(mu.quantile(.1)), median=float(mu.median()), p90=float(mu.quantile(.9)))
|
||||
out["note_2026"] = ("2026 事件的 8.85°（=atan(3210/20620)）来自‘模拟/粗网格几何’，"
|
||||
"仅作外部量级对照，不能据此断言真实事件的分位或反推 μ 范围。")
|
||||
print(f"角 of reach：n={len(ang)}，p10={ang.quantile(.1):.1f}°，中位={ang.median():.1f}°，p90={ang.quantile(.9):.1f}°")
|
||||
print("等效 tan：", {k: round(v, 3) for k, v in out["mu_eq_tan"].items()})
|
||||
with open(os.path.join(RES, "runout_angle_stats.json"), "w", encoding="utf-8") as f:
|
||||
json.dump(out, f, ensure_ascii=False, indent=2)
|
||||
@@ -0,0 +1,174 @@
|
||||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
TPDC 增量的最小审计（DeepSeek 子任务 B）。
|
||||
AI 声明：本脚本由 DeepSeek-V4.1-Flash（DeepSeek）于 2026-10-23 生成，需人工复核。
|
||||
|
||||
只做四件事，全部只读 TPDC_补充数据/下载/*.zip（嵌套包解到 /tmp/opencode 自建目录）：
|
||||
1) GFS 灾前预报：源区/口岸逐时次降水率与 6h 累计，核对 lead time、UTC、单位；
|
||||
2) TPCSM：变量、初始化时刻、网格，若可用则取源区降雨；
|
||||
3) 受灾位置点 vs 我们 20.6 km 路径终点；河网 vs 路径；
|
||||
4) 输出 JSON/CSV 到 结果/。
|
||||
不重复下载、不登录、不修改原始 ZIP。
|
||||
"""
|
||||
import os, sys, json, zipfile, shutil, struct
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import netCDF4 as nc
|
||||
|
||||
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
|
||||
TPDC = os.path.abspath(os.path.join(ROOT, "..", "TPDC_补充数据", "下载"))
|
||||
RES = os.path.join(ROOT, "结果")
|
||||
TMP = "/tmp/opencode/tpdc_increment"
|
||||
os.makedirs(RES, exist_ok=True)
|
||||
os.makedirs(TMP, exist_ok=True)
|
||||
|
||||
def unnest(zid, inner_ext):
|
||||
"""解外层 zip，再解出内层 zip 中含指定后缀的目录；返回目录路径。"""
|
||||
dst = os.path.join(TMP, zid)
|
||||
if os.path.isdir(dst) and any(f.endswith(inner_ext) for _, _, fs in os.walk(dst) for f in fs):
|
||||
return dst
|
||||
os.makedirs(dst, exist_ok=True)
|
||||
with zipfile.ZipFile(os.path.join(TPDC, zid + ".zip")) as z:
|
||||
z.extractall(dst)
|
||||
for r, _d, fs in os.walk(dst):
|
||||
for f in fs:
|
||||
if f.endswith(".zip"):
|
||||
with zipfile.ZipFile(os.path.join(r, f)) as z:
|
||||
z.extractall(os.path.join(r, "inner"))
|
||||
return dst
|
||||
|
||||
def find(dst, ext):
|
||||
for r, _d, fs in os.walk(dst):
|
||||
for f in fs:
|
||||
if f.endswith(ext):
|
||||
yield os.path.join(r, f)
|
||||
|
||||
# ---------------------------------------------------------------- 1) GFS
|
||||
print("=" * 70); print("1) GFS 灾前预报（起报 2026-08-25 00Z）")
|
||||
gdir = unnest("147456c9-f740-4259-a2ea-fd7423613510", ".nc")
|
||||
gfile = sorted(find(gdir, ".nc"))[0]
|
||||
ds = nc.Dataset(gfile)
|
||||
lat = np.asarray(ds["latitude"][:]); lon = np.asarray(ds["longitude"][:])
|
||||
fh = np.asarray(ds["forecast_hour"][:])
|
||||
init = ds.getncattr("initialization_time")
|
||||
print(f"nc={os.path.basename(gfile)}\n起报={init} 时次={list(map(int,fh))}\n"
|
||||
f"网格 lat{lat.min():.3f}-{lat.max():.3f} lon{lon.min():.3f}-{lon.max():.3f}")
|
||||
|
||||
def idx(la, lo):
|
||||
return int(np.argmin(np.abs(lat - la))), int(np.argmin(np.abs(lon - lo)))
|
||||
|
||||
SRC = (28.28, 85.52) # 源区（公开 AOI 中心）
|
||||
PORT = (28.2785, 85.3766) # 受灾位置元数据
|
||||
rows = []
|
||||
for tag, (la, lo) in [("source", SRC), ("port", PORT)]:
|
||||
i, j = idx(la, lo)
|
||||
rate = np.asarray(ds["precipitation_rate"][:, i, j], float)
|
||||
acc6 = np.asarray(ds["total_precipitation_6h"][:, i, j], float)
|
||||
acci = np.asarray(ds["total_precipitation_since_initialization"][:, i, j], float)
|
||||
cw = np.asarray(ds["cloud_water_path"][:, i, j], float)
|
||||
pw = np.asarray(ds["precipitable_water"][:, i, j], float)
|
||||
for k in range(len(fh)):
|
||||
rows.append(dict(site=tag, grid_lat=float(lat[i]), grid_lon=float(lon[j]),
|
||||
fh=int(fh[k]),
|
||||
valid_utc=str(np.datetime64("2026-08-25T00:00") + np.timedelta64(int(fh[k]), "h")),
|
||||
precip_rate_mm_h=float(rate[k]), precip_6h_mm=float(acc6[k]),
|
||||
precip_cum_mm=float(acci[k]), cloud_water_path=float(cw[k]),
|
||||
precipitable_water=float(pw[k])))
|
||||
g = pd.DataFrame(rows)
|
||||
g.to_csv(os.path.join(RES, "gfs_precip_source_port.csv"), index=False, encoding="utf-8-sig")
|
||||
print(g.to_string(index=False, float_format=lambda v: f"{v:.2f}"))
|
||||
|
||||
# 事件窗：2026-08-26 10:30 BJT = 02:30 UTC，落在 f026-f032 之间 → 用 f030 最近
|
||||
ev = g[(g.valid_utc >= "2026-08-26T00:00") & (g.valid_utc <= "2026-08-26T12:00")]
|
||||
print("\n事件当日 00-12 UTC 预报（源区/口岸）：")
|
||||
print(ev[["site", "fh", "valid_utc", "precip_rate_mm_h", "precip_6h_mm"]].to_string(index=False))
|
||||
for tag in ["source", "port"]:
|
||||
sub = g[g.site == tag]
|
||||
# 事件日 6h 累计之和（只加 6h 量，不再与 since-init 相加）
|
||||
ev6 = sub[(sub.valid_utc >= "2026-08-26T00:00") & (sub.valid_utc <= "2026-08-26T18:00")]["precip_6h_mm"].sum()
|
||||
print(f"{tag}: 2026-08-26 00-18UTC 的 6h 累计之和 = {ev6:.1f} mm（仅 6h 量求和）")
|
||||
ds.close()
|
||||
|
||||
# ---------------------------------------------------------------- 2) TPCSM
|
||||
print("\n" + "=" * 70); print("2) TPCSM 5 km 模拟")
|
||||
tdir = unnest("be1c4b80-f26e-4f43-a54a-cedbfe06203b", ".nc")
|
||||
tfiles = sorted(find(tdir, ".nc")) + sorted(find(tdir, ".nc4"))
|
||||
print("nc 文件数:", len(tfiles))
|
||||
if tfiles:
|
||||
t = nc.Dataset(tfiles[0])
|
||||
print("file:", os.path.basename(tfiles[0]))
|
||||
print("dims:", {k: len(v) for k, v in t.dimensions.items()})
|
||||
for n, v in t.variables.items():
|
||||
print(f" {n:<22}{v.dimensions} {v.shape} {v.dtype} {getattr(v,'units','')}")
|
||||
print("attrs:", {k: t.getncattr(k) for k in t.ncattrs()})
|
||||
t.close()
|
||||
else:
|
||||
docx = list(find(tdir, ".docx"))
|
||||
print("未找到 nc；docx:", docx)
|
||||
|
||||
# ---------------------------------------------------------------- 3) 点位与河网
|
||||
print("\n" + "=" * 70); print("3) 受灾位置点 / 河网 vs 我们的路径")
|
||||
def shp_points(path):
|
||||
s = open(path, "rb").read()
|
||||
st = struct.unpack("<i", s[32:36])[0]
|
||||
pos, pts = 100, []
|
||||
while pos < len(s) - 8:
|
||||
_num, clen = struct.unpack(">II", s[pos:pos+8]); clen *= 2
|
||||
body = s[pos+8:pos+8+clen]
|
||||
typ = struct.unpack("<i", body[:4])[0]
|
||||
if typ == 1:
|
||||
x, y = struct.unpack("<2d", body[4:20]); pts.append((x, y))
|
||||
pos += 8 + clen
|
||||
return st, pts
|
||||
|
||||
def shp_lines(path):
|
||||
s = open(path, "rb").read()
|
||||
pos, lines = 100, []
|
||||
while pos < len(s) - 8:
|
||||
_n, clen = struct.unpack(">II", s[pos:pos+8]); clen *= 2
|
||||
body = s[pos+8:pos+8+clen]
|
||||
typ = struct.unpack("<i", body[:4])[0]
|
||||
if typ in (3, 5):
|
||||
nparts, npts = struct.unpack("<2i", body[36:44])
|
||||
parts = struct.unpack("<%di" % nparts, body[44:44+4*nparts])
|
||||
p0 = 44 + 4*nparts
|
||||
xy = np.frombuffer(body[p0:p0+16*npts], dtype="<f8").reshape(-1, 2)
|
||||
for k in range(nparts):
|
||||
a = parts[k]; b = parts[k+1] if k+1 < nparts else npts
|
||||
lines.append(xy[a:b])
|
||||
pos += 8 + clen
|
||||
return lines
|
||||
|
||||
pdir = unnest("72aff770-2f75-459d-a5cc-77812f5b32be", ".shp")
|
||||
st, pts = shp_points(sorted(find(pdir, ".shp"))[0])
|
||||
print("受灾点位 shp 类型", st, "点数", len(pts), "坐标", pts[:3])
|
||||
|
||||
rdir = unnest("e7f8d50b-f647-4efe-8d15-a87173c13ab8", ".shp")
|
||||
lines = shp_lines(sorted(find(rdir, ".shp"))[0])
|
||||
allp = np.vstack(lines)
|
||||
print(f"河网线 {len(lines)} 条，节点 {len(allp)}，bbox {allp[:,0].min():.4f}-{allp[:,0].max():.4f}E "
|
||||
f"{allp[:,1].min():.4f}-{allp[:,1].max():.4f}N")
|
||||
|
||||
import math
|
||||
def hav(a, b):
|
||||
R = 6371.0
|
||||
p1, p2 = math.radians(a[1]), math.radians(b[1])
|
||||
dp = p2-p1; dl = math.radians(b[0]-a[0])
|
||||
return 2*R*math.asin(math.sqrt(math.sin(dp/2)**2 + math.cos(p1)*math.cos(p2)*math.sin(dl/2)**2))
|
||||
|
||||
path = pd.read_csv(os.path.join(RES, "q4_path_profile.csv"))
|
||||
end = (float(path.lon.iloc[-1]), float(path.lat.iloc[-1]))
|
||||
dmg = pts[0]
|
||||
print(f"路径终点 ({end[0]:.4f},{end[1]:.4f}) 高程 {path.elev_m.iloc[-1]:.0f} m；"
|
||||
f"受灾点 ({dmg[0]:.4f},{dmg[1]:.4f})；距离 {hav(end, dmg):.2f} km")
|
||||
# 路径每点 1 km 内是否有河网节点（粗判路径是否沿河道）
|
||||
res = []
|
||||
for lo, la in zip(path.lon, path.lat):
|
||||
dmin = min(hav((lo, la), (x, y)) for x, y in allp)
|
||||
res.append(dmin)
|
||||
path["dist_to_river_km"] = res
|
||||
path.to_csv(os.path.join(RES, "q4_path_profile_with_river.csv"), index=False, encoding="utf-8-sig")
|
||||
print(f"路径到最近河网节点：中位 {np.median(res):.3f} km，p90 {np.percentile(res,90):.3f} km，最大 {max(res):.3f} km")
|
||||
print("前 6 个路径点的最近河距 (km):", [round(x,3) for x in res[:6]])
|
||||
print("后 6 个路径点的最近河距 (km):", [round(x,3) for x in res[-6:]])
|
||||
@@ -0,0 +1,334 @@
|
||||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
MathorCup 2026 赛道B 初步试解（DeepSeek 子任务 B）
|
||||
AI 辅助说明：本脚本由 DeepSeek-V4.1-Flash（DeepSeek）于 2026-10-23 生成，
|
||||
用于竞赛初步评估；核心建模与结论需参赛队人工复核后方可采用。
|
||||
|
||||
做三件事（全部只读原始附件，结果写到本目录 结果/）：
|
||||
Q1 基线：候选气象/土壤/辐射特征在事件窗相对同季同时刻对照的百分位排序，并做按年 bootstrap 稳定性。
|
||||
Q4 基线：在真实 DEM 上从公开核定的源区出发做 D8 汇流路径 + Voellmy 型一维运动积分，
|
||||
输出路径长度/落差/到达时间/速度，并对摩擦参数与流深做敏感性。
|
||||
数据核查：附件冰川 NetCDF 的 UTM 坐标与本地 DEM、以及公开 AOI 坐标是否一致。
|
||||
|
||||
不做的：不把情景区间冒充 95% 置信区间；不把无记录日当负样本；不宣称预报能力。
|
||||
"""
|
||||
import os, re, json, struct
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
from PIL import Image
|
||||
|
||||
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) # B_DeepSeek试解
|
||||
ATT = os.path.abspath(os.path.join(ROOT, "..", "..", "附件")) # MathorCup2026/附件
|
||||
EXT = os.path.join(ROOT, "外部数据")
|
||||
RES = os.path.join(ROOT, "结果")
|
||||
os.makedirs(RES, exist_ok=True)
|
||||
|
||||
TMP = "/tmp/opencode/b_solve" # 临时解压目录
|
||||
os.makedirs(TMP, exist_ok=True)
|
||||
|
||||
# ---------------------------------------------------------------- 通用读取
|
||||
def parse_ts(s):
|
||||
s = str(s).strip().replace("/", "-")
|
||||
m = re.match(r"(\d{4})-(\d{1,2})-(\d{1,2})[ ](\d{1,2}):(\d{2})", s)
|
||||
if not m:
|
||||
return None
|
||||
y, mo, d, h, mi = m.groups()
|
||||
return f"{int(y):04d}-{int(mo):02d}-{int(d):02d}T{int(h):02d}:{mi}:00"
|
||||
|
||||
def load_table(kind, year):
|
||||
"""读取附件（1）解压后的某类表；9999.9 置为缺失。"""
|
||||
f = os.path.join(TMP, "met", f"{kind}_Gyirong_{year}.csv")
|
||||
df = pd.read_csv(f, encoding="utf-8-sig")
|
||||
df.columns = [c.strip().replace("\ufeff", "") for c in df.columns]
|
||||
t = np.array([parse_ts(x) for x in df["Timestamp"]], dtype="datetime64[m]")
|
||||
num = df.drop(columns=["Timestamp"]).apply(pd.to_numeric, errors="coerce")
|
||||
a = num.to_numpy(dtype=float)
|
||||
a[np.isclose(a, 9999.9)] = np.nan
|
||||
return t, pd.DataFrame(a, columns=num.columns)
|
||||
|
||||
def ensure_met():
|
||||
if os.path.isdir(os.path.join(TMP, "met")):
|
||||
return
|
||||
import zipfile
|
||||
z = [f for f in os.listdir(ATT) if f.startswith("（1）") and f.endswith(".zip")][0]
|
||||
with zipfile.ZipFile(os.path.join(ATT, z)) as zf:
|
||||
zf.extractall(os.path.join(TMP, "met"))
|
||||
|
||||
# ---------------------------------------------------------------- Q1 基线
|
||||
print("=" * 70)
|
||||
print("Q1 候选特征百分位基线")
|
||||
print("=" * 70)
|
||||
ensure_met()
|
||||
EV = np.datetime64("2026-08-26T10:00")
|
||||
EV_DOY = 238
|
||||
|
||||
FEATS = [
|
||||
("precip_1h", "mm/h"), ("precip_3h", "mm/3h"), ("precip_24h", "mm/24h"), ("precip_72h", "mm/72h"),
|
||||
("Ta_2m", "degC"), ("dT_2_10", "degC"), ("RH_2m", "%"), ("WS_10m", "m/s"),
|
||||
("shear_1_20", "1/s"), ("dpres_3h", "hPa/3h"), ("VWC_0.1", "m3/m3"),
|
||||
("dVWC_24h", "m3/m3/d"), ("soilT_0.1", "degC"), ("netrad", "W/m2"),
|
||||
]
|
||||
|
||||
def build_features(year):
|
||||
"""把 MET+SOIL+RADM 拼成小时特征表；全部按原始时间轴对齐。"""
|
||||
tm, m = load_table("MET", year)
|
||||
ts, s = load_table("SOIL", year)
|
||||
tr, r = load_table("RADM", year)
|
||||
assert (tm == ts).all() and (tm == tr).all()
|
||||
p = m["Precip(mm)"].to_numpy()
|
||||
out = {}
|
||||
out["precip_1h"] = p
|
||||
for k, w in [(3, 3), (24, 24), (72, 72)]:
|
||||
c = pd.Series(p).rolling(w, min_periods=w).sum().to_numpy()
|
||||
out[f"precip_{k}h"] = c
|
||||
out["Ta_2m"] = m["Ta_2m(degree C)"].to_numpy()
|
||||
out["dT_2_10"] = m["Ta_2m(degree C)"].to_numpy() - m["Ta_10m(degree C)"].to_numpy()
|
||||
out["RH_2m"] = m["RH_2m(%)"].to_numpy()
|
||||
out["WS_10m"] = m["WS_10m(m/s)"].to_numpy()
|
||||
out["shear_1_20"] = (m["WS_20m(m/s)"].to_numpy() - m["WS_1m(m/s)"].to_numpy()) / 19.0
|
||||
pres = pd.Series(m["Pressure(hPa)"].to_numpy())
|
||||
out["dpres_3h"] = pres.diff(3).to_numpy()
|
||||
vwc = s["VWC_0.1m(m3/m3)"].to_numpy()
|
||||
out["VWC_0.1"] = vwc
|
||||
out["dVWC_24h"] = pd.Series(vwc).diff(24).to_numpy()
|
||||
out["soilT_0.1"] = s["soilT_0.1m(degree C)"].to_numpy()
|
||||
out["netrad"] = (r["Rsd(W/m2)"].to_numpy() - r["Rsu(W/m2)"].to_numpy()
|
||||
+ r["Rld(W/m2)"].to_numpy() - r["Rlu(W/m2)"].to_numpy())
|
||||
return tm, out
|
||||
|
||||
data = {}
|
||||
for y in range(2021, 2027):
|
||||
data[y] = build_features(y)
|
||||
|
||||
def matched_values(year, name, doy=EV_DOY, hour=10, win=20):
|
||||
t, f = data[year]
|
||||
doy_v = (t.astype("datetime64[D]") - t.astype("datetime64[Y]").astype("datetime64[D]")).astype(int) + 1
|
||||
hh = t.astype("datetime64[h]").astype(int) % 24
|
||||
sel = (((np.abs(doy_v - doy) <= win) | (np.abs(doy_v - doy) >= 365 - win)) & (hh == hour))
|
||||
v = f[name][sel]
|
||||
return v[np.isfinite(v)]
|
||||
|
||||
def pct_rank(x, base):
|
||||
base = base[np.isfinite(base)]
|
||||
if not np.isfinite(x) or len(base) < 20:
|
||||
return np.nan
|
||||
return 100.0 * ((base < x).sum() + 0.5 * (base == x).sum()) / len(base)
|
||||
|
||||
rows = []
|
||||
for name, unit in FEATS:
|
||||
te, fe = data[2026]
|
||||
i = int(np.argmin(np.abs(te - EV)))
|
||||
x = fe[name][i]
|
||||
base = np.concatenate([matched_values(y, name) for y in range(2021, 2026)])
|
||||
p = pct_rank(x, base)
|
||||
rows.append(dict(feature=name, unit=unit, event_value=None if not np.isfinite(x) else float(x),
|
||||
pctl_vs_season=p, n_control=int(len(base)),
|
||||
extremeness=abs(p - 50) if np.isfinite(p) else np.nan))
|
||||
q1 = pd.DataFrame(rows).sort_values("extremeness", ascending=False)
|
||||
q1.to_csv(os.path.join(RES, "q1_feature_percentiles.csv"), index=False, encoding="utf-8-sig")
|
||||
print(q1.to_string(index=False, float_format=lambda v: f"{v:.3f}"))
|
||||
|
||||
# 按年 bootstrap：5 个对照年有放回重抽，看排序位次稳定性
|
||||
rng = np.random.default_rng(20261023)
|
||||
boot_rank = {f: [] for f, _ in FEATS}
|
||||
for _ in range(500):
|
||||
yrs = rng.choice([2021, 2022, 2023, 2024, 2025], size=5, replace=True)
|
||||
vals = {}
|
||||
for name, _unit in FEATS:
|
||||
te, fe = data[2026]
|
||||
i = int(np.argmin(np.abs(te - EV)))
|
||||
base = np.concatenate([matched_values(y, name) for y in yrs])
|
||||
p = pct_rank(fe[name][i], base)
|
||||
vals[name] = abs(p - 50) if np.isfinite(p) else -1
|
||||
order = sorted(vals, key=lambda k: -vals[k])
|
||||
for r, k in enumerate(order, 1):
|
||||
boot_rank[k].append(r)
|
||||
stab = pd.DataFrame([
|
||||
dict(feature=k, rank_mean=float(np.mean(v)), rank_p05=float(np.percentile(v, 5)),
|
||||
rank_p95=float(np.percentile(v, 95)), top5_freq=float(np.mean(np.array(v) <= 5)))
|
||||
for k, v in boot_rank.items()]).sort_values("rank_mean")
|
||||
stab.to_csv(os.path.join(RES, "q1_rank_stability.csv"), index=False, encoding="utf-8-sig")
|
||||
print("\n按年 bootstrap 排序稳定性（500 次，5 年有放回）：")
|
||||
print(stab.to_string(index=False, float_format=lambda v: f"{v:.2f}"))
|
||||
|
||||
# ---------------------------------------------------------------- 数据核查
|
||||
print("\n" + "=" * 70)
|
||||
print("数据核查：附件冰川 NetCDF 与本底 DEM / 公开 AOI 的坐标一致性")
|
||||
print("=" * 70)
|
||||
import zipfile as _zf
|
||||
if not os.path.isdir(os.path.join(TMP, "dem")):
|
||||
_dz = [f for f in os.listdir(ATT) if f.startswith("（7）") and f.endswith(".zip")][0]
|
||||
with _zf.ZipFile(os.path.join(ATT, _dz)) as z:
|
||||
z.extractall(os.path.join(TMP, "dem"))
|
||||
dem = np.array(Image.open(os.path.join(TMP, "dem", "DEM.tif"))).astype(float)
|
||||
dem[dem <= -32767] = np.nan
|
||||
nx, ny = dem.shape[1], dem.shape[0]
|
||||
px, py = 0.0002777467191455384, 0.00027774196409497727
|
||||
x0, y0 = 85.18975577464204, 28.564382284579324
|
||||
|
||||
import netCDF4 as nc
|
||||
gz = [f for f in os.listdir(ATT) if f.startswith("（4）") and f.endswith(".zip")][0]
|
||||
import zipfile
|
||||
if not os.path.isdir(os.path.join(TMP, "glacier")):
|
||||
with zipfile.ZipFile(os.path.join(ATT, gz)) as zf:
|
||||
zf.extractall(os.path.join(TMP, "glacier"))
|
||||
gfile = os.path.join(TMP, "glacier", [f for f in os.listdir(os.path.join(TMP, "glacier")) if f.endswith(".nc")][0])
|
||||
ds = nc.Dataset(gfile)
|
||||
gx = np.asarray(ds["x"][:]); gy = np.asarray(ds["y"][:]); gdem = np.asarray(ds["dem"][:], float)
|
||||
def utm_inv(E, N, lon0=87.0, k0=0.9996, E0=500000.0):
|
||||
a = 6378137.0; f = 1/298.257223563; e2 = 2*f - f*f; ep2 = e2/(1-e2)
|
||||
x = E - E0; M = N/k0
|
||||
mu = M/(a*(1 - e2/4 - 3*e2**2/64 - 5*e2**3/256))
|
||||
e1 = (1-np.sqrt(1-e2))/(1+np.sqrt(1-e2))
|
||||
p1 = (mu + (3*e1/2 - 27*e1**3/32)*np.sin(2*mu) + (21*e1**2/16 - 55*e1**4/32)*np.sin(4*mu)
|
||||
+ (151*e1**3/96)*np.sin(6*mu) + (1097*e1**4/512)*np.sin(8*mu))
|
||||
C1 = ep2*np.cos(p1)**2; T1 = np.tan(p1)**2
|
||||
N1 = a/np.sqrt(1-e2*np.sin(p1)**2); R1 = a*(1-e2)/(1-e2*np.sin(p1)**2)**1.5
|
||||
D = x/(N1*k0)
|
||||
lat = p1 - (N1*np.tan(p1)/R1)*(D**2/2 - (5+3*T1+10*C1-4*C1**2-9*ep2)*D**4/24
|
||||
+ (61+90*T1+298*C1+45*T1**2-252*ep2-3*C1**2)*D**6/720)
|
||||
lon = np.radians(lon0) + (D - (1+2*T1+C1)*D**3/6
|
||||
+ (5-2*C1+28*T1-3*C1**2+8*ep2+24*T1**2)*D**5/120)/np.cos(p1)
|
||||
return np.degrees(lon), np.degrees(lat)
|
||||
glo, gla = utm_inv(gx.mean(), gy.mean())
|
||||
aoi = json.load(open(os.path.join(EXT, "zenodo_22717809_gyirong_precursors", "detachment_aoi.geojson")))
|
||||
def flat(c):
|
||||
if isinstance(c[0], (int, float)):
|
||||
return [c]
|
||||
out = []
|
||||
for z in c:
|
||||
out += flat(z)
|
||||
return out
|
||||
pts = flat(aoi["features"][0]["geometry"]["coordinates"])
|
||||
alon = [p[0] for p in pts]; alat = [p[1] for p in pts]
|
||||
print(f"nc 质心 lon={glo:.5f} lat={gla:.5f}；公开 AOI bbox lon[{min(alon):.4f},{max(alon):.4f}] "
|
||||
f"lat[{min(alat):.4f},{max(alat):.4f}]；质心是否落在 AOI 内："
|
||||
f"{min(alon)<=glo<=max(alon) and min(alat)<=gla<=max(alat)}")
|
||||
# nc dem 与本地 DEM 重叠格点比对
|
||||
lons = gx; lats = utm_inv(gx[:, None]*np.ones((1, len(gy))), (gy[None, :]*np.ones((len(gx), 1))).T)[1] \
|
||||
if False else None
|
||||
# 逐格点反算（145x138=20010，可承受）
|
||||
xx, yy = np.meshgrid(gx, gy) # shape (138,145)
|
||||
llon, llat = utm_inv(xx, yy)
|
||||
col = ((llon - x0)/px).astype(int); row = ((y0 - llat)/py).astype(int)
|
||||
ok = (col >= 0) & (col < nx) & (row >= 0) & (row < ny) & np.isfinite(gdem)
|
||||
d_local = dem[row[ok], col[ok]]
|
||||
d_nc = gdem[ok]
|
||||
diff = d_nc - d_local
|
||||
print(f"重叠格点 {ok.sum()}；本地 DEM 均值 {np.nanmean(d_local):.0f} m，nc 均值 {np.nanmean(d_nc):.0f} m；"
|
||||
f"差值中位 {np.nanmedian(diff):+.1f} m，MAE {np.nanmean(np.abs(diff)):.1f} m")
|
||||
ds.close()
|
||||
|
||||
# ---------------------------------------------------------------- Q4 传播基线
|
||||
print("\n" + "=" * 70)
|
||||
print("Q4 情景传播基线：真实 DEM 上 D8 路径 + Voellmy 型一维运动积分")
|
||||
print("=" * 70)
|
||||
from scipy.ndimage import minimum_filter
|
||||
# 选路网格：4x4 块平均降到约 120 m 再 3x3 侵蚀填洼。
|
||||
# 30 m 单流向会卡在局部洼地（实测 10.6 km 处终止），120 m 网格可贯通到口岸河段。
|
||||
K = 4
|
||||
ny2, nx2 = ny//K, nx//K
|
||||
dem_route = dem[:ny2*K, :nx2*K].reshape(ny2, K, nx2, K).mean(axis=(1, 3))
|
||||
dem_route = minimum_filter(dem_route, size=3)
|
||||
dzdy, dzdx = np.gradient(dem, 30.7, 27.2)
|
||||
slope = np.degrees(np.arctan(np.hypot(dzdx, dzdy)))
|
||||
|
||||
src_lon, src_lat = float(np.mean(alon)), float(np.mean(alat))
|
||||
sc = int((src_lon - x0)/(px*K)); sr = int((y0 - src_lat)/(py*K))
|
||||
print(f"源区中心 (lon={src_lon:.4f}, lat={src_lat:.4f}) -> 行 {sr}, 列 {sc}, "
|
||||
f"DEM 高程 {dem[sr, sc]:.0f} m, 坡度 {slope[sr, sc]:.1f}°")
|
||||
|
||||
def lowest_neighbor(r, c, rad=3):
|
||||
"""D8 下沉；遇洼地在 rad 半径内找最低格点（简化 breach）。"""
|
||||
best = None
|
||||
for dr in range(-rad, rad+1):
|
||||
for dc in range(-rad, rad+1):
|
||||
if dr == 0 and dc == 0:
|
||||
continue
|
||||
rr, cc = r+dr, c+dc
|
||||
if not (0 <= rr < ny2 and 0 <= cc < nx2):
|
||||
continue
|
||||
if not np.isfinite(dem[rr, cc]):
|
||||
continue
|
||||
if dem_route[rr, cc] < dem_route[r, c] - 0.01:
|
||||
d = np.hypot(dr*30.7*K, dc*27.2*K)
|
||||
key = (dem_route[rr, cc], d)
|
||||
if best is None or key < best[0]:
|
||||
best = (key, rr, cc)
|
||||
return (best[1], best[2]) if best else (None, None)
|
||||
|
||||
path = [(sr, sc)]
|
||||
seen = {(sr, sc)}
|
||||
while len(path) < 20000:
|
||||
r, c = path[-1]
|
||||
if dem_route[r, c] < 1850: # 到达口岸河段高程量级
|
||||
break
|
||||
rr, cc = lowest_neighbor(r, c, rad=3)
|
||||
if rr is None or (rr, cc) in seen:
|
||||
# 允许一次半径外扩，仍失败则停
|
||||
rr, cc = lowest_neighbor(r, c, rad=6)
|
||||
if rr is None or (rr, cc) in seen:
|
||||
break
|
||||
seen.add((rr, cc))
|
||||
path.append((rr, cc))
|
||||
path = np.array(path)
|
||||
# 高程用真实 30 m DEM 采样（粗网格行/列映射回 30 m），沿路径 5 点中值平滑抑制 DEM 噪声
|
||||
zc = dem[path[:, 0]*K + K//2, path[:, 1]*K + K//2]
|
||||
prof_z = pd.Series(zc).rolling(5, center=True, min_periods=1).median().to_numpy()
|
||||
step = np.hypot(np.diff(path[:, 0])*30.7*K, np.diff(path[:, 1])*27.2*K)
|
||||
dist = np.concatenate([[0], np.cumsum(step)])/1000.0
|
||||
print(f"路径格点 {len(path)}，水平长度 {dist[-1]:.2f} km，落差 {prof_z[0]-prof_z[-1]:.0f} m，"
|
||||
f"末端高程 {prof_z[-1]:.0f} m")
|
||||
pd.DataFrame(dict(km=dist, elev_m=prof_z,
|
||||
lon=x0+path[:, 1]*px*K, lat=y0-path[:, 0]*py*K)).to_csv(
|
||||
os.path.join(RES, "q4_path_profile.csv"), index=False, encoding="utf-8-sig")
|
||||
|
||||
G = 9.81
|
||||
def voellmy(mu, xi, h, v0=0.0, z=None, s=None, dt=0.02):
|
||||
"""沿路径推进的 Voellmy 型一维运动（固定时间步 + 节点插值）。
|
||||
dv/dt = g sinθ - g μ cosθ - g v²/(ξ h)，θ 取当前节点到下一点的坡角。
|
||||
返回每个路径节点的到达时间 t 与速度 v。"""
|
||||
z = prof_z if z is None else z
|
||||
s = step if s is None else s
|
||||
L = float(np.sum(s))
|
||||
SN = np.concatenate([[0.0], np.cumsum(s)]) # 节点累计里程
|
||||
vt = np.zeros(len(SN)); vt[0] = v0
|
||||
ar = np.zeros(len(SN))
|
||||
v = max(v0, 0.1); t = 0.0; pos = 0.0; i = 0
|
||||
while pos < L and t < 20000:
|
||||
i = min(np.searchsorted(SN, pos, side="right") - 1, len(s)-1)
|
||||
th = np.arctan2(max(z[i]-z[i+1], 1e-6), s[i])
|
||||
a = G*np.sin(th) - G*mu*np.cos(th) - G*v*v/(abs(xi)*h)
|
||||
v = max(v + a*dt, 0.5) # 0.5 m/s 下限代表惯性，避免数值上“冻结”
|
||||
pos += v*dt
|
||||
t += dt
|
||||
while i+1 < len(SN) and SN[i+1] <= pos:
|
||||
vt[i+1] = v; ar[i+1] = t; i += 1
|
||||
arr = ar.copy()
|
||||
arr[arr == 0] = t
|
||||
return vt, arr
|
||||
|
||||
scen = []
|
||||
for mu in [0.05, 0.10, 0.20]:
|
||||
for xi in [500, 1000, 2000]:
|
||||
for h in [10.0, 25.0]:
|
||||
v, tt = voellmy(mu, xi, h)
|
||||
scen.append(dict(mu=mu, xi_m_s2=xi, h_m=h, v_max_m_s=float(np.nanmax(v)),
|
||||
t_total_s=float(tt[-1]), t_total_min=float(tt[-1]/60),
|
||||
v_end_m_s=float(v[-1]), runout_km=float(dist[-1])))
|
||||
sc = pd.DataFrame(scen).sort_values("t_total_s")
|
||||
sc.to_csv(os.path.join(RES, "q4_voellmy_sensitivity.csv"), index=False, encoding="utf-8-sig")
|
||||
print(sc.to_string(index=False, float_format=lambda v: f"{v:.1f}"))
|
||||
|
||||
# 中位参数的代表情景，逐段到达时间
|
||||
v_ref, t_ref = voellmy(0.10, 1000, 25.0)
|
||||
marks = []
|
||||
for k in np.arange(0, dist[-1]+0.001, 2.0):
|
||||
if k <= 0:
|
||||
continue
|
||||
j = int(np.argmin(np.abs(dist - k)))
|
||||
marks.append(dict(km=float(dist[j]), t_min=float(t_ref[j]/60), v_m_s=float(v_ref[j]), elev=float(prof_z[j])))
|
||||
pd.DataFrame(marks).to_csv(os.path.join(RES, "q4_arrival_marks.csv"), index=False, encoding="utf-8-sig")
|
||||
print("\n代表情景 (μ=0.10, ξ=1000, h=25 m) 每 2 km 到达时间：")
|
||||
print(pd.DataFrame(marks).to_string(index=False, float_format=lambda v: f"{v:.1f}"))
|
||||
Reference in new issue
Block a user