Aislo/common_util/common_util_contour_tiles.py

"""등고선 줌별 벡터 조각 — 웹 메르카토르 z/x/y MVT 를 공용 창고에(PLAN 615 · 요구 496 · 499).

조각 하나 = 전국 등고선(EPSG:5179)을 조각 상자(+ 여유 8 px)로 읽음 → 줌 간격 거름(SQL · 높이 색인)
→ 5179 에서 자르고 단순화(약 1 px) → 3857 → 조각 상자로 자름 → MVT(겹 `contour` · 속성 `h` 높이 ·
`k` 굵은 선 1/0 · 높이 하나 = 지물 하나). 512 px 조각 · 범위 4096.
창고 = `storage/_shared/contour_tiles/v<규칙 판>_<전국 파일 지문>/z/x/y.pbf` — 처음 볼 때
만들어 영구 · 빈 조각도 저장(다시 안 셈) · 잠금 · `.part` → 바꿈.
전국 파일 없는 PC = 있는 판 조각을 믿음 · 없으면 빈 조각.
화면 = `ui_template_gl_map.ts` `contourTileSource` ·
길 = `D03_PlanAnalysis_Router` `/plan/contour/{z}/{x}/{y}.pbf`.
"""

from __future__ import annotations

import hashlib
import math
from pathlib import Path

import numpy as np
from pyproj import Transformer

from common_util.common_util_terrain_store import NATIONAL_GPKG, _gpkg_stamp, _Lock
from config.config_system import STORAGE_BASE_DIR

TILES_ROOT = Path(STORAGE_BASE_DIR) / "_shared" / "contour_tiles"
RULE_VERSION = 1
MIN_ZOOM, MAX_ZOOM = 8, 15
LAYER = "contour"
EXTENT = 4096
TILE_PX = 512
BUFFER_PX = 8
# 줌별 담는 간격 · 굵은 선(m) — 화면 선 간격 약 6~12 px(분석 496 2장 「보통」)
STEP_M = {8: 1000, 9: 500, 10: 200, 11: 100, 12: 50, 13: 25, 14: 10, 15: 5}
BOLD_M = {8: 5000, 9: 2500, 10: 1000, 11: 500, 12: 250, 13: 25, 14: 25, 15: 25}
SIMPLIFY_PX = {14: 0.5, 15: 0.25}  # 그 밖 1 px
# 국내 경위도(밖 = 빈 조각 · 요청 안 함 — 화면 소스 `bounds` 와 같음)
BOUNDS_LONLAT = (124.0, 33.0, 132.0, 39.0)
WORLD_M = 2 * math.pi * 6378137.0
# 단순화 허용치 — 메르카토르 m/px × 위도 36.5° 줄임(실제 땅 m)
GROUND_SCALE = math.cos(math.radians(36.5))

_TO_MERC = Transformer.from_crs("EPSG:5179", "EPSG:3857", always_xy=True)
_FROM_MERC = Transformer.from_crs("EPSG:3857", "EPSG:5179", always_xy=True)


def tile_bounds(z: int, x: int, y: int) -> tuple[float, float, float, float]:
    """조각 상자(EPSG:3857 m · 서 · 남 · 동 · 북)."""
    size = WORLD_M / 2**z
    west = -WORLD_M / 2 + x * size
    north = WORLD_M / 2 - y * size
    return west, north - size, west + size, north


def in_bounds(z: int, x: int, y: int) -> bool:
    """조각이 국내 경위도 상자에 닿나."""
    lon0, lat0, lon1, lat1 = BOUNDS_LONLAT
    n = 2**z

    def tx(lon: float) -> int:
        return math.floor((lon + 180) / 360 * n)

    def ty(lat: float) -> int:
        r = math.radians(lat)
        return math.floor((1 - math.log(math.tan(r) + 1 / math.cos(r)) / math.pi) / 2 * n)

    return tx(lon0) <= x <= tx(lon1) and ty(lat1) <= y <= ty(lat0)


def _folder() -> Path | None:
    """지금 판 창고 — 전국 파일 지문 · 없으면 있는 판 중 마지막 · 그것도 없으면 None."""
    stamp = _gpkg_stamp()
    if stamp:
        tag = hashlib.sha1(stamp.encode()).hexdigest()[:8]
        return TILES_ROOT / f"v{RULE_VERSION}_{tag}"
    found = sorted(TILES_ROOT.glob(f"v{RULE_VERSION}_*"))
    return found[-1] if found else None


def _empty() -> bytes:
    import mapbox_vector_tile

    return mapbox_vector_tile.encode([{"name": LAYER, "features": []}])


def build_tile(z: int, x: int, y: int) -> bytes:
    """조각 하나 만들기(창고 안 씀) — 전국 파일에서 읽어 MVT."""
    import mapbox_vector_tile
    import pyogrio
    import shapely

    if not (MIN_ZOOM <= z <= MAX_ZOOM) or not in_bounds(z, x, y) or not NATIONAL_GPKG.is_file():
        return _empty()
    west, south, east, north = tile_bounds(z, x, y)
    pad = (east - west) * BUFFER_PX / TILE_PX
    box = (west - pad, south - pad, east + pad, north + pad)
    # 메르카토르 상자 둘레 → 5179 상자(가장자리 휨 감쌈)
    t = np.linspace(0.0, 1.0, 9)
    mx = np.concatenate(
        [box[0] + (box[2] - box[0]) * t] * 2 + [np.full(9, box[0]), np.full(9, box[2])]
    )
    my = np.concatenate(
        [np.full(9, box[1]), np.full(9, box[3])] + [box[1] + (box[3] - box[1]) * t] * 2
    )
    gx, gy = _FROM_MERC.transform(mx, my)
    grow = 10.0
    rect = (gx.min() - grow, gy.min() - grow, gx.max() + grow, gy.max() + grow)
    step, bold = STEP_M[z], BOLD_M[z]
    frame = pyogrio.read_dataframe(
        NATIONAL_GPKG,
        bbox=rect,
        columns=["CTRLN_HG"],
        where=f"CTRLN_HG % {step} = 0 OR CTRLN_HG % {bold} = 0",
    )
    if frame.empty:
        return _empty()
    tol = SIMPLIFY_PX.get(z, 1.0) * WORLD_M / (TILE_PX * 2**z) * GROUND_SCALE
    geoms = shapely.simplify(shapely.clip_by_rect(frame.geometry.values, *rect), tol)
    geoms = shapely.transform(
        geoms, lambda xy: np.column_stack(_TO_MERC.transform(xy[:, 0], xy[:, 1]))
    )
    geoms = shapely.clip_by_rect(geoms, *box)
    parts, owner = shapely.get_parts(geoms, return_index=True)
    keep = shapely.get_num_points(parts) >= 2
    parts, levels = parts[keep], frame["CTRLN_HG"].to_numpy(dtype=np.int64)[owner][keep]
    features = [
        {
            "geometry": shapely.multilinestrings(parts[levels == level]),
            "properties": {"h": int(level), "k": int(level % bold == 0)},
        }
        for level in np.unique(levels)
    ]
    return mapbox_vector_tile.encode(
        [{"name": LAYER, "features": features}],
        default_options={"quantize_bounds": (west, south, east, north), "extents": EXTENT},
    )


def ensure_tile(z: int, x: int, y: int) -> bytes:
    """조각 바이트 — 창고에 있으면 그대로 · 없으면 만들어 저장(잠금 · `.part` → 바꿈)."""
    if not (MIN_ZOOM <= z <= MAX_ZOOM) or not (0 <= x < 2**z and 0 <= y < 2**z):
        raise KeyError(f"{z}/{x}/{y}")
    if not in_bounds(z, x, y):
        return _empty()
    folder = _folder()
    if folder is None:
        return _empty()
    path = folder / str(z) / str(x) / f"{y}.pbf"
    if path.is_file():
        return path.read_bytes()
    if not NATIONAL_GPKG.is_file():
        return _empty()  # 전국 파일 없는 PC = 남은 판에 없으면 빈 조각(저장 안 함)
    with _Lock(path.with_suffix(".lock")):
        if path.is_file():
            return path.read_bytes()
        data = build_tile(z, x, y)
        part = path.with_suffix(".part")
        part.write_bytes(data)
        part.replace(path)
        return data