투영좌표 변환 (TM과 UTM)

TM(Transverse Mercator) 좌표를 변환하다가, 숫자가 계속 안맞는 안타까운 경험을 했습니다. 정말 오랜만에 TM 좌표를 써 본 것인지, 그 사이 우리 법이 바뀐 것을 몰랐습니다.  결국 검색을 하다보니 어이없게도 단 한문장을 제가 그동안 모르고 살았습니다.

측량, 수로업무 및 지적에 관한 법률(2009.12.14)에 따라 TM 투영원점의 가산수치는 600,000/200,000이 적용됩니다.

남북방향의 가산수치(FALSE NORTHING)을 500,000이라고 믿고 있던 저는, 제 기억을 확신하며 애꿎은 IDL 만 의심했네요. 국토지리정보원이 고시한 GNSS 지점 좌표가 이상하게 안맞는 겁니다. FALSE NORTHING이 600,000으로 개정된지 여러해가 지났건만, 다루던 자료들이 다들 예전 자료이다 보니 그동안 문제 없이 살아왔던 것 같습니다. 고생 좀 하던 문제가 해결되었으니, 관련된 글을 하나 쓸 텐데, 사실 IDL에서 투영좌표 변환을 하는 방법은, 문법 그 자체로는 매우 쉽습니다. 항상 Parameter 들이 문제인 거죠.

세세한 내용은 기억할 수도 없는 것이어서, 도움말을 활용해야 합니다. 다음의 몇 가지만 기억해 두십시오.

  • 투영법을 정의하는 다양한 파라메터들이 있습니다. IDL에서는 이런 파라메터들을 하나의 구조체로 관리할텐데, 이 구조체를 만들어 내는 함수가 MAP_PROJ_INIT() 입니다. 투영법을 정의하는 함수라고 생각하면 되겠습니다.
  • 아무런 투영법이 적용되지 않은 좌표는 경위도 좌표입니다. 이 좌표로부터 다른 투영좌표로 변환하는 것을 FORWARD, 반대로 투영좌표를 경위도 좌표로 되돌리는 것을 INVERSE라고 합니다. 그래서, 경위도를 투영 좌표로 변환할 때는 MAP_PROJ_FORWARD() 함수를, 투영 좌표를 경위도 좌표로 변환할 때에는 MAP_PROJ_INVERSE() 함수를 사용합니다.
  • 투영좌표 A를 투영좌표 B로 바로 변환하는 함수는 없습니다. 투영좌표 A를 경위도 좌표로, 그 경위도 좌표를 투영좌표 B로 변환해야 합니다.

시작해 볼까요.

국토지리정보원에서 고시한 GNSS 기준점 중 일부인 부산, 대구, 울산 지점의 평면직각좌표(TM 좌표)입니다.

우리나라의 TM 좌표가 서부, 중부, 동부 원점이 있는데, 위 세 지점은 모두 동부원점을 기준으로 합니다. 이는 위도 38도, 경도 129도를 원점으로 합니다. 서부 원점은 경도 125도, 중부원점은 경도 127도 점을 기준으로 합니다. 원점으로부터 북쪽으로 증가하는 미터(m) 단위의 값을 Northing, 동쪽으로 증가하는 미터 단위의 값을 Easting이라고 합니다. 원점을 기준으로 서쪽이나 남쪽으로 가게되면 값이 줄어드는데, 음수가 나오는 것을 피하기 위해(계산착오를 줄이기 위해) 애초에 큼직한 수를 원점좌표값으로 쓰는데 이것이 남북값은 600,000m(FALSE NORTHING)이고 동서값은 200,000m(FALSE EASTING) 입니다. 지구 타원체를 정의하는 여러 장반경/단반경 조합이 있는데, 요즘은 WGS84 또는 GRS80 (둘은 유사합니다) 을 주로 사용됩니다.

내용이 제법 많아서 다 읽게 되지 않는 MAP_PROJ_INIT 함수의 도움말에 나오는 내용 중, TRANSVERSE MERCATOR 투영법과 관련된 얘기는 이 정도 되는 것 같습니다. 이 내용들을 IDL에서 다음과 같이 씁니다.

국토지리정보원 고시좌표 기준이 이러한데, 예전에는 MERCATOR_SCALE이 0.9996 이었던 때도 있었고, 아주 오래 전 일본 측량기준을 끌어다 쓸 때, 동경 기준점 자체가 틀려먹은 것을 알게 되어 기준경도(CENTER_LONGITUDE)에 10.405초를 더해주기도 했고, BESSEL 타원체를 쓰기도 했습니다. 앞에 말씀드렸듯이 비교적 최근까지도 FALSE_NORTHING=500,000 으로 설정되어 있었구요. 이게 역사와 이력이 있는 좌표계라서, 저는 TM 좌표계만 등장하면 이거 저거 조합을 시도해 봅니다. 물론 정확히 명시되어 있는 기준이 있으면 간단하지요.

어쨌든 위 한줄로 우리나라의 동부원점 TM 투영 좌표계 정의가 끝났습니다. 그 결과는 이제 map_info_tm 변수(구조체)에 들어 있습니다.

이 구조체의 내용은 사실 알 필요는 없습니다. 이후 함수들에게 이 구조체를 넘겨 주기만 하면 됩니다.

TM 투영이 말 그대로 Transverse라서 동서방향을 Y, 남북방향을 X로 표기하는, 일반인들에게는 반대로 보이는 좌표 표기를 합니다. 그러다 보니, 위 함수에서 Y, X 의 순서로 인수를 넘겨주고 있습니다. 중요한 것은 MAP_STRUCTURE 키워드인데, 여기에 앞에서 MAP_PROJ_INIT 함수로 정의한 구조체를 넣어 줍니다. 위에 사용된 MAP_PROJ_INVERSE 한 줄의 의미는 이렇겠네요.

 주어진 Y, X 좌표는 map_info_tm 에 정의된 투영 좌표이니, 이를 경위도로 변환해 주세요.

MAP_PROJ_INIT가 수많은 투영법을 정의하는 옵션을 알려주기 위해 비교적 긴 도움말을 가지고 있었던 것에 비해, MAP_PROJ_INVERSE 함수의 도움말은 심플합니다. 이 기능이 거의 전부거든요.

참고로, 10진 체계의 경위도를 도분초의 60진 체계로 변환하고자 한다면, IDL Astronomy 라이브러리의 SIXTY() 함수를 이용하시면 편합니다.

MAP_PROJ_FORWARD() 도 한번 사용해 보기 위해, 이제 위경도를 UTM 좌표로 변환해 보겠습니다. UTM 투영법이 새로 등장했으니 정의해야 되겠죠.

Universal Transverse Mercator 라는 말 처럼 Universal 하게 사용하려고 만든 투영법이어서 그런지 앞의 TM에 비해 훨씬 심플합니다. 우리 나라는 52N Zone에 속해서 ZONE=52 라고 설정하면 거의 끝나는 거예요(혹시나 남쪽 Zone을 정의하려면 음수를 쓰면 됩니다).

결과의 첫 행 (부산지점)은 52N Zone의 원점을 기준으로 동쪽으로 506810.23m, 북쪽으로 3898988.9 m에 위치한 지점이란 의미가 되겠습니다.

재미삼아, 위경도 좌표를 다시 TM 좌표로 변환하면 같은 값이 나올까요?

다행히도 같은 값이 나오고 있네요.

MAP_PROJ_FORWARD와 MAP_PROJ_INVERSE는 인자로

  • X좌표 배열, Y좌표 배열 또는 경도좌표배열, 위도좌표배열을 따로 줄 수도 있고,   ;예) MAP_PROJ_INVERSE(Y, X)
  • [2, n] 크기의 두 컬럼 배열로 묶어서 줄 수도 있습니다.  ;예) MAP_PROJ_FORWARD(LL)
  • 두 함수 모두 리턴값은 위에 LL 변수와 같이 [2, n] 의 두 컬럼 배열입니다.

타원체를 정의하는 방법으로 ELLIPSOID 키워드를 사용하지 않고, SEMIMAJOR_AXIS, SEMIMINOR_AXIS 또는 SPHERE_RADIUS 키워드를 이용하여 직접 미터(m) 단위로 입력할 수 있습니다.