동네 예보 좌표 변환용 C 코드 변환

공공 데이터(동네 예보) 개발 가이드에 보면 위/경도 좌표를 격자 좌표로 변환 하는 C 코드가 있다.
바이너리로 만들어 subprocess로 호출 해서 해도 되지만 OS 환경에 따라 빌드해야 하는 번거로움이 있다.
이에, 파이썬 코드로 변환했다.

변환 과정에서 주의 사항은 다음과 같다.

math.pow()의 경우 결과를 float 타입으로 변경하는 과정이 있고, 여기서 복소수 같은 복합 연산 결과를 만났을때, 에러를 일으킬 수 있다.
So, 기본적으로 제공하는 pow 함수를 사용하면, 결과를 그대로 넘겨 준다.

math.pow 사용 예제
In [7]: math.pow(-2.6050890646938014, 0.715566847180628)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-7-7fb3d71c4e22> in <module>()
----> 1 math.pow(-2.6050890646938014, 0.715566847180628)
ValueError: math domain error

pow 사용 예제
In [8]: pow(-2.6050890646938014, 0.715566847180628)
Out[8]: (-1.2432573697690525+1.5461879897834474j)

In [9]: type(pow(-2.6050890646938014, 0.715566847180628))
Out[9]: complex

* complex(복소수)의 경우, real(실수)만 뽑아서 사용하면 된다.

코드가 하는 역활에 대해서 찾아 봤다.

지도 투영법은 위경도로 이루어진 가상의 좌표를 평면상에 옮기는 방법이다.
지도 투영법 중에 람베르트 정각원추도법(람베르트正角圓錐圖法, Lambert conformal conic projection)
을 사용한다.

마지막으로 실제 변경된 코드이다.

격자 좌표에서 위/경도 좌표로 변환하는 역방향은 제외 했다.

from math import asin, sin, cos, tan, log

NX = 149  # X축 격자점 수
NY = 253  # Y축 격자점 수


def main(lon, lat):
    # 동네예보 지도 정보
    class map:
        pass

    map.Re = 6371.00877  # 지도반경
    map.grid = 5.0  # 격자간격 (km)
    map.slat1 = 30.0  # 표준위도 1
    map.slat2 = 60.0  # 표준위도 2
    map.olon = 126.0  # 기준점 경도
    map.olat = 38.0  # 기준점 위도
    map.xo = 210 / map.grid  # 기준점 X좌표
    map.yo = 675 / map.grid  # 기준점 Y좌표

    return map_conv(lon, lat, map)  # x, y


# 좌표 변환
def map_conv(lon, lat, map):
    _x, _y = lamcproj(lon, lat, map)
    _x = int(_x + 1.5)
    _y = int(_y + 1.5)
    return _x, _y


# Lambert Conformal Conic Projection
def lamcproj(lon, lat, map):
    PI = asin(1.0) * 2.0
    DEGRAD = PI / 180.0

    re = map.Re / map.grid
    slat1 = map.slat1 * DEGRAD
    slat2 = map.slat2 * DEGRAD
    olon = map.olon * DEGRAD
    olat = map.olat * DEGRAD

    sn = tan(PI * 0.25 + slat2 * 0.5) / tan(PI * 0.25 + slat1 * 0.5)
    sn = log(cos(slat1) / cos(slat2)) / log(sn)
    sf = tan(PI * 0.25 + slat1 * 0.5)
    sf = pow(sf, sn) * cos(slat1) / sn
    ro = tan(PI * 0.25 + olat * 0.5)
    ro = re * sf / pow(ro, sn)

    ra = tan(PI * 0.25 + lat * DEGRAD * 0.5)
    ra = re * sf / pow(ra, sn)
    theta = lon * DEGRAD - olon
    if theta > PI:
        theta -= 2.0 * PI
    if theta < -PI:
        theta += 2.0 * PI
    theta *= sn
    x = ((ra * sin(theta)) + map.xo).real
    y = ((ro - ra * cos(theta)) + map.yo).real
    return x, 

커버 사진