MAP_PROJ_FORWARD, MAP_PROJ_INVERSE 함수의 이해와 활용

지난 회에서는 NG 체계의 MAP 함수에서 지원되는 MapForward, MapInverse 메서드들에 관하여 알아보았습니다. 이미 언급되었듯이 이 메서드들은 지도상에서 경위도(Longitude/Latitude) 좌표와 Cartesian(X/Y) 좌표 사이의 변환을 담당합니다. 오늘은 DG 체계의 지도에서 이와 같은 좌표 변환을 담당하는 MAP_PROJ_FORWARD, MAP_PROJ_INVERSE 함수들에 관하여 알아보고자 합니다. 전자의 경우는 경위도 좌표를 거리 기반의 Cartesian 좌표로 변환해주고, 후자의 경우는 반대로 Cartesian 좌표를 경위도 좌표로 변환하는 역할을 합니다. 설명을 위하여 먼저 지도 그림을 하나 그려놓고 시작을 해야 하므로, 지난 회와 마찬가지로 LCC(Lambert Conformal Conic) 투영법으로 지도를 먼저 그려 봅시다. 물론 오늘은 DG 체계에서 작업을 해야 하므로 이에 맞는 명령들이 사용되어야 합니다. 지도의 경위도 범위 및 전반적인 사항들은 지난 회에서 그렸던 지도와 거의 유사합니다.

DEVICE, DECOMPOSED=0

WINDOW, XSIZE=600, YSIZE=600, RETAIN=2

!P.BACKGROUND = 255

limit = [25, 115, 45, 135]

MAP_SET, 35, 125, /CONIC, LIMIT=limit, STANDARD_PARALLELS=[30, 60], $

  /GRID, COLOR=0, LABEL=0, LONDEL=1, LATDEL=1, XMARGIN=[4, 4], YMARGIN=[2, 2]

MAP_CONTINENTS, /HIRES, /COASTS, COLOR=0

MAP_GRID, /BOX_AXES, COLOR=0, LONS=[115:135], LATS=[25:45]

이 예제코드를 실행하면 다음과 같이 DG 체계를 기반으로 한 LCC 투영법의 지도가 표출됩니다. 경도의 범위는 동경 115~135도, 위도의 범위는 북위 25~45도이고, 중심의 경도 및 위도는 각각 125도와 35도입니다.

그러면 먼저 MAP_PROJ_FORWARD 함수를 사용하여 특정 지점의 경위도 좌표를 Cartesian 좌표로 변환해 봅시다. 그런데 이를 위해서는 사전에 반드시 필요한 하나의 단계가 있습니다. 그것은 MAP_PROJ_INIT라는 함수를 사용하여 지도의 투영법, 범위 등에 관한 정보를 따로생성해주는 작업입니다. 그 내용은 다음과 같습니다.

m = MAP_PROJ_INIT('Lambert Conformal Conic', LIMIT=limit, $

  CENTER_LONGITUDE=125, CENTER_LATITUDE=35, $

  STANDARD_PAR1=30, STANDARD_PAR2=60)

여기서 MAP_PROJ_INIT 함수의 내부에서 사용된 인수 및 키워드들을 보면 앞서 MAP_SET에서 설정했던 지도의 세부 정보들이 그대로 반영되어 있습니다. 즉 경위도 범위, 중심 경위도 좌표, Standard Paralle 등과 같은 세부 정보들이 동일하게 명시되어야 합니다. 이렇게 하면 좌변의 변수(여기서는 m)를 통하여 MAP_PROJ_INIT 함수에 의하여 생성된 지도 정보 전반이 담겨진 일종의 구조체(Structure)를 돌려받게 됩니다. 이 변수는 이후부터 MAP_PROJ_FORWARD 및 MAP_PROJ_INVERSE 함수의 MAP_STRUCTURE 키워드에 부여하여 사용하게 됩니다. 이제 MAP_PROJ_FORWARD 함수를 사용하여 다음과 같이 특정 지점의 경위도 좌표를 변환해 봅시다.

c = MAP_PROJ_FORWARD(125, 35, MAP_STRUCTURE=m)

PRINT, c

여기서는 앞서 MAP_PROJ_INIT 함수를 사용하여 얻은 변수 m을 MAP_STRUCTURE 키워드에 부여해야 한다는 것을 반드시 유념해야 합니다. 그리고 지도 중심점의 경위도인 125, 35를 인수로 넣어봤는데 PRINT로 출력된 값들은 다음과 같습니다.

   2.2805716e-09   9.3132257e-10

이 값들은 사실상 0, 0이라고보면 됩니다. 이번에는 다음과 같이 위도만 36도로 바꿔서 결과를 출력해 봅시다.

c = MAP_PROJ_FORWARD(125, 36, MAP_STRUCTURE=m)

PRINT, c

그러면 출력된 값들은 다음과같습니다.

   2.2460023e-09       108769.30

앞서 위도가 35도일 경우와 달리 Y좌표값이 꽤 큽니다. XY좌표값의 단위는 m(meter)입니다. 따라서 현재와 같은 지도상에서 경위도 좌표 (125, 35)와 (125, 36)인 지점들에 대한 XY 좌표값에 의하면 두 지점 사이의 거리가 약 108.8km에 해당된다는 의미입니다. 위도만 37도로 바꿔서 계산해보면 그 결과는 다음과 같습니다.

c = MAP_PROJ_FORWARD(125, 37, MAP_STRUCTURE=m)

PRINT, c

   2.2115220e-09       217258.58

지금까지 변환된 XY 좌표값들을 보면 X값은 거의 그대로인 채 Y값만 계속 변화하고 있습니다. 인수로 사용되는 값은 위와 같이 단일값으로 넣어도 되지만 배열의 형태로도 가능합니다. 다음과 같이 위도만 일정 간격으로 변화하는 5개의 지점들에 대한 경도 및 위도 값 배열을 만들어서 인수로 투입해 봅시다.

lons = [125, 125, 125, 125, 125]

lats = [35, 37, 39, 41, 43]

c = MAP_PROJ_FORWARD(lons, lats, MAP_STRUCTURE=m)

HELP, c

PRINT, c

이렇게 하면 결과 역시 배열의 형태로 얻어집니다. 출력된 내용을 보면 결과 배열인 c는 2x5의 형태를 가지며 그 값들은 다음과 같습니다.

C               DOUBLE    = Array[2, 5]

   2.2805716e-09   9.3132257e-10

   2.2115220e-09       217258.58

   2.1427901e-09       433517.40

   2.0742993e-09       649017.71

   2.0059722e-09       864003.04

이러한 결과를 그림을 표출하는 것도 가능합니다. 즉 다음과 같이 PLOTS 프로시저를 사용하여 각 지점들을 심볼로 표시하는 것입니다. 물론 여기서는 인수로 들어가는 좌표값들은 경위도 단위가 되어야 하기 때문에 앞서 우리가 정의했던 lons, lats 배열을 사용해야 합니다.

LOADCT, 34

PLOTSYM, 0, /FILL

PLOTS, lons, lats, PSYM=8, COLOR=240, SYMSIZE=1.5, /DATA

LOADCT, 0

여기서는 각 지점을 심볼로 표시하기 위하여 PLOTSYM이라는 프로시저를 사용하였는데요. 이 기능은 IDL에서 기본적으로 지원되는 기능이 아니라 외부 라이브러리에서 가져온 것입니다. DG 체계에서는 PLOT, PLOTS 등의 프로시저로 데이터 포인트를 심볼을 표시하기 위하여 PSYM 키워드를 많이 사용하는데, IDL 자체적으로 PSYM 키워드를 통하여 지원되는 심볼들의 종류가 제한적이고 색상 처리 등에 있어서여러가지 한계점들이 있기 때문에 PLOTSYM과 같은 외부 루틴을 가져와서 사용하기도 합니다. 이 루틴은 원래는 IDL Astro 라이브러리에 포함되어 있기 때문에 혹시 이 라이브러리를 설치해서 사용하는 경우에는 그냥 쓰면 됩니다. 하지만 이 라이브러리를 굳이 설치하지 않고 그냥 PLOTSYM 기능만 사용하고자 할 경우에는 구글에서 'idl plotsym'으로 검색해보면 코드를 받을 수 있는 링크가 바로 뜨므로 이를 통하여 받아서 사용해도 됩니다. 그리고 심볼의 색상을 파란색으로 처리하기 위하여 임시로 34번 컬러테이블을 불러서 그 중 48번 색상이 사용되도록 하였습니다. 이와 같은 과정에 의하여 표출된 결과 그림은 다음과 같습니다.

그러면 이번에는 다음과 같이 경도 좌표를 35로 고정한 채 위도 좌표만 변화하는 배열을 투입하여 결과를 얻어 봅시다.

lons = [125, 127, 129, 131, 133]

lats = [35, 35, 35, 35, 35]

c = MAP_PROJ_FORWARD(lons, lats, MAP_STRUCTURE=m)

HELP, c

PRINT, c

출력된 결과는 다음과 같습니다.

C               DOUBLE    = Array[2, 5]

   2.2805716e-09   9.3132257e-10

       179240.31       2238.9780

       358368.77       8954.5147

       537273.58       20142.419

       715843.11       35795.710

이 결과값들을 보면 경도가 커질 경우에는 당연히 X좌표의 값들이 그만큼 커지는 것을 볼 수 있습니다. 그런데 Y좌표의 값들을 보면 0으로 고정되지 않고 약간씩 커지는 것을 볼 수 있습니다. 그 이유는 현재와 같은 지도에서는 위도선(동일 위도를 따라가는 선)이 굴곡이 있으므로 XY좌표의 관점에서 보면 Y좌표가 일정할 수는 없기 때문입니다. 이러한 현상에 대해서는 잠시후 MapInverse 메서드를 사용한 결과에서 더 확실하게 알 수 있습니다. 앞서 했던 것처럼 5개의 지점들을 지도상에 표시해보면 그 결과는 다음 그림과 같습니다.

LOADCT, 34

PLOTSYM, 0, /FILL

PLOTS, lons, lats, PSYM=8, COLOR=48, SYMSIZE=1.5, /DATA

LOADCT, 0

이제 이번에는 MAP_PROJ_INVERSE 함수를 사용해 봅시다.이 함수를 사용하려면 MAP_PROJ_FORWARD와는 반대로 XY 좌표값들을 인수로 넣어줘야 합니다. 그래서 다음과 같이 X좌표 값들은 0으로 고정시키고 Y좌표 값들만 일정 간격으로 떨어진 5개의 지점들을 가정하였습니다.

xs = [0, 0, 0, 0, 0]

ys = [0, 2E5, 4E5, 6E5, 8E5]

c = MAP_PROJ_INVERSE(xs, ys, MAP_STRUCTURE=m)

HELP, c

PRINT, c

여기서는 Y방향으로 2E5(m) 즉 20만m 또는 200km 만큼 간격이 떨어지도록 하였습니다. 거리 단위로 값을 주었기 때문에 이러한 지점들의 경위도 좌표가 어떻게 나올지는 결과를 봐야 알 수 있습니다. 실제로 출력된 결과는 다음과 같습니다.

C               DOUBLE    = Array[2, 5]

       125.00000       35.000000

       125.00000       36.840758

       125.00000       38.689509

       125.00000       40.544578

       125.00000       42.404217

이 결과를 보면 위도 좌표 값들은 35도부터 시작해서 약 1.8도 정도의 간격으로 떨어져 있음을 알 수 있습니다. 물론 간격이 일정하지는 않은데, 이는 투영법의 특성으로 인한 것입니다. 이 5개의 지점들의 실제 위치를 표시하려면 앞서와 같이 PLOTSYM, PLOTS 프로시저들을 사용하면 됩니다. 다만 이 프로시저들을 사용할 때에는 좌표를 경위도 단위로 줘야 하기 때문에, 이번에는 다음과 같이 결과 배열인 c를 적절히 인덱싱하여 PLOTS에 인수로 투입해야 합니다.

LOADCT, 34

PLOTSYM, 0, /FILL

PLOTS, c[0, *], c[1, *], PSYM=8, COLOR=240, SYMSIZE=1.5, /DATA

LOADCT, 0

그리고 심볼 색상은 붉은색으로 처리하였습니다. 이렇게 표출된 그림은 다음과 같습니다.

이번에는 다음과 같이 Y좌표를 0으로 고정한 채 X좌표만 일정 간격으로 떨어진 5개의 지점들을 가정하여 결과를 얻어봅시다.

xs = [0, 2E5, 4E5, 6E5, 8E5]

ys = [0, 0, 0, 0, 0]

c = MAP_PROJ_INVERSE(xs, ys, MAP_STRUCTURE=m)

HELP, c

PRINT, c

먼저 출력된 결과를 보면 다음과 같습니다.

C               DOUBLE    = Array[2, 5]

       125.00000       35.000000

       127.23083       34.974416

       129.45820       34.897735

       131.67869       34.770169

       133.88893       34.592069

이 결과를 보면 5개의 지점들의 경도 값들이 증가하는 것은 충분히 이해가 됩니다. 그런데 위도 값들을 보면 처음엔 35도였다가 점점 조금씩 낮아지는 것을 볼 수 있습니다. 이 역시 앞서 언급했던 투영법의 특성으로 인한 효과입니다. 5개의 지점들을 지도상에 표시한 그림을 그려보면 그 이유를 쉽게 이해할 수 있습니다.

LOADCT, 34

PLOTSYM, 0, /FILL

PLOTS, c[0, *], c[1, *], PSYM=8, COLOR=240, SYMSIZE=1.5, /DATA

LOADCT, 0

이 그림에서 보는 것처럼 XY 좌표계 기준으로는 Y좌표가 고정된 채 X방향으로 일정 간격으로 떨어진 지점들이지만, 경위도 좌표의 관점에서 보면 경위도선의 굴곡 또는 휘어짐에 의하여 X좌표가 커질수록 위도 좌표는 조금씩 낮아진다는 것을 알 수 있습니다. 따라서 거리 기반의XY좌표가 경위도 좌표로 어떻게 변환될 것인가는 당연히 지도 자체의 투영법의 특성에 따라 달라질 수 밖에 없습니다.

오늘 소개한 지도 표출 작업은 DG 체계 기반의 관련 기능들을 사용하여 처리하고 있는데, 맨 앞부분에서 그래픽창의 배경색이 기본적으로 흰색이 되도록 하기 위하여 !P.BACKGROUND라는 시스템 변수의 값을 255로 설정했었습니다. 이제 모든 작업들이 끝났으므로 이 값을 다음과 같이 다시 0으로 되돌려놓고 마무리하겠습니다.

!P.BACKGROUND = 0

오늘은 DG 체계에서 MAP_PROJ_FORWARD 및 MAP_PROJ_INVERSE 함수들을 활용하여 경위도 좌표와 Cartesian 좌표 사이의 변환을 하는 방법을 살펴보았습니다. 앞서 언급했듯이 이 함수들을 제대로 활용하기 위해서는 MAP_PROJ_INIT 함수를 사용하여 지도에 관한 세부적인 정보들(투영법, 범위, 중심좌표 등)이 담긴 정보 구조체를 생성해야 합니다. 사실 이러한 정보들은 지도의 표출을 위하여 MAP_SET 프로시저를 사용할 때 이미 명시가 됩니다. 하지만 좌표 변환을 위해서는 이러한 정보들을 MAP_PROJ_INIT 함수로도 따로 만들어줘야 합니다. DG 체계의 지도 기반으로 작업을 할 때에는 이렇습니다. 이러한 점이 지난 회 게시물에서 소개한 NG 체계 지도의 경우보다는 좀 번거롭게 느껴질 수도 있습니다. 그리고 제가 소개해드린 예제는 LCC 투영법의 지도를 예시로 다루었지만, 다른 투영법의 지도일 경우에는 지도에관한 세부 정보들을 MAP_SET 및 MAP_PROJ_INIT에서 정의하는 방식도 조금씩 다릅니다. 따라서 이러한 부분들에 관해서는 IDL 도움말을참조하여 실제 작업에 반영해야 한다는 점을 염두에 두어야 하겠습니다.