Aislo/resources/data_global_contours/build_contour_gpkg.py

"""국가기본도 등고선(TN_CTRLN) → 단일 GeoPackage 변환.

기존 convert_to_gpkg.py는 표고 컬럼 후보를 고정 목록으로만 찾아 속성을 통째로
버렸다. 이 스크립트는 geometry 외 모든 컬럼을 그대로 보존하고, 공간 인덱스와
조회용 속성 인덱스까지 만들어 크롭 단계에서 바로 쓸 수 있게 한다.

원본은 5m 간격 그대로 유지한다. 1m 보간은 프로젝트 크롭 시점에 국소 영역에서만
수행한다(전국 1m DEM은 계산량이 비현실적).

[사전 필수 패키지]
    pip install geopandas pyogrio shapely

사용법
    python build_contour_gpkg.py <입력> [-o national_contours.gpkg]
                                 [--layer contours] [--chunk 50000]
                                 [--append] [--limit N]

    <입력>   : .zip | .shp | SHP가 든 폴더
    --layer  : 저장할 레이어명. 나중에 능선을 --append 로 덧붙일 때 사용
    --append : 기존 gpkg를 지우지 않고 레이어를 추가
    --limit  : 소스당 읽을 최대 건수 (동작 확인용)
"""

import argparse
import glob
import os
import sqlite3
import sys
import time
import zipfile

try:
    import pyogrio
except ImportError:
    print("[오류] 필수 라이브러리가 누락되었습니다.")
    print(">>> pip install geopandas pyogrio shapely")
    sys.exit(1)

# 레이어별로 조회 인덱스를 걸 컬럼. 없는 컬럼은 조용히 건너뛴다.
INDEX_COLUMNS = {
    "contours": ("CTRLN_HG", "CTRLN_SE", "TPGRPH_SE", "CONT", "DIVI", "SCLS"),
    "ridges": ("SCLS", "SERV"),
    "streams": ("DIVI", "TYPE", "STAT"),
}


def collect_sources(path: str) -> list[str]:
    """입력 경로에서 읽을 SHP 목록을 만든다 (zip은 /vsizip/ 경로로 변환)."""
    path = os.path.abspath(path)
    if path.lower().endswith(".zip"):
        with zipfile.ZipFile(path) as zf:
            inner = sorted(n for n in zf.namelist() if n.lower().endswith(".shp"))
        base = "/vsizip/" + path.replace("\\", "/")
        return [f"{base}/{n}" for n in inner]
    if path.lower().endswith(".shp"):
        return [path]
    if os.path.isdir(path):
        found = sorted(glob.glob(os.path.join(path, "**", "*.shp"), recursive=True))
        if found:
            return found
        out: list[str] = []
        for z in sorted(glob.glob(os.path.join(path, "**", "*.zip"), recursive=True)):
            out.extend(collect_sources(z))
        return out
    return []


def write_source(
    src: str,
    out_path: str,
    layer: str,
    encoding: str,
    chunk: int,
    limit: int | None,
    first: bool,
) -> int:
    """SHP 하나를 청크 단위로 읽어 GeoPackage에 이어쓴다.

    전체를 한 번에 메모리에 올리면 1.8GB SHP에서 수 GB를 쓰게 되므로 나눠 읽는다.
    """
    info = pyogrio.read_info(src, encoding=encoding)
    total = info["features"] if limit is None else min(info["features"], limit)
    written = 0
    started = time.monotonic()

    while written < total:
        take = min(chunk, total - written)
        gdf = pyogrio.read_dataframe(
            src, encoding=encoding, skip_features=written, max_features=take
        )
        if not len(gdf):
            break
        # 첫 청크만 새 레이어를 만들고 이후는 append 한다.
        pyogrio.write_dataframe(
            gdf,
            out_path,
            layer=layer,
            driver="GPKG",
            append=not (first and written == 0),
            spatial_index=True,
            promote_to_multi=True,
        )
        written += len(gdf)
        print(
            "    %s  %s / %s  (%.0fs)"
            % (
                os.path.basename(src),
                format(written, ","),
                format(total, ","),
                time.monotonic() - started,
            ),
            flush=True,
        )
        if len(gdf) < take:
            break
    return written


def create_attribute_indexes(out_path: str, layer: str) -> list[str]:
    """GeoPackage는 SQLite이므로 조회용 속성 인덱스를 직접 만든다."""
    made: list[str] = []
    wanted = INDEX_COLUMNS.get(layer, ())
    if not wanted:
        return made
    with sqlite3.connect(out_path) as con:
        cols = {row[1] for row in con.execute(f'PRAGMA table_info("{layer}")')}
        for col in wanted:
            if col not in cols:
                continue
            name = f"idx_{layer}_{col.lower()}"
            con.execute(f'CREATE INDEX IF NOT EXISTS "{name}" ON "{layer}" ("{col}")')
            made.append(col)
        con.commit()
    return made


def summarize(out_path: str) -> None:
    """생성 결과를 요약 출력한다."""
    print("\n[결과]")
    for name, geom in pyogrio.list_layers(out_path):
        info = pyogrio.read_info(out_path, layer=name)
        print("  레이어 %-12s %s  %s건" % (name, geom, format(info["features"], ",")))
        print("     컬럼 %s" % list(info["fields"]))
        box = info.get("total_bounds")
        if box is not None:
            print("     범위 X %.0f~%.0f / Y %.0f~%.0f" % (box[0], box[2], box[1], box[3]))
    size = os.path.getsize(out_path) / (1024**3)
    print("  파일 %s  (%.2f GB)" % (out_path, size))


def main() -> None:
    parser = argparse.ArgumentParser(description="국가 등고선 원본 → GeoPackage 변환")
    parser.add_argument("input", help=".zip | .shp | SHP가 든 폴더")
    parser.add_argument("-o", "--output", default=None, help="출력 gpkg 경로")
    parser.add_argument("--layer", default="contours", help="레이어명 (기본 contours)")
    parser.add_argument("--encoding", default="cp949", help="DBF 인코딩 (기본 cp949)")
    parser.add_argument("--chunk", type=int, default=50000, help="청크 크기 (기본 50000)")
    parser.add_argument("--limit", type=int, default=None, help="소스당 최대 건수 (확인용)")
    parser.add_argument("--append", action="store_true", help="기존 gpkg에 레이어 추가")
    args = parser.parse_args()

    out_path = args.output or os.path.join(
        os.path.dirname(os.path.abspath(__file__)), "national_contours.gpkg"
    )

    sources = collect_sources(args.input)
    if not sources:
        print("[오류] 읽을 SHP를 찾지 못했습니다: %s" % args.input)
        sys.exit(1)

    if os.path.exists(out_path) and not args.append:
        print("[중단] 출력 파일이 이미 있습니다: %s" % out_path)
        print("       덮어쓰려면 파일을 먼저 옮기거나 --append 로 레이어를 추가하세요.")
        sys.exit(1)

    info = pyogrio.read_info(sources[0], encoding=args.encoding)
    print("=" * 74)
    print("입력    : %s" % args.input)
    print("소스    : SHP %d개" % len(sources))
    print("컬럼    : %s" % list(info["fields"]))
    print("좌표계  : %s" % info["crs"])
    print("출력    : %s  (레이어 %s)" % (out_path, args.layer))
    print("=" * 74)

    started = time.monotonic()
    total = 0
    for idx, src in enumerate(sources):
        print("\n[%d/%d] %s" % (idx + 1, len(sources), os.path.basename(src)))
        total += write_source(
            src,
            out_path,
            args.layer,
            args.encoding,
            args.chunk,
            args.limit,
            first=(idx == 0 and not args.append),
        )

    print("\n총 %s건 기록 (%.1f분)" % (format(total, ","), (time.monotonic() - started) / 60))

    made = create_attribute_indexes(out_path, args.layer)
    print("속성 인덱스: %s" % (", ".join(made) if made else "없음"))

    summarize(out_path)


if __name__ == "__main__":
    main()