구면상의 지점 분포 데이터에 대한 격자화 및 유의사항

오늘은 글로벌 맵(Global Map) 즉 세계 전체 지도 상에서 분포하는 지점별 데이터들을 규칙 격자 데이터로 변환하는 작업에 대하여 소개하고 이와 관련된 유의사항에 관해서도 함께 언급해보고자 합니다. 사실 불규칙하게 지점별로 분포하는 데이터를 규칙 격자 데이터로 변환하는 방법에 관해서는 제가 이 블로그에서 여러 차례 관련 게시물들을 통하여 소개한 바 있습니다. 이 게시물들에서는 격자화를 위한 내삽 연산을 담당하는  GRIDDATA, KRIG2D, SPH_SCAT 등의 함수들의 사용법 및 관련 예제들이 상세히 설명되어 있습니다. 그런데 이 게시물들에서 주로 예제로 사용되었던 지점 분포 데이터들의 경우 그 지점들이 주로 한반도 및 그 주변 영역에 제한되었습니다. 아무래도 국내 사용자들 입장에서는 그런 데이터를 사용하게 될 가능성이 높다고 판단했기 때문에 저도 한반도 및 주변에만 분포하는 지점들로 가정하여 예제를 준비하여 소개했던 것입니다.

 

그런데 연구 분야나 작업의 특성에 따라서는 한반도와 같은 국지적인 영역에 국한되지 않고 전 세계에 걸쳐 퍼져있는 지점 분포 데이터를 사용하게 되는 경우도 꽤 자주 있습니다. 그래서 오늘은 이러한 글로벌 지점 분포 데이터를 예제로 사용해보고자 합니다. 그리고 글로벌 지점 분포 데이터를 규칙 격자 데이터로 변환하는 작업에서 특별히 유의해야 할 사항도 있습니다. 이 부분은 글로벌 분포 지점들이라는 것이 결국 지구라는 구체(Sphere)의 표면 상에 존재하는 지점들이기 때문에 고려되어야 하는 기하학적 특성과 관련이 있는데, 함께 언급하도록 하겠습니다. 오늘 소개될 내용은 대략 다음과 같은 사항들을 따를 것입니다.

 

1. 글로벌 맵 상에 분포하는 것으로 가정하여 생성된 가상의 지점별 데이터 사용

2. 격자화 변환을 위해서는 GRIDDATA 함수를 사용하고 세부적으로는 Kriging 기법 적용

 

그러면 먼저 예제로 사용할 가상 데이터를 생성합니다. 다음과 같은 내용을 실행하여 data, lons, lats 배열을 생성하면 됩니다.

 

n = 300

data = RANDOMU(-1, n, /NORMAL)*12+50

lons = RANDOMU(-3, n)*360-180

lats = RANDOMU(-5, n)*180-90

HELP, lons, lats, data

PRINT, MIN(data), MAX(data)

 

여기서 생성된 data, lons, lats는 각각 300개의 값들로 구성됩니다. 즉 전 세계에 걸쳐 분포하는 300개의 지점들에 대한 데이터로 가정한 것입니다. 지점들이 전 세계에 분포하도록 하기 위하여 lons 배열의 값들은 경도의 전체 범위인 -180~+180의 범위 내에 있도록 하였고 lats 배열의 값들은 위도의 전체 범위인 -90~+90의 범위 내에 있도록 하였습니다. 그리고 지점별 값에 해당되는 data 배열의 값들은 대략 0~100의 범위 내에 있도록 하였습니다. 이를 위하여 RANDOMU 함수를 사용하여 난수를 발생시키고 이러한 값들이 각각의 범위에 들어가도록 하기 위한 약간의 연산이 수행되도록 하였습니다. HELP와 PRINT로 출력된 내용을 보면 의도에 맞게 배열들이 생성되었음을 확인할 수 있습니다. 출력된 내용은 다음과 같습니다.

 

LONS            FLOAT     = Array[300]

LATS            FLOAT     = Array[300]

DATA            FLOAT     = Array[300]

      12.3302      81.3593

 

이 내용을 보면 우리가 의도했던 대로 가상 데이터가 생성되었음을 알 수 있습니다. 그러면 이번에는 이 지점별 데이터를 있는 그대로 글로벌 맵 상에 먼저 표출해봅시다. 다음과 같이 MAP 함수 및 MAPCONTINENTS 함수를 사용하여 글로벌 맵을 표출한 다음, SCATTERPLOT 함수를 사용하여 지점별 데이터를 원형 기호로 중첩하여 표출합니다.

 

win1 = WINDOW(DIMENSIONS=[800, 400], /NO_TOOLBAR)

m1 = MAP('geographic', LIMIT=[-90, -180, 90, 179.99], $

  MARGIN=0, font_size=10, /CURRENT)

m1.MapGrid.GRID_LONGITUDE = 30

m1.MapGrid.GRID_LATITUDE = 30

m1.MapGrid.LINESTYLE = 2

mc1 = MAPCONTINENTS()

ct = COLORTABLE(74, /REVERSE)

spl = SCATTERPLOT(lons, lats, RGB_TABLE=ct, MAGNITUDE=BYTSCL(data), $

  SYMBOL='circle', /SYM_FILLED, /OVERPLOT)

 

여기서는 세계 전체 지도가 표출되도록 하기 위하여 경도 및 위도의 범위를 그에 맞게 설정하였습니다. 그리고 각 지점별로 데이터 값의 크기에 따라 색상이 차등적으로 반영되도록 하기 위하여 SCATTERPLOT 함수에서 RGB_TABLE 및 MAGNITUDE 속성을 사용하였습니다. 이렇게 표출된 그림은 다음과 같습니다.

 

 

이제는 이와 같이 지점별로 불규칙하게 분포하는 데이터를 규칙 격자 데이터로 변환해봅시다. 여기서는 앞서 언급했듯이 GRIDDATA 함수를 사용하고 세부적으로는 Kriging 기법이 사용되도록 할 것입니다. 이 과정은 다음과 같습니다.

 

data_gridded = GRIDDATA(lons, lats, data, DIMENSION=[360, 180], /KRIGING)

HELP, data_gridded

PRINT, MIN(data_gridded), MAX(data_gridded)

 

여기서는 lons, lats, data 배열로 존재하는 지점별 데이터를 GRIDDATA 함수에 투입하여 Kriging 기법을 적용하여 360x180의 격자 구조를 갖는 2차원 데이터로 변환하여 그 결과를 data_gridded라는 배열로 얻은 것입니다. HELP 및 PRINT에 의하여 출력된 내용을 보면 그러한 사실을 확인할 수 있습니다.

 

DATA_GRIDDED    FLOAT     = Array[360, 180]

      13.8664      80.9980

 

그러면 이번에는 이 2차원 배열을 글로벌 맵 상에 중첩하여 표출해봅시다. 이를 위하여 또 다른 그래픽 창을 띄워서 앞서와 동일한 방식으로 글로벌 맵을 표출하되 이번에는 data_gridded 배열을 이미지의 형태로 중첩하여 표출하면 됩니다. 그 내용은 다음과 같습니다.

 

win2 = WINDOW(DIMENSIONS=[800, 400], /NO_TOOLBAR)

m2 = MAP('geographic', LIMIT=[-90, -180, 90, 179.99], $

  MARGIN=0, font_size=10, /CURRENT)

m2.MapGrid.GRID_LONGITUDE = 30

m2.MapGrid.GRID_LATITUDE = 30

m2.MapGrid.LINESTYLE = 2

i = IMAGE(data_gridded, RGB_TABLE=ct, $

  IMAGE_LOCATION=[-180, -90], IMAGE_DIMENSION=[360, 180], /OVERPLOT)

mc2 = MAPCONTINENTS()

 

여기서는 data_gridded 배열을 이미지의 형태로 표출하기 위하여 IMAGE 함수를 사용하였습니다. 또한 이미지가 지도 범위에 맞게 중첩되도록 하기 위하여 IMAGE 함수 내에서 IMAGE_LOCATION 및 IMAGE_DIMENSION 속성이 사용되었습니다. 이렇게 표출된 그림은 다음과 같습니다.

 

 

사실 이 그림을 보면 결과가 무난히 잘 나온 것으로 보이며, 아마도 더 이상 할 일은 없을 것 같기도 합니다. 그런데 이 그림을 다시 한번 잘 살펴봅시다. 특히 맵 상에서 우측 끝부분과 좌측 끝부분을 잘 보시기 바랍니다. 즉 경도상으로 +180인 부분과 -180인 부분입니다. 어차피 경도 +180과 -180은 구체(Sphere) 형태의 지구 표면 상에서는 서로 동일한 지점이기 때문에 이 두 부분은 구면 상에서 서로 이어집니다. 그런데 중첩된 격자 이미지를 보면 그 형태가 불연속적입니다. 이러한 문제점은 위의 그림을 구체 폴리곤 상에 입혀서 표출해보면 확실하게 확인이 가능합니다. 위의 그림을 이미지로 캡쳐하고 구체 폴리곤을 구현하여 그 위에 입히는 과정은 다음과 같습니다.

 

img = win2.CopyWindow(WIDTH=800)

MESH_OBJ, 4, vertices, polygons, REPLICATE(0.5, 101, 101)

oImage = OBJ_NEW('IDLgrImage', img)

vector = FINDGEN(101)/100.

texure_coordinates = FLTARR(2, 101, 101)

texure_coordinates[0, *, *] = vector # REPLICATE(1., 101)

texure_coordinates[1, *, *] = REPLICATE(1., 101) # vector

oSphere = OBJ_NEW('IDLgrPolygon', SHADING=1, DATA=vertices, $

  POLYGONS=polygons, COLOR=[255, 255, 255], $

  TEXTURE_COORD=texure_coordinates, $

  TEXTURE_MAP=oImage, /TEXTURE_INTERP)

XOBJVIEW, oSphere, XSIZE=600, YSIZE=600, SCALE=0.8

 

* 이와 같이 구체 형태의 폴리곤을 구현하고 그 위에 이미지를 입히는 방법에 관해서는 제가 예전에 올렸던 관련 게시물들이 있습니다. 따라서 관련된 더 자세한 내용이 궁금하시다면 이 게시물들(게시물 1, 게시물 2)을 참조하시기 바랍니다.

 

위와 같은 내용을 실행하면 이미지가 입혀진 구체 폴리곤이 구현되어 나타나며 마우스 클릭 & 드래그에 의하여 지구 구체를 이리저리 돌려볼 수 있습니다. 문제가 되는 부분을 보면 다음과 같습니다.

 

 

이 그림을 보면 경도가 +180 (또는 -180)인 부분에서 자오선을 따라 이미지의 불연속이 확연하게 드러납니다. 분명히 이것은 간과할 수 없는 문제입니다. 그리고 또 다른 문제도 있는데, 이 그림에서 북극 지역을 잘 보면 뭔가 길쭉한 형태의 패턴들이 북극점으로 수렴하는 것 같은 모습이 보입니다. 사실 이것 역시도 매우 부자연스러운 패턴입니다. 첫번째 문제는 결국 격자화 내삽의 과정에서 경도 +180과 -180이 서로 동일한 지점이라는 구체(Sphere) 표면의 특수성이 반영되지 않았기 때문에 발생한 것입니다. 그리고 두번째 문제도 맥락은 비슷합니다. 즉 경도선(수직 방향의)들의 간격이 극지방으로 갈수록 점점 좁아진다는 구체 표면의 특수성이 전혀 반영되지 않았기 때문에 발생한 문제입니다. 결국 우리가 얻은 내삽 결과는 사각형 평면에 가까운 공간에서라면 아무 문제가 없었겠지만, 지구 구체라는 입체적 형태의 표면이라는 기하학적 특성이 전혀 고려되지 않은 결과인 셈입니다.

 

사실 데이터의 지점들이 이와 같이 글로벌 맵 전체에 분포하는 경우가 아닌 국지적인 영역(예 : 한반도 영역) 내에서만 제한적으로 분포할 경우에는 거의 평면에 가까운 공간상에 존재하는 것이나 마찬가지이기 때문에 위와 같은 문제에 대하여 굳이 고민할 필요는 없습니다. 하지만 지금의 예와 같이 지구 구체 표면 상의 글로벌 맵 전체에 분포하는 지점들에 대한 데이터를 다룰 경우에는 분명히 큰 문제가 됩니다. 물론 다행히도 해결책은 존재합니다. 앞서 GRIDDATA 함수가 사용된 내용에서 /SPHERE 및 /DEGREES  두 속성을 추가적으로 사용하는 것입니다. 즉 다음과 같은 방식입니다.

 

data_gridded = GRIDDATA(lons, lats, data, DIMENSION=[360, 180], /KRIGING, /DEGREES, /SPHERE)

 

이와 같은 수정사항을 반영하여 다시 얻은 결과는 다음과 같습니다.

 

 

이 그림을 보면 경도 +180 및 -180인 부분에서 발생하던 이미지 패턴의 불연속 현상 그리고 북극 지역에서 보이던 이상한 패턴들은 완전히 사라지고 구체 표면이라는 특성이 잘 반영된 것을 확인할 수 있습니다. 구면상의 데이터라는 특성이 내삽 과정이 제대로 적용되어 산출된 결과이기 때문입니다. 여기서 해결책으로 사용된 /DEGREES 및 /SPHERES 속성은 GRIDDATA 함수에서 Kriging 외에 다른 세부 기법을 적용할 경우에도 얼마든지 사용이 가능합니다. 따라서 글로벌 맵 상에 분포한 지점 데이터를 규칙 격자 데이터로 변환해야 하는 경우에는 이러한 방법을 사용하는 것이 좋습니다. 그리고 글로벌 맵의 규모까지는 아니더라도 지점들의 분포 영역이 상당히 넓을 경우, 지점들의 분포 위치가 주로 고위도 지역일 경우, 그리고 지점들의 분포 위치가 경도 +180(즉 -180)인 부분과 가까울 경우에도 이와 같은 방법을 적용하는 것이 필요할 수도 있습니다. 제가 테스트 삼아서 한반도 및 주변 영역으로 지점들의 분포를 한정하여 얻어본 결과를 보여드리면 다음 그림과 같습니다.

 

 

여기서는 지점들이 경도 120~135, 위도 30~45의 범위 내에 존재하는 것으로 가정하였고, GRIDDATA 함수를 이용하여 격자화 내삽을 하는데 있어서 앞서 소개한 /DEGREES 및 /SPHERE 속성을 사용하지 않은 경우와 사용한 경우를 차례로 표출한 것입니다. 이와 같이 국지적인 영역으로 한정될 경우에는 지구 구면이라는 고려를 굳이 안해도 되긴 합니다. 실제로 위의 두번째 및 세번째 그림을 서로 비교해봐도 극명하게 차이가 나지는 않습니다. 하지만 잘 보면 미세하게 차이가 있긴 합니다. 아무래도 구면이라는 특성을 고려한 경우와 그렇지 않은 경우 사이에는 계산상으로는 약간의 차이는 존재할 수 밖에 없습니다. 어떤 방향으로 갈 것이냐는 결국 데이터의 특성이나 제반 상황들을 고려해서 연구자가 최종적으로 판단해야 할 문제라고 봐야 할 것입니다.