Aislo/B04_wf1_Surface/B04_wf1_Surface_Engine_MapSheet.py

# B04_wf1_Surface_Engine_MapSheet.py
# 1:5,000 수치지형도 도엽번호 산출 엔진.
# 국토지리정보원 도엽번호 체계:
#   - 경위도 1°x1° 구획: 위도 2자리 + 경도 끝 1자리 (예: 369 = 위도 36~37°, 경도 129~130°)
#   - 1:50,000 = 1° 구획을 4x4(15'x15') 분할, NW부터 행 우선 01~16 (예: 36906)
#   - 1:5,000  = 1:50,000 도엽을 10x10(1.5'x1.5') 분할, NW부터 행 우선 001~100
#     (예: 36906042 = 36906 도엽의 042 구획)

from __future__ import annotations

# 1:5,000 도엽 한 변 크기 (도 단위): 15' / 10 = 1.5'
SHEET5K_SIZE_DEG = 0.025
# 1:50,000 도엽 한 변 크기 (도 단위): 15'
SHEET50K_SIZE_DEG = 0.25


def _lon_digit_to_degree(digit: int) -> int:
    """경도 끝자리 → 경도 정수부 복원. 한반도 범위(124~132°E)에서 유일."""
    if digit >= 4:
        return 120 + digit  # 4~9 → 124~129
    return 130 + digit  # 0~2 → 130~132


def latlon_to_sheet5k(lat: float, lon: float) -> str:
    """WGS84 위경도 → 1:5,000 도엽번호(8자리 문자열)."""
    lat_deg = int(lat)
    lon_deg = int(lon)
    cell = f"{lat_deg:02d}{lon_deg % 10}"

    # 1° 구획 내 오프셋 (북서 원점 기준: 북쪽에서 남쪽으로 행 증가)
    lat_frac = (lat_deg + 1) - lat  # 구획 상단으로부터 남쪽 거리
    lon_frac = lon - lon_deg

    row50 = min(int(lat_frac / SHEET50K_SIZE_DEG), 3)
    col50 = min(int(lon_frac / SHEET50K_SIZE_DEG), 3)
    idx50 = row50 * 4 + col50 + 1  # 01~16

    # 1:50,000 도엽 내 오프셋
    lat_in_50 = lat_frac - row50 * SHEET50K_SIZE_DEG
    lon_in_50 = lon_frac - col50 * SHEET50K_SIZE_DEG
    row5 = min(int(lat_in_50 / SHEET5K_SIZE_DEG), 9)
    col5 = min(int(lon_in_50 / SHEET5K_SIZE_DEG), 9)
    idx5 = row5 * 10 + col5 + 1  # 001~100

    return f"{cell}{idx50:02d}{idx5:03d}"


def sheet5k_to_bounds(sheet_no: str) -> tuple[float, float, float, float]:
    """1:5,000 도엽번호 → WGS84 도곽 범위 (lon_min, lat_min, lon_max, lat_max)."""
    s = str(sheet_no).strip()
    if len(s) != 8 or not s.isdigit():
        raise ValueError(f"1:5,000 도엽번호는 8자리 숫자여야 함: {sheet_no!r}")

    lat_deg = int(s[0:2])
    lon_deg = _lon_digit_to_degree(int(s[2]))
    idx50 = int(s[3:5])
    idx5 = int(s[5:8])
    if not (1 <= idx50 <= 16):
        raise ValueError(f"1:50,000 구획번호(01~16) 범위 밖: {sheet_no!r}")
    if not (1 <= idx5 <= 100):
        raise ValueError(f"1:5,000 구획번호(001~100) 범위 밖: {sheet_no!r}")

    row50, col50 = divmod(idx50 - 1, 4)
    row5, col5 = divmod(idx5 - 1, 10)

    lat_max = (lat_deg + 1) - row50 * SHEET50K_SIZE_DEG - row5 * SHEET5K_SIZE_DEG
    lon_min = lon_deg + col50 * SHEET50K_SIZE_DEG + col5 * SHEET5K_SIZE_DEG
    return (lon_min, lat_max - SHEET5K_SIZE_DEG, lon_min + SHEET5K_SIZE_DEG, lat_max)


def sheet5k_center(sheet_no: str) -> tuple[float, float]:
    """도엽 중심 위경도 (lat, lon)."""
    lon_min, lat_min, lon_max, lat_max = sheet5k_to_bounds(sheet_no)
    return ((lat_min + lat_max) / 2.0, (lon_min + lon_max) / 2.0)


def neighbors_3x3(sheet_no: str) -> list[str]:
    """중심 도엽 + 인접 8매 = 3x3 도엽번호 목록 (NW부터 행 우선 9매).

    구획(1:50,000, 1°) 경계를 넘어가는 인접 도엽은 중심 좌표에
    도엽 크기만큼 오프셋을 더해 재산출하므로 번호 이월이 자동 처리된다.
    """
    lat_c, lon_c = sheet5k_center(sheet_no)
    result = []
    for dr in (1, 0, -1):  # 북 → 남
        for dc in (-1, 0, 1):  # 서 → 동
            result.append(
                latlon_to_sheet5k(lat_c + dr * SHEET5K_SIZE_DEG, lon_c + dc * SHEET5K_SIZE_DEG)
            )
    return result


def neighbors_for_points(points: list[tuple[float, float]]) -> list[str]:
    """기준 좌표들이 속한 도엽 + 각각의 주변 8매를 합친 목록(중복 제거, 순서 유지).

    계획노선 시점·종점이 같은 도엽이면 9매, 이웃한 두 도엽에 걸치면 12매가 된다
    (3×3 두 벌이 한 줄을 공유하므로 3×4). 도엽 하나가 늘 때마다 병합 산출물도 늘어나므로
    기준 좌표는 노선의 양 끝만 쓴다(2026-08-01 사용자 지시).
    """
    ordered: list[str] = []
    seen: set[str] = set()
    for lat, lon in points:
        for sheet_no in neighbors_3x3(latlon_to_sheet5k(lat, lon)):
            if sheet_no not in seen:
                seen.add(sheet_no)
                ordered.append(sheet_no)
    return ordered