GRIB 포맷의 파일을 읽고 데이터를 표출하기 [2]

* 지난 회 게시물에서 바로 이어집니다.

지난 회 게시물에서는 확장자가 gb2인 GRIB2 포맷의 파일을 읽어서 파일 내에 수록된 레코드 및 세부 정보 등을 확인하는 작업을 진행하였습니다. 여기서는 '2 meter temperature'에 해당되는 21번 레코드로부터 Ordered Hash 타입의 데이터를 data라는 이름으로 추출하였고, 이렇게 추출된 data 내부에 204개의 key-value 쌍들이 존재함을 확인하였습니다. 아직은 이 204개의 키(Key)들 중 몇몇 항목들만 조금 살펴본 상태입니다. 이제 오늘은 표출 작업을 위하여 필요한 키(Key)들을 골라내고 각각에 대한 실제 데이터를 추출하는 작업부터 진행해보겠습니다. 일단 제가 보기에 후보가 될 만한 항목들을 다음과 같은 방식으로 추출해보았습니다.

dstr = STRTRIM(STRING(data['dataDate']), 2)

dname = data['name']

lon0 = data['longitudeOfFirstGridPointInDegrees']

lat0 = data['latitudeOfFirstGridPointInDegrees']

dx = data['DxInMetres']

dy = data['DyInMetres']

nx = data['Nx']

ny = data['Ny']

lonc = data['LoVInDegrees']

latc = data['LaDInDegrees']

latsp1 = data['Latin1InDegrees']

latsp2 = data['Latin2InDegrees']

lons = data['longitudes']

lats = data['latitudes']

values = data['values']

여기서는 data라는 Ordered Hash 타입의 데이터로부터 각각의 키(Key) 명칭으로 인덱싱(Indexing)을 하는 방식으로 데이터를 추출하였습니다. 지난 회 게시물 막바지에서 언급했던 'dataDate', 'name', 'Nx', 'Ny' 외에도 꽤 많은 종류의 데이터들을 추출한 것을 볼 수 있습니다. 여기서 2차원 격자 데이터에 해당되는 것은 'vaules' 키로 인덱싱한 values입니다. 그런데 이 values에 대하여 HELP 명령을 사용하여 그 실체를 확인해보면 결과는 다음과 같습니다.

HELP, values

VALUES          DOUBLE    = Array[470162]

그런데 제가 분명히 "2차원" 격자 데이터라고 했지만, 이와 같이 values의 실제 구조는 2차원이 아닌 1차원 배열임을 확인할 수 있습니다. 사실 이러한 상황에 대한 실마리는 nx,ny에서 찾을 수 있는데요. 지난 회차에서 이미 언급했지만 여기서도 다음과 같이 확인해볼 수 있습니다.

PRINT, nx, ny

이 때 nx, ny는 각각 602, 781로 나옵니다. 이 두 숫자를 곱하면 바로 470162가 나옵니다. 따라서 values는 겉모습 자체는 1차원이긴 한데, 내부적으로는 602x781 구조의 배열에 포함된 모든 값들을 1차원 배열의 형태로 갖고 있는 셈입니다. 따라서 향후의 표출 작업을 위해서는 다음과 같이 values에 대하여 REFORM 함수를 사용하여 2차원으로 변형된 배열을 만들어두는 것이 필요합니다.

img = REFORM(values, nx, ny)

HELP, img

IMG             DOUBLE    = Array[602, 781]

사실 이러한 상황은 values 뿐만 아니라 lons, lats의 경우도 마찬가지입니다. 참고로 lons는 470162개의 데이터 포인트들의 경도 좌표값들로 구성된 1차원 배열이고, lats는 위도에 해당되는 동일한 형태의 배열입니다. 그래서 values, lons, lats를 그대로 SCATTERPLOT 함수에 투입하는 표출 방식도 생각해 볼 수 있겠으나, 이런 방식의 작업은 매우 비효율적일 가능성이 큽니다. 47만여개나 되는 포인트들을 하나하나 찍는 표출 방식 자체가 시간이 꽤 걸릴 뿐더러 메모리 부담도 상당할 가능성이 높기 때문입니다. 따라서 그 대신 602x781 구조의 2차원 배열인 상태로 IMAGE 함수를 사용하여 표출하는 것이 여러모로 더 효율적입니다.

또한 이 데이터와 관련된 중요한 특성은, 현재 우리가 다루고 있는 img라는 2차원 배열의 기반 격자 구조가 경도/위도 기반이 아니라 거리 기반이라는 것입니다. 이것은 LDAPS 데이터의 특성 자체가 LCC(Lambert Conformal Conic) 투영법의 지도 상에서 거리 기반의 격자 분포를 하고 있다는 점을 근거로 한 것입니다. 따라서 표출 과정에서도 이러한 점을 반영해야 하며, 이러한 표출을 위해서 필요한 인자값들이 바로 위에서 추출한 lon0, lat0, dx, dy, lonc, latc, latsp1, latsp2입니다. 각각의 의미는 다음과 같습니다.

lon0 : 2차원 격자 데이터의 좌측 하단 꼭지점의 경도 좌표

lat0 : 2차원 격자 데이터의 좌측 하단 꼭지점의 위도 좌표

dx : 격자와 격자 사이의 X 방향 거리

dy : 격자와 격자 사이의 Y 방향 거리

lonc : LCC 투영법으로 지도 표출시 중심 경도

latc : LCC 투영법으로 지도 표출시 중심 위도

latsp1 : Standard Parallel 1 위도 (LCC 투영법에서 필요)

latsp2 : Standard Parallel 2 위도 (LCC 투영법에서 필요)

그래서 이러한 인자들을 사용하여 지도를 LCC 투영법으로 먼저 표출해야 합니다. 이 과정은 다음과 같습니다.

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

limit = [32, 120, 44, 134]

m = MAP('Lambert Conformal Conic', LIMIT=limit, $

  STANDARD_PAR1=latsp1, STANDARD_PAR2=latsp2, $

  CENTER_LONGITUDE=lonc, CENTER_LATITUDE=latc, $

  HORIZON_LINESTYLE=0, HORIZON_THICK=2, $

  ASPECT_RATIO=0, POSITION=[0.06, 0.21, 0.94, 0.92], $

  Font_Size=11, CLIP=0, /CURRENT)

m.MapGrid.GRID_LONGITUDE = 2

m.MapGrid.GRID_LATITUDE = 2

m.MapGrid.LINESTYLE = 2

m.MapGrid.LABEL_POSITION = 0

이 과정에서는 LCC 투영법으로 지도를 표출하는데 있어서 lonc, latc, latsp1, latsp2가 핵심적으로 사용되었습니다. 그리고 limit라고 되어 있는 것은 지도를 표출할 때 경위도 범위에 해당되는데, 이 범위는 나중에 2차원 격자 자료의 중첩까지 고려하여 적절하게 정해보았습니다. 그러면 이제는 지도상에 중첩할 2차원 격자 자료에 관한 처리 작업이 필요한데, 앞서 values라는 1차원 배열을 2차원으로 변형한 img라는 배열을 이미 만들어놓은 상태입니다. 하지만 여기서 필요한 작업이 하나 더 있습니다. 이 2차원 데이터가 위치하는 영역에 해당되는 X 방향 및 Y 방향 격자들의 좌표값들로 구성된 배열의 생성입니다. 이 과정은 다음과 같습니다.

xy = m.MapForward(lon0, lat0)

HELP, xy

PRINT, xy

xs = xy[0]+FINDGEN(nx)*dx

ys = xy[1]+FINDGEN(ny)*dy

여기서는 먼저 MapForward라는 메서드를 사용하여 2차원 격자 자료의 좌측 하단 꼭지점의 경위도 좌표인 (lon0, lat0)를 거리 기반의 좌표로 환산하였습니다. 환산된 결과는 xy라는 배열로 생성되었는데, HELP와 PRINT로 출력된 xy에 관한 정보는 다음과 같습니다.

XY              DOUBLE    = Array[2]

      -2761659.0      -2594770.3

꼭지점의 위치에 대한 거리 기반의 XY 좌표값들이 이렇게 나옵니다. 이렇게 (-)값들로 나온 이유는 거리 기반 XY 좌표의 경우 원점은 지도의 중심점인 (lonc, latc) 위치가 되기 때문입니다. 그리고 바로 이어서 xs, ys라는 두개의 1차원 배열을 생성하였습니다. 거리 기반의 X 방향 격자점 좌표값들끼리 모은 1차원 배열이 xs이며, 각 격자간의 거리 간격의 값이 담긴 dx라는 변수도 함께 사용되었음을 알 수 있습니다. Y 방향에 대해서도 유사한 요령으로 ys라는 배열을 생성하였습니다. 이렇게 img, xs, ys가 모두 있어야 2차원 격자 자료의 중첩을 제대로 할 수 있게 됩니다. 그래서 표출 작업을 바로 진행하면 좋겠지만, 조금만 참으시고 일단 img에 대하여 다음과 같은 처리 작업 하나만 더 해줍시다.

ww = WHERE(img GT 9000, count)

IF count NE 0 THEN img[ww] = !values.f_nan

이러한 작업을 굳이 하는 이유는 2차원 격자 자료 배열인 img 안에는 9999라는 비정상적인 값들도 포함되어 있기 때문입니다. 그래서 이러한 값들은 데이터가 없는 격자로 간주하여, 이러한 격자들에 대해서는 의도적으로 NaN(Not a Number)값으로 대체하였습니다. 이제 다음과 같이 격자 자료를 지도상에 중첩하여 표출해봅시다.

ctnum = COLORTABLE(74, /REVERSE)

im = IMAGE(img, xs, ys, RGB_TABLE=ctnum, GRID_UNITS=1, /OVERPLOT)

mc = MAPCONTINENTS(/HIRES)

cb = COLORBAR(TARGET=im, TITLE='Value', FONT_SIZE=12, $

  TAPER=0, POSITION=[0.2, 0.07, 0.8, 0.10])

tstr = dname + ' (' + dstr + ')'

tx = TEXT(0.5, 0.96, tstr, ALIGNMENT=0.5, FONT_SIZE=16, COLOR='blue', /NORMAL)

여기서는 IMAGE 함수를 사용하여 2차원 데이터의 중첩을 하였습니다. IMAGE 함수에서 img, xs, ys가 차례로 인수로 사용된 것을 주목해야 합니다. 또한 GRID_UNITS라는 키워드의 값이 1로 설정되어 있는데, 이것은 격자들의 좌표 기반이 경위도가 아니라 거리이기 때문에 필요한 설정입니다. 그 외에도 이미지를 컬러로 표출하기 위하여 IDL에 내장되어 있는 74번 컬러테이블의 색상 순서를 뒤집은 컬러테이블을 COLORTABLE 함수로 생성하여 사용하였고, COLORBAR 함수를 사용하여 적절한 위치에 컬러바를 삽입하였습니다. 또한 표출 데이터의 종류 및 날짜 정보를 문자로 삽입하기 위하여 TEXT 함수를 사용하여 상단에 적절한 문구를 삽입하였습니다. 이러한 과정들을 거쳐 표출된 그림은 다음과 같습니다.

지금 여기서 표출된 2차원 격자 자료는 이미 언급했듯이 "2 meter temperature"입니다. 그런데 값의 단위는 실제로는 절대온도(K)로 되어 있습니다. 만약 섭씨온도 단위로 바꾸려면 절대온도 기준값인 273.15를 빼줘야 합니다. 그래서 img에 대해서만 다음과 같이 처리할 수도 있습니다.

IF rnum EQ 21 THEN img = img - !const.T0

여기서는 IDL이 자체적으로 내장하고 있는 절대온도 기준값에 해당되는 !const.T0라는 시스템 상수를 사용하였습니다. 물론 273.15라고 직접 숫자로 적어줘도 됩니다. 그리고 굳이 IF문을 사용한 이유는 레코드 번호 21에 해당되는 "2 meter temperature"에 대해서만 이러한 처리가 필요하기 때문이라고 판단했기 때문입니다. 이렇게 한 다음에 그림을 다시 그려보면 그 모습은 다음과 같습니다. 하단의 컬러테이블의 값 범위가 변한 것을 확인할 수있습니다.

지난 회 게시물에서 이미 확인했듯이 지금 사용중인 예제용 gb2 파일 내에는 총 136개의 레코드들이 있습니다. 그 중 21번 레코드에 해당되는 "2 meter temperature" 데이터를 지도상에 중첩 표출하는 작업을 앞서 여러가지 과정들을 거쳐 진행해보았는데요. 그 외에도 하나 주목해볼만한 것이 133번 레코드에 해당되는 "land-sea mask"입니다. 이것은 육지와 바다를 구분지어 보여주는 역할을 하는데, 다른 데이터들을 표출하는데 있어서 육지와 바다 부분에 대한 구분이 제대로 된 상태인지를 확인하는 차원에서 한번 표출해보는 것도 좋을 것 같습니다. 이전 게시물에 있던 예제 코드의 내용에서 rnum이라는 변수의 값만 기존의 21에서 133으로 바뀌주기만 하면 나머지 과정들은 동일합니다. 결과로 표출된 그림의 모습은 다음과 같습니다.

이 그림을 보면 육지와 바다가 제대로 구분되어 있다는 것을 충분히 알 수 있습니다. 사실 다른 레코드에 해당되는 데이터들을 표출할 경우에는 육지와 바다 부분이 제대로 구분되어 있는지를 가늠하기가 힘든데, 이렇게 "land-sea mask"에 해당되는 격자 자료를 표출해보면 쉽게 판단이 가능합니다. 참고로 이 그림에서만큼은 컬러바는 별 의미가 없다고 봐도 무방합니다. 어차피 0 아니면 1이란 값만 존재하므로 육지는 붉은색, 바다는 보라색으로 두가지 색상으로만 표출되었기 때문입니다.

아마 기상 분야에서 수치모델 자료를 다뤄보신 분들은 아시겠지만 LDAPS(국지예보모델) 말고도 RDAPS(지역예보모델)라고 하는 자료도 있습니다. RDAPS 데이터 파일 역시 비슷한 요령으로 읽고 처리하는 것이 가능합니다. 다만 내부적으로 레코드 갯수나 번호별 데이터 등은 다르기 때문에, 이러한 부분만 확인해서 데이터 특성에 맞게 조정해주면 지금까지 소개된 예제 코드의 내용을 거의 그대로 활용하여 파일을 읽고 데이터를 표출하는 것이 가능합니다. 지금까지 2회에 걸쳐서 GRIB2 포맷인 LDAPS 데이터를 읽고 그 안의 레코드 및 각종 메타 정보들을 추출하여 이를 바탕으로 지도 및 데이터 중첩 표출하는 과정들을 종합적으로 살펴보았습니다. GRIB/GRIB2 포맷의 데이터를 다뤄야 하는 IDL 유저들께서는 이 내용을 한번 잘 참조해두시면 좋을 것 같습니다.