netCDF 포맷의 데이터를 읽고 표출하기 [2]

(지난 회 게시물에서 바로 이어집니다)

지난 회 게시물에서 언급했듯이 netCDF 파일에 대하여 NCDF_LIST 명령을 다음과 같이 /VARIABLES와 /VATT 키워드를 함께 사용하면 이 파일 내에 수록되어 있는 모든 변수들의 속성(Variable Attribute)에 대한 정보를 확인할 수 있습니다.

NCDF_LIST, file, /VARIABLES, /VATT

그리고 현재 우리가 다루고 있는 rhum.2003.nc라는 파일에 대하여 이러한 명령을 사용하여 출력된 결과를 보면, 이 파일 내에 수록된 총 5종의 변수들(level, lat, lon, time, rhum) 각각에 대한 좀 더 구체적인 정보를 볼 수 있습니다. 예를 들어 lat이란 변수에 관하여 출력된 내용을 보면 다음과 같습니다.

           1  lat:   FLOAT(73) = FLOAT(lat)

                 0  units: degrees_north

                 1  actual_range:       90.0000

                 1  actual_range:      -90.0000

                 2  long_name: Latitude

이 내용을 보면 lat는 -90부터 +90의 범위를 갖는 73개의 실수값들로 구성된 1차원 배열임을 알 수 있습니다. 위도 격자값들에 해당되며 값들 사이의 간격은 2.5입니다. 경도 격자값들에 해당되는 lon에 대한 정보를 확인해보면 다음과 같습니다.

           2  lon:   FLOAT(144) = FLOAT(lon)

                 0  units: degrees_east

                 1  long_name: Longitude

                 2  actual_range:       0.00000

                 2  actual_range:       357.500

이 내용을 보면 0부터 357.5의 범위를 갖는 144개의 실수값들로 구성되어 있는데, 값들 사이의 간격은 위도와 마찬가지로 2.5입니다. 그리고 level이라는 변수에 대한 정보를 확인해보면 다음과 같습니다.

           0  level: FLOAT(8) = FLOAT(level)

                 0  units: millibar

                 1  actual_range:       1000.00

                 1  actual_range:       300.000

                 2  long_name: Level

                 3  positive: down

                 4  GRIB_id:      100

                 5  GRIB_name: hPa

이 level은 기압고도에 해당되는 8개의 실수값들로 구성된 배열입니다. 최하층은 1000hPa이고 최상층은 300hPa에 해당되는 것으로 확인됩니다. 값들 사이의 간격은 일정하지는 않습니다. 이것은 변수를 추출하여 직접 값들을 확인해보면 알 수 있습니다. 이제 rhum이라는 변수에 대한 정보를 확인해봅시다. 출력된 내용을 보면 다음과 같습니다.

           4  rhum:  INT(144,73,8,365) = INT(lon,lat,level,time)

                 0  long_name: mean Daily relative humidity

                 1  valid_range:      -25.0000

                 1  valid_range:       125.000

                 2  actual_range:       0.00000

                 2  actual_range:       100.000

                 3  units: %

                 4  add_offset:       302.660

                 5  scale_factor:     0.0100000

                 6  missing_value:    32766

실제로는 이보다 좀 더 많지만, 핵심이 되는 부분들만 추려보았습니다. 이 변수는 상대습도에 해당되는 % 단위의 값들로 구성된 4차원 배열임을 알 수 있습니다. 유효 범위는 %이므로 0~100이라고 나와있습니다. 그러면 이 변수부터 먼저 추출해서 직접 확인해봅시다. 다음과 같이 NCDF_VARGET 명령을 사용하여 rhum을 추출하고, HELP와 PRINT를 이용하여 관련 정보들을 출력해봅시다.

NCDF_VARGET, id, 'rhum', rhum

HELP, rhum

PRINT, MIN(rhum), MAX(rhum)

이렇게 하여 출력된 정보는 다음과 같습니다.

RHUM            INT       = Array[144, 73, 8, 365]

  -30266  -20266

HELP에 의하여 출력된 정보는 이미 확인한 것과 다르지 않습니다. 그런데 바로 이어서 출력된 rhum의 최소값 및 최대값은 좀 특이합니다. 분명히 상대습도에 해당되는 % 단위의 값이어야 할텐데, (-)로 큰 값들인데다가 자료형도 정수형입니다. 사실 이 값 자체가 바로 의미있는 데이터값이 되지는 않습니다. 앞서 확인했던 정보에 의하면 add_offset과 scale_factor라는 값들이 각각 302.66과 0.01으로 주어진 것을 확인할 수 있습니다. 따라서 rhum에 대해서는 다음과 같이 0.01을 곱한 다음 302.66을 더해줘야 실제로 의미있는 상대습도 데이터가 된다는 것을 유의해야 합니다.

rhum = rhum*0.01+302.66

HELP, rhum

PRINT, MIN(rhum), MAX(rhum)

이렇게 하여 출력된 내용을 보면 이제 rhum은 0~100의 범위를 갖는 실수값들로 구성된 4차원 배열로 변환되었음을 알 수 있습니다.

RHUM            FLOAT     = Array[144, 73, 8, 365]

      0.00000      100.000

그러면 이제부터는 표출 과정으로 들어가보겠습니다. 이 과정에서는 rhum라는 4차원 배열로부터 특정 시간 및 층에 대한 데이터를 2차원 배열의 형태로 뽑아낸 다음 이 2차원 데이터를 지도상에 중첩 표출하는 작업을 진행하고자 합니다. 먼저 다음과 같이 rhum에 대하여 3, 4번째 차원에 대해서 0, 0이란 값으로 특정지어서 2차원 배열을 추출합니다. 즉 최하층에 해당되는 맨 첫 시간대의 2차원 데이터를 추출하는 셈입니다.

data = rhum[*, *, 0, 0]

HELP, data

이렇게 하면 data는 144x73의 구조를 갖는 2차원 배열임을 알 수 있습니다. 일단 이 2차원 데이터의 모습을 확인해보기 위하여 다음과 같이 data의 가로/세로 크기 비율과 유사하면서 크기가 더 큰 그래픽 창을 띄운 후 IMAGE 함수를 사용하여 data를 간단하게 이미지의 형태로 먼저 표출해봅시다.

win1 = WINDOW(DIMENSIONS=[144, 73]*6, /NO_TOOLBAR)

im1 = IMAGE(data, RGB_TABLE=34, MARGIN=0, /CURRENT)

이렇게 표출된 모습은 다음 그림과 같습니다.

앞서 이 데이터의 경도 및 위도 범위를 확인해봤듯이, 이 데이터는 경도로는 0~357.5 그리고 위도로는 -90~90의 범위를 갖습니다. 즉 경위도 전 범위를 커버하는 데이터임을 알 수 있습니다. 따라서 이 데이터를 중첩할 바탕 지도를 그린다면 그 지도는 당연히 전세계 지도가 되어야 합니다. 그래서 이번에는 그래픽 창을 별도로 하나 더 띄운 후 다음과 같이 바탕 지도를 표출하고 이미지를 중첩하고 대륙경계선을 그리는 과정을 진행해봅시다. 이 과정은 NG 체계에서 지도상에 2차원 데이터를 중첩하는 전형적인 과정으로 제가 이 블로그에서 여러번 다룬 바 있는 방법입니다.

win2 = WINDOW(DIMENSIONS=[900, 500], /NO_TOOLBAR)

m = MAP('Geographic', MARGIN=[0.05, 0.20, 0.03, 0.05], $

  ASPECT_RATIO=0, CLIP=0, FONT_SIZE=10, /CURRENT)

m.MapGrid.LABEL_POSITION = 0

im2 = IMAGE(data, IMAGE_LOCATION=[-180, -90], IMAGE_DIMENSION=[360, 180], $

  GRID_UNITS=2, RGB_TABLE=34, ASPECT_RATIO=0, /OVERPLOT)

mc = MAPCONTINENTS(THICK=2)

여기서는 2차원 데이터가 지구 전영역을 커버한다고 가정하고 IMAGE 함수의 IMAGE_LOCATION 및 IMAGE_DIMENSION 속성들을 사용하였습니다. 이미지의 시작점은 경도 -180도 및 위도 -90도에 해당되는 지점이 되고, 여기서부터 경도 방향으로 360도 그리고 위도 방향으로 180도 만큼을 커버하는 이미지가 되도록 중첩한 것입니다. 결과는 다음 그림과 같습니다.

이 그림을 보면 제법 그럴싸한 모습으로 표출이 된 것처럼 보입니다. 하지만 이 그림은 잘못된 그림입니다. 표출 과정 특히 IMAGE 함수를 사용한 내용에서 이미지가 커버하는 범위를 설정하는데 있어서 매우 중대한 실수를 범하고 있습니다. 사실 이 과정에 있어서는 앞서 기술했던 변수 lon, lat를 다음과 같이 추출하여 함께 사용해야 합니다. 왜냐하면 이 변수들이 2차원 데이터를 구성하는 격자들의 경도 및 위도 방향 위치에 관한 정보를 갖고 있기 때문입니다.

NCDF_VARGET, id, 'lon', lon

NCDF_VARGET, id, 'lat', lat

그래서 이렇게 추출된 lon, lat를 IMAGE 함수에서 다음과 같은 방식으로 사용해야 합니다. 위의 예제코드에서 IMAGE 함수가 사용된 부분만 다음과 같은 내용으로 대체합시다.

im2 = IMAGE(data, lon, lat, GRID_UNITS=2, RGB_TABLE=34, $

  ASPECT_RATIO=0, /OVERPLOT)

이렇게 하여 표출된 그림은 다음과 같습니다.

이 그림을 보면 경도가 동쪽에 해당되는 영역은 제대로 이미지가 중첩되어 있습니다. 그런데 서쪽 영역은 전혀 이미지가 보이지 않고 있습니다. 따라서 표출 과정이 아직도 제대로 진행되지 않고 있음을 알 수 있습니다. 이것은 이유가 있습니다. 가장 먼저 lon이란 배열에 있는 경도값들이 문제입니다. 앞서 확인했듯이 lon은 0~357.5의 범위의 144개의 값들을 갖고 있습니다. 그런데 IDL에서 지도상에 데이터나 annotation 등을 중첩 표시하고자 할 경우에는 경도의 좌표값이 반드시 -180~+180의 범위내에 있어야 합니다. 이 범위를 벗어난 값을 사용할 경우 제대로 표시가 안됩니다. 따라서 경도값들을 담은 lon에 대하여 일제히 180만큼 빼줘야 합니다. 그러면 경도 방향의 범위가 0~357.5가 아닌 -180~177.5로 인식되어 중첩이 정상적으로 됩니다. 그리고 이 상황은 경도 방향으로 데이터가 왼쪽으로 180도만큼(격자 크기로는 144의 절반인 72만큼) 이동하는 것이나 마찬가지이므로, data에 대해서 다음과 같이 SHIFT 함수를 사용하여 X축 방향으로 그만큼 이동을 시켜줘야 합니다. 따라서 IMAGE 함수를 사용하는 부분은 다음과 같은 내용이 되어야 맞습니다.

im2 = IMAGE(SHIFT(data, 72, 0), lon-180, lat, GRID_UNITS=2, $

  RGB_TABLE=34, ASPECT_RATIO=0, /OVERPLOT)

그리고 그래픽 창의 하단에 컬러바도 함께 표시하는 것이 더 좋을 것 같아서, 다음과 같이 COLORBAR 함수를 사용하여 컬러바까지 표출하였습니다.

cb = COLORBAR(TARGET=im2, TITLE='Relative Humidity', FONT_SIZE=10, $

  POSITION=[0.2, 0.07, 0.8, 0.10], /BORDER)

이렇게 하여 표출된 그림은 다음과 같습니다.

이 정도면 중첩 표출이 제대로 된 모습이라 할 수 있겠습니다. 그런데 이 그림을 아까 맨 처음에 지도 없이 이미지만 봤던 모습과 비교해보면 엄청난 차이점이 있습니다. 즉 win1과 win2 그래픽 창에 표출된 두 그림을 서로 비교해보면 위도 방향 즉 세로 방향이 서로 반대인 것을 확인할 수 있습니다. 이러한 이유에 대해서는 위도 방향 격자점 좌표들을 담고 있는 lat라는 배열의 값들을 출력해보면 알 수 있습니다. 출력을 해보면 값들이 -90으로 시작되는 것이 아니라 +90으로 시작되어 -90까지 -2.5의 간격으로 감소하는 방향으로 들어가 있습니다. 그래서 rhum에서도 위도에 해당되는 차원에서는 값들의 순서가 이 방향에 맞도록 되어 있기 때문입니다. 물론 IMAGE 함수에서는 격자점들의 좌표값 하나하나에 맞게 이미지 격자를 표시하기 때문에, lat 배열 내 값의 순서 자체에는 영향을 받지 않습니다. 어쨌든 전체 데이터를 구성하는 정보들을 면밀히 파악하지 않으면, 표출 및 분석 과정에서 뜻하지 않은 실수가 나올 가능성이 항상 있음을 다시 한번 알 수 있습니다.

그리고 한가지 더 언급한다면, 위의 IMAGE 함수가 사용된 내용에서는 다음과 같이 MIN_VALUE, MAX_VALUE 속성을 추가적으로 사용하는 것도 좋습니다.

im2 = IMAGE(SHIFT(data, 72, 0), lon-180, lat, GRID_UNITS=2, $

  MIN_VALUE=0, MAX_VALUE=100, RGB_TABLE=34, ASPECT_RATIO=0, /OVERPLOT)

이렇게 하면 이미지를 표출할 때 값의 스케일을 항상 0~100으로 고정하여 맞추게 됩니다. 이것은 혹시라도 다른 시간이나 층의 데이터를 표출할 때 그 데이터에서 최소 및 최대값이 0~100이 아닐 경우라 하더라도 그림상에서는 항상 0~100의 범위를 가정하여 컬러테이블의 색상들을 부여하도록 하기 위함입니다. 이렇게 하면, 값의 범위가 다른 데이터가 표출될 경우에도 서로간의 동등한 비교가 가능합니다. 이러한 설정을 안할 경우 IMAGE 함수에서는 항상 그 데이터의 최소 및 최대값을 직접 구하여 그대로 색상으로 반영합니다. 즉 데이터 자체의 값 범위가 달라질 경우에는 거기에 연동하여 컬러테이블의 모든 색상들을 사용하는 것이 디폴트 설정이기 때문입니다. 예를 들어, 값의범위가 0~100인 데이터이든 아니면 0~50인 데이터이든간에 컬러테이블의 가장 뒤쪽 색상이 최대값에 대응된다면, 눈으로 봐서는 값 범위의 차이를 느낄 수가 없습니다. 따라서 기본적으로 이 데이터가 상대습도라는 특성상 0~100이 유의미한 범위라고 한다면, 표출 대상 데이터의 값 범위가 어떻든간에 0이 가장 앞쪽 색상으로 그리고 100이 가장 뒤쪽 색상으로 표시되도록 고정되는 것이 더 필요할 수도 있다고 봅니다.

지금까지 2회에 걸쳐 netCDF 파일의 내용을 파악하여 추출하고 그 데이터의 특성에 맞는 표출 과정도 함께 살펴보았습니다. 당연한 얘기이지만, 파일 형식이 netCDF라 하더라도 실제 데이터 파일마다 그 안에 담겨 있는 데이터들은 천차만별이기 때문에, 파일 내에 수록된 각종 정보들을 정확히 파악하여 데이터의 추출 뿐 아니라 이후의 처리 및 분석 과정에 있어서 제대로 활용될 수 있도록 하는 것이 중요하다는 점을 다시 한번 느끼게 됩니다.