불규칙 분포 2차원 데이터 격자화

Q. x좌표, y좌표에서의 실험 측정값을 z라고 하여 x, y, z 배열이 있습니다. 그런데 x, y 좌표의 간격이 일정하지 않고 불규칙하게 측정이 되어 있습니다. 이를 일정한 간격으로 재배치 하고자 합니다.
A. 실험실에서 측정하는 결과라면 일정 간격으로 측정을 하여 그 값 자체를 2차원 배열에 담을 수 있겠지만, 지형이나 기후, 수온 등의 측정은 일정한 간격으로 측정한다는 것이 거의 불가능합니다. IDL에서는 Triangulation과 TriGrid라는 빠르고 안정적인 명령문을 통해, 불규칙하게 측정된 데이터라도 일정 간격의 데이터로 변환할 수 있습니다.
2차원 배열의 형태로 저장되는 데이터는, 이미지, 3차원 표면 그래프, 등고선 등으로 디스플레이 될 수 있습니다. 2차원의 측정 데이터(예를 들면, 어떤 경도, 위도에서 측정한 어떤 값)라도 배열처럼 일정 간격으로 측정된 데이터가 아니면, 보간법을 이용해서, 일단은 2차원 배열 형태의 데이터로 만들어야 합니다. (여기서 배열의 간격은 0, 1, 2, 3, …. 으로 일정하다는 점을 기억해 주시기 바랍니다. 그러므로 일정 간격의 데이터로 다시 만든다는 의미입니다).

데이터를 모의로 만들어 보겠습니다

IDL> seed=1L
IDL> x=RandomU(seed, 32)
IDL> y=RandomU(seed, 32)
IDL> z=exp(-3*((x-0.5)^2+(y-0.5)^2))

여기서는 서른 두 지점에서의 측정 결과를 마구잡이수를 이용해 모의로 만들었지만, 그 의미는 서른 두개의 임의의 위치 x, y에서 측정된 z값을 의미합니다(실험, 관측에서 흔히 등장하는 상황입니다).

측정 위치 x와 y의 분포를 확인해 봅시다

IDL> plot, x, y, psym=1

불규칙 분포 좌표

불규칙 분포 좌표

원리는 이렇습니다

평면상에 여러 개의 점이 있을 때 이를 세 개씩 조합하여 여러개의 삼각형을 만들 수 있습니다. 들로네의 삼각형(Delaunay triangulation)이란 각각의 세 점을 꼭지점으로 하는 삼각형들이 서로 간의 영역을 침범하지 않도록 조합하는 것입니다(IDL에는 triangulate라는 해당 프로시저가 있습니다). 이 후에 이 평면 상에 일정한 간격의 격자를 치면, 각각의 격자 위치가 포함되는 삼각형이 있을 것입니다. 이 삼각형 세 점의 데이터 값(z값)으로 부터, 격자 위치의 데이터 값을 추정해 내면 됩니다(IDL에는 trigrid라는 해당 함수가 있습니다).

삼각형 조합(Triangulation)

x와 y가 어떤 평면 상의 점들의 위치 정보를 담고 있는 배열이라면, triangulate 프로시저를 이용하여 들로네의 삼각형 조합을 만들 수 있습니다.

IDL> triangulate, x, y, tr

이 때, tr 배열은 3*?의 2차원 배열이 되는데, 3은 각각 삼각형의 세 꼭지점을 가리키는 x와 y의 index(배열 번호)를 의미하고, ?는 x와 y의 개수 및 분포에 따라 달라집니다. 각 위치들을 어떻게 삼각형으로 조합했는지 확인해 보려면 다음과 같이 해 볼 수 있습니다.

IDL> plot, x, y, psym=1, TITLE=’Triangulation’
IDL> FOR i=0, N_ELEMENTS(tr)/3-1 do begin & $
t=[tr[*,i], tr[0,i]] & plots, x[t], y[t] & ENDFOR

Triangulate 결과

Triangulate 가 하는 일

일정 간격 위치에서의 값을 계산(Trigrid)

새로운 데이터 포인트에서의 값은 어떻게 계산할까요? 그 데이터 포인트를 포함하고 있는 삼각형의 세 꼭지점에서 측정된 값으로부터 적당히 내삽(interpolation)하여 계산할 것입니다(대충 계산한다는 뜻은 아닙니다. 거리가 가까운 쪽에 비중을 높게 두어 계산합니다). TRIANGULATE는 이런 계산을 위해 미리 삼각형 구역을 할당해 놓은 것입니다. 어쨌거나 이런 계산은 TRIGRID가 알아서 합니다.

IDL> Surface, TRIGRID(x, y, z, tr)

[x, y, z, tr], 우리가 알고 있는 모든 값들이 들어갔습니다. 이렇게 그릴 수 있는 것은 TRIGRID가 2차원 배열을 리턴하는 함수이기 때문입니다.

IDL> help, TRIGRID(x, y, z, tr)
FLOAT     = Array[51, 51]

TRIGRID로 생성된 2차원 배열은 등고선, 3차원 그래프, 이미지 등 어떤 방법으로도 디스플레이 할 수 있음은 물론입니다. 꽤 복잡한 일을 한 것 같지만, 실제로는 두 줄입니다.

요약하면 두 줄입니다.

일정하지 않은 간격으로 x, y 위치에서 측정된 값을 z라고 하면(x, y, z는 같은 길이의 1차원 배열), 이 데이터를 2차원 배열로 재구성하기 위한 과정은 다음과 같습니다.

IDL> TRIANGULATE, x, y, tr
IDL> array=TRIGRID(x, y, z, tr)

TRIGRID의 옵션

1) 출력하는 배열의 간격 및 범위 지정
IDL 온라인 도움말에 나와 있는 TRIGRID의 사용 문법은 다음과 같습니다.

Result = TRIGRID( X, Y, Z, Triangles [, GS, Limits] )

이 때 옵션인 GS는 두 개의 요소를 가지는 배열로서, 각각 x간격, y간격을 의미합니다. Limits는 네 개의 요소를 가지는 배열로서 [x최소, y최소, x최대, y최대]입니다.

다음과 같이 사용할 수 있습니다.

New_grid=trigrid(x, y, z, tr, [0.01, 0.02], [0.0, 0.0, 1.0, 1.0])

x, y 각각 0~1의 범위에서 x간격은 0.01, y간격은 0.02인 2차원 배열을 리턴하라는 의미입니다.

2) 출력하는 배열 크기 지정

TRIGRID함수는 디폴트로 51*51의 2차원 행렬을 리턴합니다. NX, NY 키워드를 사용하여 리턴하는 배열의 크기를 임의 설정할 수 있습니다.

IDL> Surface, TRIGRID(x, y, z, tr, NX=100, NY=100)

3) Smooth Fitting

배열의 크기를 늘려도, 데이터는 매끄럽게 보이지 않는다는 것을 알 수 있습니다. 이는 디폴트로 선형 보간(Linear Interpolation)에 의한 계산을 하기 때문인데, /QUINTIC 키워드를 사용하여 데이터가 매끄럽게(smooth) 되도록 할 수 있습니다. 다음의 두 결과를 비교해 보세요.

IDL> contour, TRIGRID(x, y, z, tr), nlevels=10
IDL> contour, TRIGRID(x, y, z, tr, /QUINTIC), nlevels=10

Linear Interpolation vs Quintic

Linear Interpolation vs Quintic

4) 외삽(Extrapolate)

TRIGRID는 기본적으로 내삽만 하기 때문에, 아예 측정한 자료가 없는 곳(그림에서 좌측 아래쪽 대부분과 우측 아래쪽 일부)은 계산을 하지 못합니다. 이 영역에는 어떠한 삼각형도 없기 때문입니다. 그렇지만, 외삽(Extrapolate)을 통해 이 영역의 값도 추정할 수 있습니다(내삽에 비해 신뢰도는 떨어집니다).
외삽을 하기 위해서는 TRIANGULATE 단계에서 변수를 하나 더 지정해 주어야 합니다.

IDL> triangulate, x, y, tr, b                ;변수 이름은 꼭 b가 아니어도 됩니다.

이렇게 하면, 변수 b에는 모든 가장 바깥쪽에 있는 x, y 값의 배열 번호(index)들이 들어갑니다. 순서는 반시계 방향인데, 이 번호의 x,y 좌표를 연결하면, 모든 x, y 좌표를 포함하는 폐곡선이 됩니다(boundary). 그리고 나서,

IDL> contour, TRIGRID(x, y, z, tr, EXTRAPOLATE=b)

와 같이 EXTRAPOLATE 키워드 값으로 b배열을 넘겨 주면, 데이터가 없는 부분도 외삽을 하여 추정합니다. EXTRAPOLATE 키워드가 설정되면, /QUINTIC은 자동적으로 설정이 되므로, 두가지 키워드를 모두 써 줄 필요는 없습니다.

외삽을 적용한 결과

외삽을 적용한 결과

One Comment

Comments are closed.