Aislo/resources/tester/test_589_tm_grid.py

# -*- coding: utf-8 -*-
"""589-1 B04 · B05 지도 — 경위도 겹(도엽 · 지적 · 유역 · 위성 그림)을 TM 겹(노선 · 관) 자리에.

- 서버 `tm_grid_for_meta` 마디 = pyproj · 왼위 · 오른아래 마디 = 메타 귀퉁이 1 mm 안(Bessel · GRS80)
- 메타 귀퉁이와 다른 좌표계면 격자 안 줌
- TS `normOfLonLat`(격자 쌍선형) = pyproj TM 자리 0.02 m 안 · 격자 없으면 옛 같은 네모 식 그대로
- 위성 그림 행렬(화면 가운데 1차) = 그림 px → 경위도 → pyproj TM 자리와 화면 안 0.5 px 안
- 옛 같은 네모 꼴은 귀퉁이에서 용화 ≥ 3 m · D ≥ 60 m 어긋남(고친 까닭)

Node 헬퍼(`helper_589_tm_grid.cjs`)가 진짜 TS 를 돌린다.
"""

import json
import math
import os
import subprocess

import pytest
from pyproj import Transformer

from B04_PreProcess.B04_PreProcess_Engine_VWorld import tm_grid_for_meta

HERE = os.path.dirname(os.path.abspath(__file__))

BESSEL = (
    'PROJCS["Korean 1985 / Modified East Belt",GEOGCS["Korean 1985",DATUM["Korean Datum 1985",'
    'SPHEROID["Bessel 1841",6377397.155,299.1528128],TOWGS84[-114.61998513841,475.96296986986,'
    '675.01832868573,1.162,-2.347,-1.592,6.342]],PRIMEM["Greenwich",0],UNIT["degree",'
    '0.01745329251994328]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",38],'
    'PARAMETER["central_meridian",129.002890277778],PARAMETER["scale_factor",1],'
    'PARAMETER["false_easting",200000],PARAMETER["false_northing",500000],UNIT["metre",1]]'
)
# 검증 B(용화 · Bessel) · D(5186) 위성 상자(lon_min, lon_max, lat_min, lat_max) · 그림 px
SITES = {
    "B": (
        BESSEL,
        (129.0728759765625, 129.100341796875, 36.798288873837045, 36.82687474287728),
        (2560, 3328),
    ),
    "D": (
        "EPSG:5186",
        (127.9193115234375, 127.979736328125, 36.270850207239036, 36.30184530368433),
        (2816, 1792),
    ),
}
VIEWS = [(1400, 560, 1, 0, 0), (1400, 560, 2.33, -120, 80), (1400, 560, 300, -60000, 40000)]


def _meta(crs, box):
    """서버 식 그대로 — 왼위 · 오른아래 두 귀퉁이만 TM 으로."""
    lon_min, lon_max, lat_min, lat_max = box
    fwd = Transformer.from_crs("EPSG:4326", crs, always_xy=True)
    x_min, y_max = fwd.transform(lon_min, lat_max)
    x_max, y_min = fwd.transform(lon_max, lat_min)
    return {
        "x_min": x_min,
        "x_max": x_max,
        "y_min": y_min,
        "y_max": y_max,
        "width_meters": x_max - x_min,
        "height_meters": y_max - y_min,
        "center_x": (x_min + x_max) / 2,
        "center_y": (y_min + y_max) / 2,
        "lon_min": lon_min,
        "lon_max": lon_max,
        "lat_min": lat_min,
        "lat_max": lat_max,
    }


def _merc(lat):
    return math.log(math.tan(math.pi / 4 + math.radians(lat) / 2))


def _lonlat_of_px(meta, size, u, v):
    """그림 px → 경위도(가로 경도 선형 · 세로 메르카토르 y 선형)."""
    top, bottom = _merc(meta["lat_max"]), _merc(meta["lat_min"])
    lon = meta["lon_min"] + u / size[0] * (meta["lon_max"] - meta["lon_min"])
    m = top - v / size[1] * (top - bottom)
    return lon, math.degrees(2 * math.atan(math.exp(m)) - math.pi / 2)


def _run(site, with_grid=True):
    crs, box, size = SITES[site]
    meta = _meta(crs, box)
    if with_grid:
        meta["tm_grid"] = tm_grid_for_meta(meta, crs)
    lon_min, lon_max, lat_min, lat_max = box
    lonlat = [
        (lon_min + a / 10 * (lon_max - lon_min), lat_min + b / 10 * (lat_max - lat_min))
        for a in range(11)
        for b in range(11)
    ]
    payload = {"meta": meta, "lonlat": lonlat, "image": size, "views": VIEWS}
    proc = subprocess.run(
        ["node", os.path.join(HERE, "helper_589_tm_grid.cjs")],
        input=json.dumps(payload),
        capture_output=True,
        text=True,
        encoding="utf-8",
        timeout=60,
    )
    assert proc.returncode == 0, proc.stderr
    return meta, lonlat, json.loads(proc.stdout)


@pytest.fixture(scope="module", params=["B", "D"])
def site(request):
    return request.param, *_run(request.param)


def test_서버_격자는_pyproj_마디(site):
    name, meta, _, _ = site
    grid = meta["tm_grid"]
    n = grid["n"]
    assert n == 9 and len(grid["x"]) == n * n
    assert abs(grid["x"][n * (n - 1)] - meta["x_min"]) < 1e-3
    assert abs(grid["y"][n * (n - 1)] - meta["y_max"]) < 1e-3
    assert abs(grid["x"][n - 1] - meta["x_max"]) < 1e-3
    assert abs(grid["y"][n - 1] - meta["y_min"]) < 1e-3
    fwd = Transformer.from_crs("EPSG:4326", SITES[name][0], always_xy=True)
    i, j = 3, 5
    lon = meta["lon_min"] + i / (n - 1) * (meta["lon_max"] - meta["lon_min"])
    lat = meta["lat_min"] + j / (n - 1) * (meta["lat_max"] - meta["lat_min"])
    x, y = fwd.transform(lon, lat)
    assert math.hypot(grid["x"][j * n + i] - x, grid["y"][j * n + i] - y) < 1e-3


def test_다른_좌표계_메타면_격자_안_줌():
    meta = _meta("EPSG:5186", SITES["D"][1])
    assert tm_grid_for_meta(meta, "EPSG:5187") is None


def test_TS_격자_자리는_pyproj_TM(site):
    name, meta, lonlat, out = site
    assert out["hasToNorm"]
    fwd = Transformer.from_crs("EPSG:4326", SITES[name][0], always_xy=True)
    worst = 0.0
    for (lon, lat), (nx, ny) in zip(lonlat, out["tm"]):
        x, y = fwd.transform(lon, lat)
        gx = meta["x_min"] + nx * meta["width_meters"]
        gy = meta["y_min"] + (1 - ny) * meta["height_meters"]
        worst = max(worst, math.hypot(gx - x, gy - y))
    assert worst < 0.02, worst


def test_격자_없으면_옛_식_그대로():
    meta, lonlat, out = _run("B", with_grid=False)
    assert not out["hasToNorm"]
    for (lon, lat), tm, old, screen in zip(lonlat, out["tm"], out["old"], out["oldScreen"]):
        want = [
            (lon - meta["lon_min"]) / (meta["lon_max"] - meta["lon_min"]),
            1 - (lat - meta["lat_min"]) / (meta["lat_max"] - meta["lat_min"]),
        ]
        assert tm == want and old == want and screen == want


def test_위성_그림_행렬은_TM_자리(site):
    name, meta, _, out = site
    crs, _, size = SITES[name]
    fwd = Transformer.from_crs("EPSG:4326", crs, always_xy=True)
    for (w, h, scale, ox, oy), got in zip(VIEWS, out["views"]):
        rect = got["mapRect"]
        a, b, c, d, e, f = got["matrix"]
        ax, ay = rect["width"] * scale, rect["height"] * scale
        bx = rect["x"] * scale + w / 2 * (1 - scale) + ox
        by = rect["y"] * scale + h / 2 * (1 - scale) + oy
        if scale < 10:
            # 그림 귀퉁이 · 가장자리 가운데 · 가운데
            probes = [(u, v) for u in (0, size[0] / 2, size[0]) for v in (0, size[1] / 2, size[1])]
        else:
            # 크게 당김 — 화면 귀퉁이 · 가운데에 놓인 그림 px
            det = a * d - b * c
            probes = []
            for sx, sy in ((0, 0), (w, 0), (0, h), (w, h), (w / 2, h / 2)):
                du, dv = sx - e, sy - f
                probes.append(((d * du - c * dv) / det, (a * dv - b * du) / det))
        for u, v in probes:
            x, y = fwd.transform(*_lonlat_of_px(meta, size, u, v))
            want = (
                (x - meta["x_min"]) / meta["width_meters"] * ax + bx,
                (1 - (y - meta["y_min"]) / meta["height_meters"]) * ay + by,
            )
            got_px = (a * u + c * v + e, b * u + d * v + f)
            gap = math.hypot(got_px[0] - want[0], got_px[1] - want[1])
            assert gap < 0.5, (name, scale, u, v, gap)


@pytest.mark.parametrize("name,least", [("B", 3.0), ("D", 60.0)])
def test_옛_같은_네모_꼴은_귀퉁이에서_어긋남(name, least):
    crs, box, _ = SITES[name]
    meta = _meta(crs, box)
    x, y = Transformer.from_crs("EPSG:4326", crs, always_xy=True).transform(box[1], box[3])
    # 옛 식 = 오른위 귀퉁이 경위도를 상자 오른위(x_max, y_max) 에 둠
    assert math.hypot(meta["x_max"] - x, meta["y_max"] - y) >= least