벡터장(Vector Field)과 입자 추적(Particle Tracing) [1]

이번에 제가 소개하고자 하는 주제는 벡터장과 입자 추적입니다. 즉 바람이나 해류와 같은 유체의 속도 값들이 격자형으로 분포하는 데이터가 있을 때, 특정 위치에 있던 무질량의 입자가 유체의 흐름을 타고 어떻게 이동하는가를 계산하는 작업에 관한 것입니다. 규칙 격자의 형태로 분포하는 벡터 자료를 벡터장(Vector Field)의 형태로 표출하기 위하여 VECTOR 함수를 사용하는 방법 및 예제에 관하여 제가 예전에 이 블로그에서 소개했던 바 있습니다. 해당 게시물들은 다음과 같습니다.

 

http://blog.daum.net/swrush/437

http://blog.daum.net/swrush/438

 

이 게시물들에서는 VECTOR 함수를 사용하여 벡터장을 표출하는 방법에 관한 기능적 소개에 중점을 두었었는데요. 이번 주제에서는 벡터장의 표출 뿐 아니라 그 안에서의 입자의 이동을 추적하는 방법에 관한 내용까지 다뤄보고자 합니다. 이 과정에서는 주로 VECTOR 함수 뿐만 아니라 PARTICLE_TRACE 프로시저 및 STREAMLINE 함수도 함께 사용될 예정입니다. 따라서 내용이 좀 많을 것 같으므로 회차를 나누어서 소개하도록 하겠습니다. 그러면 먼저 예제 데이터부터 가져옵시다. 이번에 사용할 예제 데이터가 수록된 파일은 IDL의 설치와 함께 딸려오는 globalwinds.dat라는 파일입니다. 이 파일로부터 다음과 같이 u, v, x, y 배열들을 얻으면 됩니다.

 

file = FILEPATH('globalwinds.dat', SUBDIR=['examples','data'])

RESTORE, file

HELP, u, v, x, y

 

그런데 위의 내용이 약간 이상하다고 생각하실 수도 있습니다. 파일의 확장자는 .dat인데 막상 이 파일을 열 때 사용한 명령은 RESTORE이기 때문입니다. 사실 RESTORE 명령은 .sav 확장자의 파일에 담긴 정보들을 추출할 때 사용하는 것이 일반적입니다. 그리고 .sav 파일은 IDL에서 생성되고 IDL끼리만 통용되는 데이터 저장용 형식입니다. 사실 globalwinds.dat 파일 자체는 실제로는 .sav 형식이 맞습니다. 다만 무슨 이유에서인지 제작자가 파일명의 확장자를 그냥 .dat로 바꿔놓은 것 같습니다. 하여간 본질적으로는 .sav 파일이기 때문에 위와 같이 RESTORE 명령을 사용하여 열면 됩니다. 그러면 u, v, x, y 4종의 배열들을 추출하게 됩니다. HELP에 의하여 출력된 각 배열에 관한 정보들은 다음과 같습니다.

 

U               FLOAT     = Array[128, 64]

V               FLOAT     = Array[128, 64]

X               FLOAT     = Array[128]

Y               FLOAT     = Array[64]

 

이 데이터는 파일명에서 짐작할 수 있듯이 바람(Wind) 자료인데, 구체적으로는 글로벌 맵(Global Map) 상에서 [128, 64]의 2차원 격자 형태로 분포하는 바람 벡터 자료입니다. 풍향 및 풍속에 대한 U-성분 및 V-성분 값들은 각각 [128, 64]의 2차원 구조를 갖는 u 및 v 배열에 수록되어 있습니다. 그리고 X 방향 격자들의 경도 좌표값들 128개 및 Y 방향 격자들의 위도 좌표값들 64개가 1차원의 형태로 각각 x 및 y 배열에 수록되어 있습니다. 실제로 배열 x의 값들은 대략 -180~+180 그리고 배열 y의 값들은 -90~+90 정도의 범위에 있습니다.

 

그러면 이 데이터를 벡터장의 형태로 표출해봅시다. 먼저 바탕 지도를 표출하고 그 위에 벡터장이 중첩 표출되도록 할 것입니다. 이를 위하여 우선 다음과 같이 MAP, MAPCONTINENTS 등의 함수들을 사용하여 바탕이 될 글로벌 맵을 표출합니다. 그 과정은 다음과 같습니다.

 

win = WINDOW(DIMENSIONS=[800, 480], /NO_TOOLBAR)

m = MAP('Geographic', LIMIT=limit, FILL_COLOR='light blue', $

  MARGIN=0.05, /CURRENT)

mc = MAPCONTINENTS(FILL_COLOR='gold')

m.MapGrid.LABEL_POSITION = 0

m.MapGrid.GRID_LONGITUDE = 30

m.MapGrid.LINESTYLE = 2

m.MapGrid.HORIZON_LINESTYLE = 0

 

이 상태에서 표출된 모습은 다음 그림과 같습니다.

 

 

이제 VECTOR 함수를 사용하여 이 지도 위에 벡터장을 중첩하여 표출해봅시다. 그 과정은 다음과 같습니다.

 

vec = VECTOR(u, v, x, y, AXIS_STYLE=2, ARROW_THICK=3, $

  HEAD_SIZE=0.5, /HEAD_PROPORTIONAL, LENGTH_SCALE=0.5, $

  COLOR='green', /OVERPLOT)

 

여기서는 벡터장의 벡터들이 초록색의 화살표로 표시되도록 하였습니다. 그 결과는 다음 그림과 같습니다.

 

 

일단 표출 자체는 제대로 된 상태입니다. 다만 문제는 벡터 화살표들의 가독성이 좀 떨어진다는 것인데, 아무래도 많은 수의 벡터들이 오밀조밀 모여있는 상태이기 때문에 어느 정도 한계는 있는 것 같습니다. 이러한 가독성 문제를 극복하는 방법은 몇가지가 있겠지만, 어차피 나중에는 지도의 경위도 범위를 더 좁혀서 표출할 예정이고 그러면 가독성은 더 나아질 것이기 때문에 일단 현 시점에서는 이대로 두겠습니다.

 

그러면 이제부터는 입자 추적(Particle Tracing)이라는 작업을 진행해보겠습니다. 입자 추적이라는 것은 위와 같은 벡터장 내에서 어떤 입자가 특정 위치에 있다고 가정하고 이 입자가 유체의 흐름에 따라 이동하는 궤적을 계산하는 작업입니다. 다소 복잡해 보이는 작업이긴 한데, IDL에서는 이를 위한 기능이 탑재되어 있습니다. 바로 PARTICLE_TRACE 프로시저입니다. IDL 도움말에 의하면 그 사용 문법은 다음과 같습니다.

 

PARTICLE_TRACE, Data, Seeds, Verts, Conn, [각종 키워드들]

 

여기서 첫번째 및 두번째 인수들인 data, seeds은 input 인수이고 세번째 및 네번째 인수들인 verts, conn은 output 인수입니다. 즉 data, seeds 항목들을 인수로 제공해주면 그 결과를 verts, conn 항목들을 통하여 전달받게 됩니다. 첫번째 인수 data의 경우는 반드시 [2, nx, ny]의 구조를 갖는 3차원 배열의 형태가 되어야 합니다. 따라서 이미 우리가 갖고 있는 u, v 배열들을 다음과 같이 재구성할 필요가 있습니다.

 

data = FLTARR(2, 128, 64)

data[0, *, *] = u

data[1, *, *] = v

 

두번째 인수인 seeds의 경우는 입자의 시작 위치 좌표 인덱스에 해당됩니다. 따라서 다음과 같이 정의해봅시다.

 

x_ind0 = 92

y_ind0 = 32

seed = [x_ind0, y_ind0]

PRINT, 'lon/lat of the starting point :', x[x_ind0], y[y_ind0]

 

여기서는 시작 위치를 정의하는데 있어서 X 방향 인덱스가 92 그리고 Y 방향 인덱스가 32가 되도록 하였습니다. 다만 유의해야 할 점은, 92와 32는 좌표값 자체가 아니라 "인덱스"라는 것입니다. 따라서 시작 위치의 경도 좌표값은 배열 x를 x_ind0로 인덱싱해야 확인이 가능하고, 위도 좌표값은 배열 y를 y_ind0로 인덱싱해야 확인이 가능합니다. 여기서는 PRINT 명령에 의하여 이러한 내용들도 출력하도록 하였는데, 출력된 내용을 보면 다음과 같습니다.

 

lon/lat of the starting point :      78.7500      1.00000

 

즉 인덱스 [92, 32]에 해당되는 실제 좌표값은 동경 78.75도 및 북위 1.0입니다. 이것은 바로 입자의 출발 위치의 경도 및 위도 좌표입니다. 즉 이 지점에서 입자가 출발하는 것으로 가정한 셈입니다. 이제 이 인수들을 투입하여 PARTICLE_TRACE 프로시저를 다음과 같이 실행합시다. 그러면 바로 결과를 얻게 됩니다.

 

n_iter = 30

PARTICLE_TRACE, data, seed, verts, conn, $

  MAX_ITERATIONS=n_iter

HELP, verts, conn

 

이와 같이 input 인수들인 data, seed를 PARTICLE_TRACE 프로시저에 투입하여 그 계산 결과를 verts, conn으로 전달받게 됩니다. 그리고 MAX_ITERATIONS라는 키워드를 사용하였는데, 이것은 입자 추적을 몇 단계에 걸쳐서 진행할 것인가를 설정하는 역할을 합니다. 여기서는 30으로 설정하였습니다. HELP에 의하여 출력된 내용을 보면 다음과 같습니다.

 

VERTS           FLOAT     = Array[2, 30]

CONN            LONG64    = Array[31]

 

이 내용을 보면 verts는 [2, 30]의 구조를 갖는 2차원 배열인데, 실제로는 입자의 이동 궤적을 구성하는 30개의 점들의 [X, Y] 좌표 인덱스들로 구성됩니다. 점들의 갯수가 30개인 것은 앞서 MAX_ITERATIONS 키워드에 30이란 값을 설정하였기 때문입니다. 즉 입자의 이동 궤적을 몇 회에 걸쳐서 추적할 것인가를 MAX_ITERATIONS 키워드에서 설정한다고 보면 됩니다. 그리고 conn의 경우는 31개의 값들로 구성된 1차원 배열인데, 이것은 30개의 점들의 연결성(connectivity)을 나타내는 정보입니다. 이 정보는 지금과 같이 입자를 한 개만 가정할 경우에는 별로 필요가 없습니다. 하지만 다수의 입자들을 한꺼번에 정의할 경우에는 필요합니다. 일단 우리는 하나의 입자만 가정한 상태이기 때문에 conn에 대해서는 더 이상 신경쓰지 않기로 하겠습니다. 그리고 이후의 작업의 편의성을 위하여, verts 배열로부터 X 및 Y 방향 데이터를 각각 분리해두는 것이 좋을 것 같아서 다음과 같은 처리 과정을 추가합니다.

 

x_inds = REFORM(verts[0, *])

y_inds = REFORM(verts[1, *])

PRINT, 'lon/lat of the ending point :', x[x_inds[-1]], y[y_inds[-1]]

 

이렇게 하면 입자의 이동 궤적에 대한 X 방향 인덱스 좌표들만으로 구성된 1차원 배열 x_inds, 그리고 Y 방향 인덱스 좌표들만으로 구성된 1차원 배열 y_inds를 정의할 수 있습니다. 그리고 각 배열에 대하여 마지막 인덱스 값에 대한 좌표값도 출력하였습니다. 즉 이동 궤적의 최종 도달점의 경도 및 위도 좌표를 확인하기 위한 것입니다. 출력된 내용은 다음과 같습니다.

 

lon/lat of the ending point :     -5.62500     -4.62500

 

이제 바로 위에서 얻은 x_inds, y_inds 배열들을 이용하여 입자의 이동 궤적을 중첩하여 표출해봅시다. 이를 위해서는 다음과 같이 POLYLINE 함수를 사용하면 됩니다. 그리고 이왕이면 궤적의 최종 도달점에 화살표를 표시해주고 싶은데, POLYLINE 함수로는 궤적상의 점들을 선으로 이어주는 것만 가능합니다. 따라서 화살표 표시를 위하여 ARROW 함수도 추가적으로 사용하겠습니다. 이 과정은 다음과 같습니다.

 

line = POLYLINE(x[x_inds], y[y_inds], THICK=3, COLOR='red', /DATA)

arr = ARROW(x[x_inds[-2:-1]], y[y_inds[-2:-1]], $

  THICK=3, HEAD_SIZE=0.8, COLOR='red', /DATA)

 

그리고 표출된 그림은 다음과 같습니다.

 

 

이 그림을 보면 입자가 출발점으로부터 어떤 경로를 통해서 이동하는가를 가늠해볼 수 있습니다. 물론 앞서 이미 언급하였듯이 글로벌 맵 상에서 바람 벡터들이 촘촘하게 붙어있는 상황이기 때문에, 위와 같은 이동 궤적이 정말 그럴듯한 결과인지를 시각적으로 확인하는 것이 다소 어려운 상태입니다. 따라서 이번에는 지도의 경위도 범위를 좁혀봅시다. 이를 위해서는 앞서 바탕 지도를 표출하는 과정에 대하여 다음과 같이 약간의 수정만 해주면 됩니다.

 

win = WINDOW(DIMENSIONS=[800, 480], /NO_TOOLBAR)

limit = [-20, -10, 30, 90]

m = MAP('Geographic', LIMIT=limit, FILL_COLOR='light blue', $

  MARGIN=0.05, /CURRENT)

mc = MAPCONTINENTS(FILL_COLOR='gold');, /HIRES)

m.MapGrid.LABEL_POSITION = 0

;m.MapGrid.GRID_LONGITUDE = 30

m.MapGrid.LINESTYLE = 2

m.MapGrid.HORIZON_LINESTYLE = 0

 

여기서 볼드체로 표시된 부분들이 바로 새롭게 수정되어야 할 부분들입니다. MAP 함수에서 LIMIT 키워드를 사용하여 경위도 범위가 제한되도록 하였고, 경도 눈금값들의 간격이 30으로 설정되지 않도록 그 라인은 주석처리하였습니다. 이러한 수정사항들만 반영하여 전체적으로 다시 실행해보면 됩니다. 이렇게 표출된 그림은 다음과 같습니다.

 

 

이와 같이 입자의 이동 궤적을 중점적으로 더 크게 나타내기 위하여 경위도 범위를 제한하여 표출했더니 이제는 주변의 바람 벡터들이 더 잘 구분되어 보입니다. 이 모습을 보면 아마도 입자의 이동 궤적이 주변의 바람 벡터들에 의하여 그럴듯하게 산출되었음을 확인할 수 있습니다. 입자의 이동 궤적을 산출하기 위하여 사용한 PARTICLE_TRACE 프로시저에서는 MAX_ITERATIONS 외에도 INTEGRATION, MAX_STEPSIZE 등의 키워드들도 있습니다. 원래 PARTICLE_TRACE 프로시저에서 내부적으로 사용하는 알고리즘은 2nd-order Runge-Kutta 기법인데, 이 경우는INTEGRATION 키워드가 디폴트 값인 0인 경우에 해당됩니다. 만약 이 키워드의 값을 1로 설정하면 4th-order Runge-Kutta 기법이 적용됩니다. 두 기법 사이의 차이는 상황에 따라 클 수도 있고 작을 수도 있습니다. 위의 예제에서 PARTICLE_TRACE 프로시저에 INTEGRATION 키워드를 추가하여 그 값을 1로 설정해서 얻은 결과는 다음 그림과 같습니다.

 

 

이 그림에서 나타나는 궤적은 앞서 디폴트 설정으로 얻었던 결과와는 살짝 다르긴 합니다. 어쨌든 입자의 시작 위치, iteration의 갯수, 세부 기법, 스텝 크기 등을 해당 키워드에 의하여 조정해가면서 여러가지 다양한 결과들을 얻을 수 있다는 점을 참조하시기 바랍니다. 오늘은 여기까지 하고, 다음 회차에서는 이 내용의 연장으로서 다수의 입자들의 이동 궤적을 동시에 추적하는 방법에 관하여 다뤄보고자 합니다. 이 과정에서는 STREAMLINE 함수를 사용하게 될 예정입니다.