Kriging 기법의 이용 방법

Kriging 기법은 불규칙한 분포를 하는 데이터들을 규칙적인 격자 데이터로 변환하는데 있어서 자주 사용되는 알고리즘입니다. 불규칙 분포 데이터를 규칙 격자 데이터로 변환하는 방법에 관해서는 제가 예전에 이 블로그에서 "불규칙 분포 데이터를 규칙 격자화된 데이터로 만들기"라는 제목으로 3회에 걸쳐 올렸던 게시물들(1, 2, 3)에서도 소개했던 바 있습니다. 다만 이 게시물들에서는 GRIDDATA 함수를 위주로 하여 이 함수 내에서 지원되는 다양한 내삽 기법들을 조금씩 소개하였고, 그 중 하나로서 Kriging 기법의 적용 방법에 관하여 매우 간단하게 소개를 했었습니다. 오늘은 IDL에서 Kriging 기법을 사용할 수 있도록 해주는 전용 내장함수인 KRIG2D 함수의 사용법을 위주로 하여 조금만 더 심층적으로 소개를 해보고자 합니다. 먼저 예제 데이터가 필요한데, 불규칙한 분포를 하는 가상의 데이터를 다음과 같이 생성해 봅시다.

 

n = 30

seed = -4

lons = 126.5+RANDOMU(seed, n)*2.7

lats = 35+RANDOMU(seed, n)*2.8

data = 20+RANDOMU(seed, n)*10

HELP, lons, lats, data

PRINT, MIN(data), MAX(data)

 

비록 가상의 데이터이긴 하지만, 나름대로는 한반도 내에서 불규칙하게 분포하는 30개의 지점들을 가정하고 각 지점의 경위도 좌표 및 관측값을 비교적 그럴듯하게 생성하기 위하여 RANDOMU 함수를 사용하여 난수로 발생시킨 것입니다. 지점들의 경도는 동경 126.5도에서 128.2도 사이가 되도록 하였고 위도는 35.0도에서 37.8도 사이가 되도록 하였습니다. 그리고 관측값은 20~30의 범위를 갖도록 하였습니다. 이렇게 가상으로 생성된 데이터를 원래의 데이터라고 가정하겠습니다. HELP 명령에 의하여 출력된 내용을 먼저 보면 다음과 같습니다.

 

LONS            FLOAT     = Array[30]

LATS            FLOAT     = Array[30]

DATA            FLOAT     = Array[30]

 

이와 같이 lons, lats, data는 30개의 지점들에 대한 경도, 위도, 관측값에 대응되는 1차원 배열들입니다. 어차피 불규칙한 분포를 하는 지점별 데이터는 이러한 형태로 얻어질 가능성이 큽니다. 그러면 우선 이 데이터의 분포 양상을 시각적으로 표출해보는 것이 필요하기 때문에 다음과 같이 한반도 지도를 먼저 그리고 그 위에 30개의 지점별 데이터들이 표시되도록 합시다.

 

win = WINDOW(DIMENSIONS=[1200, 600], /NO_TOOLBAR)

m1 = MAP('geographic', LIMIT=[34, 125, 39, 130], $

  MARGIN=[0.12, 0.2, 0.05, 0.1], /CURRENT, LAYOUT=[2, 1, 1])

mc1 = MAPCONTINENTS(/HIRES)

ct = COLORTABLE(74, /REVERSE)

sp = SCATTERPLOT(lons, lats, MAGNITUDE=BYTSCL(data, MIN=20, MAX=30), $

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

m1.MapGrid.LABEL_POSITION = 0

m1.MapGrid.LINESTYLE = 1

m1.MapGrid.HORIZON_LINESTYLE = 0

cb1 = COLORBAR(TARGET=sp, RANGE=[20, 30], /BORDER, $

  POSITION=[0.14, 0.95, 0.4, 0.98])

 

여기서는 MAP 함수를 사용하여 바탕지도를 먼저 표출하면서 경위도 범위를 한반도 영역으로 제한하였습니다. 그리고 각 지점별 관측값들이 각각의 위치에 표시되도록 하기 위하여 SCATTERPLOT 함수를 사용하였습니다. 각 지점별 관측값이 20~30의 범위에 대하여 크기에 따라 74번 컬러테이블의 색상들이 반영되도록 하기 위하여 MAGNITUDE 키워드 및 BYTSCL 함수를 사용하였습니다. 그리고 그림의 상단에 컬러바도 함께 표시되도록 하였습니다. WINDOW 함수의 내용을 보면 그래픽창이 가로 방향으로 긴 형태로 정의되었고 LAYOUT 키워드도 사용되었습니다. 이렇게 하면 SCATTERPLOT 함수로 표출되는 그림은 그래픽창의 왼쪽에 배치됩니다. 오른쪽 부분은 일단 빈 공간으로 놔두었는데, 나중에 Kriging 기법을 적용한 결과 그림을 이 부분에 배치하려는 의도입니다. 일단 여기까지의 내용을 실행하면 다음과 같은 그림을 얻게 됩니다.

 

그러면 이제 Kriging 기법을 적용하여 규칙적인 격자 분포를 하는 데이터로 변환하는 작업을 해봅시다. 원래 IDL에서는 Kriging 기법을 사용할 수 있도록 해주는 내장함수가 두가지가 있습니다. 바로 KRIG2D 함수와 GRIDDATA 함수입니다. 하지만 앞서 언급했듯이 오늘은 KRIG2D 함수를 이용하는 방법 위주로 소개하겠습니다. 그러면 KRIG2D 함수를 다음과 같이 사용해 봅시다.

 

data_new = KRIG2D(data, lons, lats, $

  EXPONENTIAL=[1.0, 0], $

  GS=[0.01, 0.01], BOUND=[126, 35, 130, 38])

HELP, data_new

PRINT, MIN(data_new), MAX(data_new)

 

여기서 처음 세개의 인수들은 원래 데이터인 data, lons, lats를 그대로 순서대로 사용하면 됩니다. 그리고 우리의 목표는 규칙격자화된 결과 데이터를 얻는 것이기 때문에, 규칙 격자들의 분포 형태를 우리가 직접 정의해야 합니다. 즉 규칙격자들이 점유하게 될 범위와 격자 해상도를 경도 및 위도 방향으로 정의하는 것입니다.이것은 GS 키워드와 BOUND 키워드를 사용하여 설정하면 됩니다. 위의 내용에서는 경도 방향으로 126도부터 130도까지 0.01도 간격으로 그리고 위도 방향으로 35도부터 38도까지 0.01도 간격으로 격자들이 분포하도록 정의한 것입니다. 즉 GS 키워드는 격자 해상도, BOUND 키워드는 격자들의 범위에 해당됩니다. 그리고 EXPONENTIAL이라는 키워드가 있고 이 키워드에 대하여 [1.0, 0]이라고 설정되어 있습니다. 이 부분에 관하여 설명하자면 Kriging 기법의 기본적인 이론에 대한 약간의 이해가 필요한데, 일단 IDL 도움말에서 KRIG2D 함수에 관한 내용을 보면 다음과 같이 소개되어 있습니다.

 

이것은 Kriging 기법을 적용하여 각 격자점에 대한 계산값을 산출하기 위하여 주변의 다른 지점들의 값들에 반영하는데 있어서 각 지점별로 어느 정도의 Weight를 줄 것이냐를 결정하는 Variogram 모형이라고 부릅니다. 기본적인 가정은 거리가 가까운 지점들이 Weight가 크고 거리가 먼 지점들은 Weight가 작다는 것입니다. 다만 거리가 멀어짐에 따라 Weight가 작아지는 패턴을 어떻게 정의하느냐의 문제인데, KRIG2D 함수에서는 총 4종의 모형들(Linear, Exponential, Gaussian, Spherical)이 지원되고 있습니다. 그래서 이들 중 어떤 모형을 적용할 것인지를 결정하여 해당 키워드 및 관련 인자들을 사용해야 합니다. 위의 예제 코드에서는 Exponential 모형을 적용하기 위하여 EXPONENTIAL 키워드를 사용하였고, 이 키워드에 대하여 [1.0, 0]이라는 배열이 부여되어 있습니다. 이러한 배열은 기본적으로 [R, N]의 형태로 구성되는데, R은 영향반경(Radius)에 해당되고 N은 Nugget이라는 개념에 해당됩니다. 그래서 [R, N]과 같은 배열을 키워드에 부여하는 방식입니다. 경우에 따라서는 Scale에 해당되는 값을 추가하여 [R, N, S]의 형태로 사용하기도 합니다. Kriging 기법에서 R, N, S의 의미에 관한 더 자세한 내용이 궁금하시다면, 검색을 하여 찾아보시길 권장합니다. 제가 여기서 더 깊이 설명드리기는 어려울 것 같습니다. 그냥 좀 쉽게 생각해본다면, R의 값이 클수록 더 많은 주변 지점들을 계산에 반영하게 되어 결과의 전반적인 패턴이 다소 부드러워지는 효과가 있는 반면, R의 값이 작을수록 가까운 주변 지점들만을 계산에 반영하게 되어 결과의 전반적인 패턴이 다소 거칠어지는 효과가 있다는 정도로 이해하면 어떨까 합니다(제 생각입니다). 어쨌든 위와 같은 과정에 의하여 격자화된 결과 데이터를 얻을 수 있습니다. 일단 HELP 및 PRINT 명령에 의하여 출력된 배열의 구조 및 배열 내 최소/최대값은 다음과 같습니다.

 

DATA_NEW        FLOAT     = Array[401, 301]

      20.3849      29.4814

 

즉 401x301의 규칙적인 격자 구조를 갖도록 변환된 결과 데이터가 data_new의 이름으로 얻어진 것입니다. 물론 401x301의 구조를 갖는 것은 격자들의 범위 및 해상도에 따른 것이고, 만약 격자 범위나 해상도가 다를 경우에는 결과 데이터 배열의 구조도 달라질 것입니다. 그러면 이 데이터를 표출하여 원래의 지점별 데이터와 비교를 해봅시다. 앞서 정의했던 그래픽창의 오른쪽 공간에 결과가 표출되도록 하기 위하여 다음과 같은 표출 과정을 추가합시다.

 

m2 = MAP('geographic', LIMIT=[34, 125, 39, 130], $

  MARGIN=[0.12, 0.2, 0.05, 0.1], /CURRENT, LAYOUT=[2, 1, 2])

ct = COLORTABLE(74, /REVERSE)

im = IMAGE(data_new, RGB_TABLE=ct, MIN_VALUE=20, MAX_VALUE=30, $

  IMAGE_LOCATION=[126, 35], IMAGE_DIMENSION=[4, 3], GRID_UNITS=2, /OVERPLOT)

mc2 = MAPCONTINENTS(/HIRES)

m2.MapGrid.LABEL_POSITION = 0

m2.MapGrid.LINESTYLE = 1

m2.MapGrid.HORIZON_LINESTYLE = 0

cb2 = COLORBAR(TARGET=im, /BORDER, $

  POSITION=[0.64, 0.95, 0.9, 0.98])

 

이 과정은 기본적으로는 왼쪽 그림을 표출할 때와 유사합니다. 다만 바탕지도를 표출한 상태에서 그 위에 2차원 격자 데이터를 중첩하기 위하여 IMAGE 함수를 사용하였고, 이 때 IMAGE_LOCATION 및 IMAGE_DIMENSION 속성을 함께 사용하여 2차원 이미지의 위치 및 범위를 적절하게 설정한 것입니다. 즉 지도상에서 이미지의 위치는 이미지의 좌측하단 꼭지점의 위치(경위도 좌표)가 되어야 하므로 [126, 35]가 되고, 이미지의 경도 및 위도 방향의 크기는 각각 4도 및 3도가 되도록 하였습니다. 이렇게 하면 결과는 다음 그림과 같습니다.

 

비교를 위하여 영향반경의 값을 1.0 대신 2.0으로 대체해 봅시다. 위의 예제 코드에서 다음과 같이 해당 부분에 대해서만 값을 바꾸고 다시 실행하면 됩니다.

 

data_new = KRIG2D(data, lons, lats, $

  EXPONENTIAL=[2.0, 0], $

  GS=[0.01, 0.01], BOUND=[126, 35, 130, 38])

 

이렇게 변경하여 다시 실행하면 결과는 다음 그림과 같습니다.

 

아무래도 영향반경이 클수록 각 격자점별로 주변의 더 많은 지점들이 계산에 반영되기 때문에 전반적인 패턴에 있어서 약간의 변화가 있음을 확인할 수 있습니다(물론 변화가 아주 크지는 않습니다). 그러면 이번에는 Variogram 모형을 Exponential에서 Spherical로 변경해 봅시다. 영향반경은 일단 1.0으로 설정합니다.

 

data_new = KRIG2D(data, lons, lats, $

  SPHERICAL=[1.0, 0], $

  GS=[0.01, 0.01], BOUND=[126, 35, 130, 38])

 

이렇게 하면 결과는 다음 그림과 같습니다.

 

아무래도 모형이 변경된만큼 결과의 모습도 Exponential의 경우와는 나름대로 다르게 나타납니다. 그리고 이 상태에서 다음과 같이 영향반경만 2.0으로 변경해 봅시다.

 

data_new = KRIG2D(data, lons, lats, $

  SPHERICAL=[2.0, 0], $

  GS=[0.01, 0.01], BOUND=[126, 35, 130, 38])

 

이렇게 하면 결과는 다음 그림과 같습니다.

이 모습을 보면 Spherical 모형의 경우에는 영향반경에 따른 차이가 앞서 Exponential에 비하면 약간 더 두드러지게 나타나는 것 같습니다. 그 외에 Gaussian 모형도 적용해볼 수 있는데, 앞의 두 모형과는 차이가 꽤 납니다. 다음과 같이 해봅시다.

 

data_new = KRIG2D(data, lons, lats, $

  GAUSSIAN=[1.0, 0], $

  GS=[0.01, 0.01], BOUND=[126, 35, 130, 38])

 

다만 여기서는 오른쪽 그림에서 이미지 표출을 담당하는 IMAGE 함수의 내용 일부를 다음과 같이 변경하는 것이 좋습니다.

 

im = IMAGE(data_new, RGB_TABLE=ct, MIN_VALUE=10, MAX_VALUE=60, $

  IMAGE_LOCATION=[126, 35], IMAGE_DIMENSION=[4, 3], GRID_UNITS=2, /OVERPLOT)

 

그 이유는 결과 데이터인 data_new의 값 범위가 원래 데이터의 20~30과는 확연히 다르게 10~60 정도의 범위로 나타나기 때문입니다. 실제로 위의 변경사항들만 반영하여 다시 실행해보면 결과는 다음 그림과 같습니다.

 

이와 같이 Gaussian 모형에 의한 결과는 앞서 Exponential, Spherical 모형의 경우와는 꽤 다르게 나타납니다. 이와 같이 Kriging 기법의 적용에 있어서 Variogram의 모형을 무엇으로 정의하느냐 그리고 관련 인자들의 값을 어떻게 설정하느냐에 따라서 결과는 서로 다르게 나타날 수 밖에 없습니다. 물론 정답이라는 것은 사실상 존재하지 않는다고 봐야 합니다. 이러한 작업에 있어서는 다양한 기법들과 인자값들을 시도하면서 최적의 결과로 향해가는 과정에서 여러번의 시행착오를 겪게 될 가능성이 높습니다. 그리고 이러한 우여곡절을 거쳐서 어떤 결과가 나오든간에 그 결과가 과연 의미있고 적절한 것인가에 대한 판단은 프로그래머 및 그 분야의 연구자들의 몫입니다. 그리고 결과로 얻어진 격자화 데이터의 값들, 특히 원래 지점이 존재하지 않던 부분의 격자에 대한 계산값이 어느 정도로 신뢰할만하냐에 대해서도 매우 신중한 판단이 필요하다는 점도 반드시 유념해야 합니다.