Aislo/tools/precompute_terrain.py

# -*- coding: utf-8 -*-
"""전국 미리 계산 일괄 도구(PLAN 612-1 · 요구 493) — 공용 지형 창고에 도엽마다 미리 셈.

실행:  ./venv/Scripts/python.exe tools/precompute_terrain.py --sheets 36901 --workers 4
       --steps dem,zero,e07  단계(차례 고정 = 지형 → 0차곡 누적 → E07 · 뒤에 tiles)
       --sheets 36901,36712089  도엽 번호 또는 앞자리(1:50,000 묶음) · 주면 높이 차 거름 안 씀
       --box lon0,lat0,lon1,lat1  도엽 가운데가 이 경위도 상자 안
       --min-relief 150  도엽 안 높이 차(m) 이상만(산지) · --limit N 앞 N 장만
       --workers 4  일꾼 수 · --status 진행 기록 요약만 · --rebuild-index 도엽 목록 다시

- 도엽 목록 = 전국 등고선 파일 꼭짓점이 닿는 도엽(0.025°) — `_precompute_sheets.json`(전국 파일
  지문이 같으면 다시 안 셈) · 차례 = 등고선 길이 큰 것(산지)부터 · 같으면 높이 차.
- 좌표계 = 도엽 가운데 경도의 TM 띠(`tm_epsg_for_lon` · 5185~5187) · 격자 2 m.
- 셈 = 공용 함수만 부름(창고 지형 · D 0차곡 누적 · E07 도엽) — 저장 열쇠가 같으면 함수가 건너뜀 ·
  쓰기 = 함수 쪽 원자 쓰기(끊겨도 반쯤 쓴 파일 없음).
- 이어 하기 = `_precompute_progress.json`(단계 · 도엽마다 성공 · 시간) — 성공한 도엽은
  다시 안 부름 · 실패는 다음에 다시 · 전국 파일이 바뀌면 기록을 새로 시작.
- 시작 때 죽은 셈의 잠금(주인 프로세스 없음) · 오래된 임시 파일을 치움 ·
  낮은 우선순위(BELOW_NORMAL).
"""

from __future__ import annotations

import argparse
import json
import math
import os
import sqlite3
import sys
import time
from concurrent.futures import ProcessPoolExecutor, as_completed
from pathlib import Path
from typing import Any

ROOT = Path(__file__).resolve().parent.parent
if str(ROOT) not in sys.path:
    sys.path.insert(0, str(ROOT))

import numpy as np  # noqa: E402

from B04_PreProcess.B04_PreProcess_Engine_MapSheet import (  # noqa: E402
    SHEET5K_SIZE_DEG,
    latlon_to_sheet5k,
    neighbors_3x3,
    sheet5k_center,
)
from common_util import common_util_terrain_store as store  # noqa: E402
from common_util.common_util_crs import tm_epsg_for_lon  # noqa: E402
from common_util.common_util_json import atomic_write_json  # noqa: E402

STEPS = ("dem", "zero", "e07")
CELL_M = 2.0
PROGRESS_NAME = "_precompute_progress.json"
INDEX_NAME = "_precompute_sheets.json"
INDEX_CHUNK = 50_000  # 가닥 묶음(약 0.5 GB 읽기)
INDEX_STRIDE = 8  # 꼭짓점 솎음(도엽 2 km 에 비해 넉넉)
SAVE_EVERY_S = 10.0
# 단계별 창고 파일 앞머리(크기 셈)
STEP_FILES = {"dem": ("dem_",), "zero": ("zero_acc_",), "e07": ("e07_", "e07si_")}
_ENVELOPE = {0: 0, 1: 32, 2: 48, 3: 48, 4: 64}


def below_normal() -> None:
    """이 프로세스를 낮은 우선순위로(Windows) — 다른 일 · 서버를 막지 않게."""
    import psutil

    if os.name == "nt":
        psutil.Process().nice(psutil.BELOW_NORMAL_PRIORITY_CLASS)


def epsg_of(sheet_no: str) -> int:
    return tm_epsg_for_lon(sheet5k_center(sheet_no)[1])


# ── 도엽 목록 ────────────────────────────────────────────────────────────────


def aggregate(wkbs: list[bytes], heights: np.ndarray) -> dict[int, list[float]]:
    """EPSG:5179 선(WKB) · 높이 → {칸 열쇠: [길이 m, 최저, 최고]} · 칸 = 0.025°(위도 · 경도)."""
    import shapely
    from pyproj import Transformer

    geoms = shapely.from_wkb(wkbs)
    coords, owner = shapely.get_coordinates(geoms, return_index=True)
    coords, owner = coords[::INDEX_STRIDE], owner[::INDEX_STRIDE]
    if not len(coords):
        return {}
    lon, lat = Transformer.from_crs("EPSG:5179", "EPSG:4326", always_xy=True).transform(
        coords[:, 0], coords[:, 1]
    )
    keys = np.floor(lat / SHEET5K_SIZE_DEG).astype(np.int64) * 100_000 + np.floor(
        lon / SHEET5K_SIZE_DEG
    ).astype(np.int64)
    seg = np.zeros(len(coords))
    same = owner[1:] == owner[:-1]
    seg[:-1] = np.where(same, np.hypot(*(coords[1:] - coords[:-1]).T), 0.0)
    z = np.asarray(heights, dtype=np.float64)[owner]
    uniq, inv = np.unique(keys, return_inverse=True)
    length = np.bincount(inv, weights=seg, minlength=len(uniq))
    lo = np.full(len(uniq), np.inf)
    hi = np.full(len(uniq), -np.inf)
    np.minimum.at(lo, inv, z)
    np.maximum.at(hi, inv, z)
    return {int(k): [float(a), float(b), float(c)] for k, a, b, c in zip(uniq, length, lo, hi)}


def _merge(total: dict[int, list[float]], part: dict[int, list[float]]) -> None:
    for k, (length, lo, hi) in part.items():
        if k in total:
            t = total[k]
            t[0] += length
            t[1] = min(t[1], lo)
            t[2] = max(t[2], hi)
        else:
            total[k] = [length, lo, hi]


def _gpkg_rows(gpkg: Path):
    """전국 등고선 가닥 묶음 — (WKB 목록 · 높이) · 읽기 전용."""
    con = sqlite3.connect(f"file:{gpkg.resolve().as_posix()}?mode=ro", uri=True)
    try:
        last = -1
        while True:
            rows = con.execute(
                "SELECT fid, geom, CTRLN_HG FROM contours WHERE fid > ? AND geom IS NOT NULL"
                " ORDER BY fid LIMIT ?",
                (last, INDEX_CHUNK),
            ).fetchall()
            if not rows:
                return
            last = rows[-1][0]
            wkbs = [g[8 + _ENVELOPE[(g[3] >> 1) & 7] :] for _, g, _ in rows]
            yield wkbs, np.array([h if h is not None else np.nan for _, _, h in rows])
    finally:
        con.close()


def sheet_index(rebuild: bool = False) -> dict[str, list[float]]:
    """{도엽: [등고선 길이 km, 최저 m, 최고 m]} — 전국 파일 지문이 같으면 저장본."""
    path = store.TERRAIN_ROOT / INDEX_NAME
    stamp = store._gpkg_stamp()
    if not rebuild:
        try:
            saved = json.loads(path.read_text(encoding="utf-8"))
            if saved.get("gpkg") == stamp or stamp is None:
                return saved["sheets"]
        except (OSError, ValueError, KeyError):
            pass
    if stamp is None:
        raise SystemExit(f"전국 등고선 파일 없음: {store.NATIONAL_GPKG}")
    started = time.monotonic()
    cells: dict[int, list[float]] = {}
    for n, (wkbs, heights) in enumerate(_gpkg_rows(store.NATIONAL_GPKG), 1):
        _merge(cells, aggregate(wkbs, heights))
        elapsed = time.monotonic() - started
        print(f"도엽 목록 {n * INDEX_CHUNK:,} 가닥 · {elapsed:.0f} s", flush=True)
    step = SHEET5K_SIZE_DEG
    sheets = {
        latlon_to_sheet5k((k // 100_000 + 0.5) * step, (k % 100_000 + 0.5) * step): [
            round(v[0] / 1000, 2),
            v[1],
            v[2],
        ]
        for k, v in cells.items()
        if math.isfinite(v[1])
    }
    atomic_write_json(path, {"gpkg": stamp, "sheets": sheets})
    return sheets


def pick_sheets(
    index: dict[str, list[float]],
    *,
    sheets: list[str] | None = None,
    box: tuple[float, float, float, float] | None = None,
    min_relief: float = 150.0,
    limit: int | None = None,
) -> list[str]:
    """고른 도엽 — 등고선 길이 큰 것부터(같으면 높이 차) · 전국 파일이 안 덮는 도엽은 뺌."""
    if sheets:
        chosen = [no for no in index if any(no.startswith(p) for p in sheets)]
    else:
        chosen = [no for no in index if index[no][2] - index[no][1] >= min_relief]
    if box is not None:
        lon0, lat0, lon1, lat1 = box
        inside = []
        for no in chosen:
            lat, lon = sheet5k_center(no)
            if lon0 <= lon <= lon1 and lat0 <= lat <= lat1:
                inside.append(no)
        chosen = inside
    chosen.sort(key=lambda no: (-index[no][0], -(index[no][2] - index[no][1]), no))
    return chosen[:limit] if limit else chosen


# ── 셈(일꾼) ─────────────────────────────────────────────────────────────────


def _mtime(path: Path) -> int | None:
    try:
        return path.stat().st_mtime_ns
    except OSError:
        return None


def run_step(step: str, sheet_no: str) -> tuple[str, str, dict[str, Any]]:
    """도엽 하나 · 단계 하나 — 공용 함수만 부름 · fresh = 이번에 새로 셈."""
    epsg = epsg_of(sheet_no)
    started = time.monotonic()
    try:
        if step == "dem":
            fresh = store._stored_ok(sheet_no, epsg, CELL_M) is None
            empty = store.ensure_sheet_terrain(sheet_no, epsg, CELL_M) is None
        elif step == "zero":
            from D03_PlanAnalysis import D03_PlanAnalysis_ZeroOrder as zero

            for no in neighbors_3x3(sheet_no):  # 여유 200 m 이웃 지형(누적 열쇠가 나중에 안 바뀌게)
                store.ensure_sheet_terrain(no, epsg, CELL_M)
            target = zero._acc_names(sheet_no, epsg, CELL_M)[0]
            before = _mtime(target)
            empty = zero.sheet_accumulation(sheet_no, epsg, CELL_M) is None
            fresh = _mtime(target) != before
        elif step == "e07":
            from E07_StabArea import E07_StabArea_Router as e07_router
            from E07_StabArea import E07_StabArea_Sheet as e07

            paths = e07._paths(sheet_no, epsg, CELL_M)
            before = (_mtime(paths["flow"]), _mtime(paths["si"]))
            empty = e07.ensure_sheet(sheet_no, epsg, CELL_M, e07_router.default_params()) is None
            fresh = (_mtime(paths["flow"]), _mtime(paths["si"])) != before
        else:
            raise ValueError(f"모르는 단계: {step}")
    except Exception as exc:  # noqa: BLE001 — 한 장이 죽어도 나머지는 감 · 기록에 까닭
        return step, sheet_no, {"ok": False, "error": f"{type(exc).__name__}: {exc}"[:300]}
    seconds = round(time.monotonic() - started, 2)
    return step, sheet_no, {"ok": True, "s": seconds, "fresh": fresh, "empty": empty}


# ── 진행 기록 · 정리 ─────────────────────────────────────────────────────────


def load_progress() -> dict[str, Any]:
    path = store.TERRAIN_ROOT / PROGRESS_NAME
    stamp = store._gpkg_stamp()
    try:
        saved = json.loads(path.read_text(encoding="utf-8"))
    except (OSError, ValueError):
        saved = None
    if not saved or saved.get("gpkg") != stamp or saved.get("cell_m") != CELL_M:
        saved = {"gpkg": stamp, "cell_m": CELL_M, "steps": {}}
    for step in STEPS:
        saved["steps"].setdefault(step, {})
    return saved


def save_progress(progress: dict[str, Any]) -> None:
    atomic_write_json(store.TERRAIN_ROOT / PROGRESS_NAME, progress)


def _alive(pid: int) -> bool:
    import psutil

    return psutil.pid_exists(pid)


def clean_dead(root: Path | None = None) -> int:
    """주인이 죽은 잠금 · 오래된 임시 파일을 지움(끊긴 일괄 셈이 남긴 것) → 지운 수."""
    root = root or store.TERRAIN_ROOT
    removed = 0
    if not root.is_dir():
        return 0
    now = time.time()
    for folder in root.iterdir():
        if not folder.is_dir():
            continue
        for path in folder.glob(".*"):
            try:
                if path.suffix == ".lock":
                    pid = int(path.read_text(encoding="utf-8").split()[0])
                    if _alive(pid):
                        continue
                elif path.suffix != ".tmp" or now - path.stat().st_mtime < store.LOCK_STALE_S:
                    continue
                path.unlink(missing_ok=True)
                removed += 1
            except (OSError, ValueError, IndexError):
                continue
    return removed


def run(
    sheets: list[str],
    steps: list[str],
    workers: int,
    progress: dict[str, Any],
    say=print,
) -> dict[str, Any]:
    """단계 차례대로 · 단계 안은 일꾼 여럿 · 성공 기록된 도엽은 건너뜀 → 단계별 벽시계 초."""
    walls: dict[str, float] = {}
    for step in [s for s in STEPS if s in steps]:
        done = progress["steps"][step]
        todo = [no for no in sheets if not done.get(no, {}).get("ok")]
        say(f"[{step}] {len(todo)} 장 (건너뜀 {len(sheets) - len(todo)})")
        started = time.monotonic()
        saved_at = started
        count = 0

        def record(result: tuple[str, str, dict[str, Any]]) -> None:
            nonlocal count, saved_at
            _, no, info = result
            done[no] = info
            count += 1
            state = f"{info['s']} s" if info["ok"] else f"실패 {info['error']}"
            say(f"[{step}] {count}/{len(todo)} {no} {state}")
            if time.monotonic() - saved_at > SAVE_EVERY_S:
                save_progress(progress)
                saved_at = time.monotonic()

        try:
            if workers <= 1:
                for no in todo:
                    record(run_step(step, no))
            else:
                with ProcessPoolExecutor(workers, initializer=below_normal) as pool:
                    futures = [pool.submit(run_step, step, no) for no in todo]
                    try:
                        for future in as_completed(futures):
                            record(future.result())
                    except BaseException:
                        pool.shutdown(wait=False, cancel_futures=True)
                        raise
        finally:
            save_progress(progress)
        walls[step] = round(time.monotonic() - started, 1)
    return walls


def summary(sheets: list[str], progress: dict[str, Any]) -> dict[str, Any]:
    """단계별 — 성공 · 새로 셈 · 실패 · 새로 셈 장당 초(가운데 · 평균 · 최대) · 장당 크기 MB."""
    out: dict[str, Any] = {}
    for step in STEPS:
        done = progress["steps"][step]
        infos = [done[no] for no in sheets if no in done]
        failed = {no: done[no]["error"] for no in sheets if done.get(no, {}).get("ok") is False}
        fresh = [i["s"] for i in infos if i.get("ok") and i.get("fresh")]
        sizes = []
        for no in sheets:
            if not done.get(no, {}).get("ok"):
                continue
            stem = f"{epsg_of(no)}_{CELL_M:g}m"
            folder = store.TERRAIN_ROOT / no
            total = 0
            for head in STEP_FILES[step]:
                total += sum(p.stat().st_size for p in folder.glob(f"{head}{stem}*"))
            if total:
                sizes.append(total / 1e6)
        out[step] = {
            "ok": sum(1 for i in infos if i.get("ok")),
            "fresh": len(fresh),
            "failed": failed,
            "s_median": round(float(np.median(fresh)), 1) if fresh else None,
            "s_mean": round(float(np.mean(fresh)), 1) if fresh else None,
            "s_max": round(float(np.max(fresh)), 1) if fresh else None,
            "mb_mean": round(float(np.mean(sizes)), 2) if sizes else None,
            "mb_total": round(float(np.sum(sizes)), 1) if sizes else 0.0,
        }
    return out


def main(argv: list[str] | None = None) -> int:
    parser = argparse.ArgumentParser(description="공용 지형 창고 전국 미리 계산(PLAN 612)")
    parser.add_argument("--steps", default="dem,zero,e07")
    parser.add_argument("--sheets", default="")
    parser.add_argument("--box", default="")
    parser.add_argument("--min-relief", type=float, default=150.0)
    parser.add_argument("--limit", type=int, default=0)
    parser.add_argument("--workers", type=int, default=4)
    parser.add_argument("--status", action="store_true")
    parser.add_argument("--rebuild-index", action="store_true")
    args = parser.parse_args(argv)
    steps = [s.strip() for s in args.steps.split(",") if s.strip()]
    unknown = [s for s in steps if s not in STEPS]
    if unknown:
        parser.error(f"모르는 단계 {unknown} — 지금 되는 단계 {STEPS}(tiles 는 다음)")
    box = tuple(float(v) for v in args.box.split(",")) if args.box else None
    if box is not None and len(box) != 4:
        parser.error("--box = lon0,lat0,lon1,lat1")

    below_normal()
    index = sheet_index(args.rebuild_index)
    wanted = [s.strip() for s in args.sheets.split(",") if s.strip()]
    sheets = pick_sheets(
        index, sheets=wanted or None, box=box, min_relief=args.min_relief, limit=args.limit or None
    )
    print(f"도엽 {len(sheets)} 장 (목록 {len(index)} 장) · 단계 {steps} · 일꾼 {args.workers}")
    progress = load_progress()
    if not args.status:
        print(f"치운 잠금 · 임시 파일 {clean_dead()}")
        started = time.monotonic()
        walls = run(sheets, steps, args.workers, progress)
        print(f"벽시계 {walls} · 전체 {time.monotonic() - started:.1f} s")
    print(json.dumps(summary(sheets, progress), ensure_ascii=False, indent=1))
    return 0


if __name__ == "__main__":
    raise SystemExit(main())