GRIDDATA 함수를 사용한 2차원 데이터의 근사 [1]

오늘 소개해드릴 GRIDDATA는 2차원 공간상에 분포하는 데이터값들을 내삽하여 우리가 원하는 형태의 새로운 데이터로 만들어주는 역할을 하는 내장함수입니다. 사실 IDL에서 꽤 오래전부터 지원되어오던 함수임에도 불구하고, 그 사용법에 대한 설명(특히 우리말로 된)을 찾아보기가 힘든 편입니다. 그래서 기본적인 내용이나마 한번 정리를 해보고자 합니다.

IDL 사용자들 중 이 함수의 존재나 사용법에 대하여 한번이라도 관심을 가져본 유저들은 대개 동기가 비슷합니다. 주로 2차원 공간상에서 불규칙한 분포를 하는 데이터 포인트들을 규칙 격자화하는 작업을 해야 하는 경우입니다. 그리고 여기서 말하는 2차원 공간은 주로 지구 표면상의 지리적 공간인 경우가 많습니다. 사실 이러한 목적이라면 GRIDDATA 말고도 쓸만한 내장함수가 또 있기는 합니다. 바로 SPH_SCAT이라는 함수인데요. 이 함수의 사용법에 대해서는 제가 예전에 한번 다뤘던 적이 있습니다(링크 참조).

그러면 SPH_SCAT과 GRIDDATA는 어떤 차이가 있는가에 대한 궁금증이 생깁니다. 일단 단순하게 얘기한다면, SPH_SCAT은 GRIDDATA의 특수한 하나의 경우에 해당된다고 볼 수 있습니다. 즉, GRIDDATA가 좀 더 일반적인 기능을 수행합니다. 따라서 지구 구체상에서 경위도 기반의 분포를 하는 데이터 포인트들을 다루는 경우에 한해서는 SPH_SCAT을 사용해도 됩니다. 그러나 구체상에서 뿐 아니라 그냥 평면상에 존재하는 데이터 포인트들에 대하여 좀 더 다양한 내삽 알고리즘을 사용하여 자유로운 격자 변환이 필요한 경우에는 GRIDDATA를 사용하는 것이 더 좋습니다. 어쨌든 오늘은 GRIDDATA에 촛점을 맞춰서 소개를 해보겠습니다.

IDL 도움말에서 GRIDDATA에 관한 내용을 찾아보면 양도 많고 내용도 좀 복잡합니다. 아무래도 다양한 알고리즘 및 사용방법들이 존재하기 때문에 어쩔 수가 없는데요. 가장 기본적인 사용 문법은 다음과 같습니다.

result = GRIDDATA(X, Y, F)

여기서 X, Y, F는 각각 데이터 포인트들의 X좌표, Y좌표, 값에 해당되는 배열입니다. 만약 N개의 데이터 포인트들이 있는 상태라면 X, Y, F 각각은 N개의 값들로 이루어진 1차원 배열이 되어야 합니다. 이렇게 1차원 배열이 되어야 한다는 것이 불규칙한 분포를 하는 데이터의 경우는 자연스럽습니다. 그런데 오늘 제가 예제로 사용할 데이터는 규칙 격자 분포를 하는 데이터입니다. 단, 유의미한 데이터값이 존재하지 않는 격자들이 군데군데 존재하는, 즉 여기저기 구멍들이 좀 뚫린 그런 데이터를 예제로 사용하고자 합니다. 그래서 GRIDDATA를 사용함으로써 이 구멍난 격자점들에 대하여 그럴듯한 값들을 채워넣어보고자 합니다.

예제 데이터는 다음과 같이 50X50의 격자구조를 갖는 배열을 SIN, COS, RANDOMU 등의 함수를 사용하여 다음과 같이 임의로 만들어 보았습니다. 대략 0~10 정도의 범위의 값을 가지며, 중간중간 NaN값을 일부러 집어넣어서 구멍난 데이터가 되도록 하였습니다. 그리고 X, Y축 방향에 해당되는 격자의 경도 및 위도 값들로 구성된 배열을 각각 lons, lats라는 이름으로 만들었습니다. 그야말로 격자점의 경도 및 위도 방향에 해당되는 각 50개씩의 값들로 이루어진 1차원 배열입니다. 따라서 lons, lats는 50X50=250개가 아닌 50개의 값들로 이루어진 배열이라는 점을 일단 염두에 두어야 합니다. 여기서는 경도가 동경 120~130도, 위도는 북위 30~40도의 범위에 해당되는 것으로 가정해보았습니다. 어쨌든 이렇게 만들어진 가상 데이터를 이미지의 형태로 표출해보면 그 모습은 대략 다음 그림과 같습니다.

xgz = 50

ygz = 50

lons = FINDGEN(xgz)/5 + 120

lats = FINDGEN(ygz)/5 + 30

data = SIN(lons/3) # COS(lats/3)

data = data*5+5+RANDOMU(seed, xgz, ygz)

HELP, data

PRINT, MIN(data), MAX(data)

PRINT, MIN(lons), MAX(lons)

PRINT, MIN(lats), MAX(lats)

data[15:18, 12:13] = !values.f_nan

data[35:38, 29:33] = !values.f_nan

data[20:22, 41:44] = !values.f_nan

win0 = WINDOW(DIMENSIONS=[500, 500])

im0 = IMAGE(data, RGB_TABLE=59, MARGIN=0.1, /CURRENT)

이 그림은 NG 체계의 IMAGE 함수를 사용하여 표출한 것입니다. 컬러테이블은 59번을 사용하였습니다. 그리고 구멍난 부분에 해당되는 격자점은 NaN값을 갖지만, 이미지로 표출할 경우 이러한 부분은 아예 표출이 되지 않습니다. 따라서 그림과 같이 NaN값인 격자는 아예 비어있는 형태로 표출이 된 것입니다. 그러면 이제는 GRIDDATA를 사용하여 저 비어있는 격자점들에 그럴듯한 값들을 채워넣어봅시다. 이를 위해서는 먼저 GRIDDATA 함수에 투입할 X, Y, F를 준비해야 합니다. 그런데 여기서 일단 생각을 해야 할 것은,X, Y, F가 모두 모든 격자점들을 커버하는 1차원 배열이 되어야 한다는 것입니다. 현재 우리가 갖고 있는 lons, lats, data는 이러한 조건에 아직은 부합되지 않습니다. 현재의 상태는 다음과 같습니다.

lons : X축 방향 50개의 경도값들로 이루어진 1차원 배열

lats : Y축 방향 50개의 위도값들로 이루어진 1차원 배열

data : 50X50의 구조를 갖는 2차원 배열

따라서 이 세 개의 항목들을 GRIDDATA 함수에 투입 가능한 형태로 바꿔줘야 합니다. 일단 이 과정은 다음과 같습니다.

data_tmp = REFORM(data)

xys = ARRAY_INDICES(data, INDGEN(xgz, ygz))

xxs = REFORM(xys[0, *])

yys = REFORM(xys[1, *])

여기서는 가장 먼저 data를 1차원 배열로 바꿔주기 위하여 REFORM이란 함수가 사용되었습니다. REFORM은 어떤 배열을 총 구성원소 갯수가 변하지 않는 한도내에서 내가 원하는 구조로 변형하는데 사용되는 함수입니다. 대상 배열외 별도의 인수가 붙지않는 지금과 같은 경우에는 50X50의 2차원이던 구조를 250의 1차원으로 변형해줍니다. 따라서 data_tmp는 250개의 값들로 이루어진 1차원 배열이 됩니다. 그리고 나머지 세 줄은 설명하자면 약간 복잡합니다. 앞서 얻은 data_tmp안에 있는 250개의 값 각각에 대한 X좌표(경도)값의 인덱스(index)들로만 새롭게 구성된 배열이 xxs, Y좌표(위도)값의 인덱스(index)들로만 새롭게 구성된 배열이 yys라는 정도로 이해하시기 바랍니다. 여기서는 사실 ARRAY_INDICES라는 함수도 사용이 되었는데, 이 과정에 대한 세세한 설명은 여기서는 일단 생략하기로 하겠습니다. 하여간 이와 같은 과정을 거쳐서 얻어진 세 개의 항목들을 다음과 같이 정리해보았습니다.

xxs : 250개 데이터 포인트들의 X좌표(경도)값 인덱스들로만 구성된 1차원 배열

yys : 250개 데이터 포인트들의 Y좌표(위도)값 인덱스들로만 구성된 1차원 배열

data_tmp : 50X50의 구조를 250의 1차원 구조로 변형하여, 250개의 데이터값들로 구성된 1차원 배열

그러면 이 세 항목들을 바로 GRIDDATA 함수에 투입하면 될 것 같은데요. 안타깝게도 아직 하나의 과정이 더 남아있습니다. 그것은 NaN값들을 제외시켜주는 과정입니다. 왜냐하면 사실 GRIDDATA는 자체적으로 NaN값을 걸러주는 기능이 포함되어 있지 않기 때문입니다. 마치 /NAN과 같은 키워드가 있을 것도 같지만 실제로는 그렇지가 않습니다. 따라서 좀 수고스럽지만, NaN값을 제외한 나머지 유효한 값들로만 이루어진 배열들을 얻어야만 합니다. 이 과정은 다음과 같이 FINITE, WHERE 등의 함수를 사용하여 해결하였습니다.

ww = WHERE(FINITE(data, /NAN) EQ 0, count)

xxs = xxs[ww]

yys = yys[ww]

data_tmp = data_tmp[ww]

참고로 FINITE라는 함수는 어떤 값이 NaN이냐를 판단하는 역할을 합니다. 위와 같이 /NAN이란 키워드와 함께 사용하면, 인자값이NaN일 경우에는 1을, 그렇지 않을 경우에는 0을 돌려받게 됩니다. 잘 알아두면 유용하게 활용이 가능한 함수라는 점 참조해두면 좋습니다. 어쨌든 이제는 정말로 GRIDDATA 함수에 투입할 준비가 다 끝났습니다. 실제 GRIDDATA 함수의 사용 과정은 몇가지 경우들로 나누어 소개해보고자 합니다. 먼저 X, Y축 방향 좌표인 경도와 위도는 일단 무시하고 그냥 인덱스 기반으로만 시도해봅시다. 물론 경위도를 무시했기 때문에 그냥 평면 데이터인 상태로 처리하는 셈입니다. 이 경우에는 GRIDDATA 함수를 다음과 같은 방식으로 사용하면 됩니다.

data_f = GRIDDATA(xxs, yys, data_tmp, DIMENSION=[xgz, ygz])

여기서는 xxs, yys, data_tmp가 인자로 투입이 되었고, DIMENSION이라는 키워드가 추가적으로 사용이 되었습니다. 이 키워드에 부여된 값은 결국 [50, 50]이 되는데, 결과인 data_f를 50X50인 배열의 형태로 얻겠다는 의미입니다. 물론 이 결과를 보면, 원본 데이터인 data에서 NaN이었던 구멍들에 어떤 값들이 채워져 있게 됩니다. 비교를 위하여 원본인 data, 처리결과인 data_f를 다음과 같은 과정으로 나란히 표출하여 보았습니다. 그 결과는 다음 그림과 같습니다.

win1 = WINDOW(DIMENSIONS=[500*2, 500])

im1o = IMAGE(data, RGB_TABLE=59, MARGIN=0.1, TITLE='Original', $

  /CURRENT, LAYOUT=[2, 1, 1])

im1f = IMAGE(data_f, RGB_TABLE=59, MARGIN=0.1, TITLE='Filtered', $

  /CURRENT, LAYOUT=[2, 1, 2])

그림에서 왼쪽이 원본인 data이고 오른쪽은 처리결과인 data_f입니다. 보이는 것처럼 원래 비어있던 부분에 값들이 채워져 있는 것을 눈으로 확인할 수 있습니다. 원래 비어있던 부분이 세 개였는데, 왼쪽의 두 군데는 그럴듯한 값들로 채워진 것 같지만 가장 오른쪽 것은 뭔가 주변과 괴리가 좀 느껴집니다. 그 원인은 GRIDDATA에서 내부적으로 사용된 알고리즘의 문제입니다. IDL 도움말에서 이함수에 관한 내용을 찾아보면, 지원되는 알고리즘의 종류가 총 10개로 나옵니다. 일단 여기서는 기본 알고리즘인 Inverse Distance가 사용된 상태인데, 이 알고리즘이 지금의 경우에 대해서는 그리 만족스럽지는 않은 결과를 준 것으로 보입니다.

따라서 이런 경우에는 다른 알고리즘을 시도해봐야 합니다. 아직 내용이 더 남아있고 이번 게시물의 내용이 많이 길어진 관계로, 이 얘기는 이어지는 후속편 [2]에서 다루도록 하겠습니다.