HDF-EOS 형식의 파일 읽기

HDF-EOS(Hierarchical Data Format - Earth Observing System)는 Terra, Aqua, Aura 등과 같은 EOS(Earth Observing System) 위성들의 자료들의 표준 포맷입니다. 사실 HDF-EOS 포맷도 근본적으로는 HDF의 범주에 속합니다. 다만 HDF-EOS의 경우는 자유로운 HDF 포맷의 일부 형식을 강제하는 대신 EOS의 특성에 맞게 지구과학적 자료를 저장하고 접근하는 데 편리하도록 특화되어 설계되어 있습니다. HDF-EOS 포맷도 근본적으로는 HDF 포맷이기 때문에 IDL에서도 이러한 파일을 다루는데 있어서 HDF 관련 루틴들을 사용하는 것도 가능은 합니다. 하지만 IDL에서는 HDF-EOS 포맷의 데이터 처리를 위한 전용 라이브러리가 별도로 지원됩니다. 따라서 HDF-EOS 포맷의 파일들을 IDL에서 취급하기 위해서는 HDF-EOS 관련 루틴들을 사용하는 것이 좋습니다. 이러한 기능에 대한 전체적인 내용은 IDL 도움말에서 Contents -> Input/Output -> HDF-EOS Routines 섹션에서 볼 수 있으며, 여기서는 개별 루틴별로 자세히 소개되어 있습니다.

 

실제로 이 내용을 보면 HDF-EOS 관련 루틴들은 그 이름이 모두 'EOS_'로 시작합니다. 그리고 세부적으로는 EOS_GD, EOS_PT, EOS_SW 등의 세가지 범주로 나눠집니다. 이것은 HDF-EOS 파일 내에 수록될 수 있는 세가지 종류의 데이터 타입(Grid, Point, Swath) 각각에 대한 루틴들로 나눠져있는 것입니다. 즉 HDF-EOS 데이터는 GRID, SWATH, POINT의 세가지 형식들로 나눠지도록 정의되어 있고, 각 종류의 데이터에 접근할 수 있도록 해주는 기능 API가 형식별로 따로 존재합니다. 세가지 형식에 대한 개요는 다음과 같습니다.

 

1. GRID는 2차원 배열 형태로, 잘 정의된 지도투영법에 의해 저장된 자료입니다. 이름이 'EOS_GD'로 시작하는 루틴들로 처리합니다.
2. POINT는 비정형으로 관측데이터와 관측점의 정보가 함께 저장되는 형태입니다.  이름이 'EOS_PT'로 시작하는 루틴들로 처리합니다.
3. SWATH는 위성이 쓸고 지나간 정확한 시간 정보가 중요할 때 사용됩니다. 실제 관측의 스캔 작업과 관련된 정보가 유지되는 데이터 형태입니다.  이름이 'EOS_SW'로 시작하는 루틴들로 처리합니다.

 

사실 세가지 종류들 모두에 대하여 자세하게 다루기는 어렵기 때문에, 오늘은 일단 GRID 형식의 데이터가 수록된 HDF-EOS 파일을 대상으로 예제를 소개해보고자 합니다. 그래서 MODIS 관측 자료를 Level 3로 가공하여 배포하는 GRID 데이터가 수록된 파일을 예제용으로 사용하겠습니다. 제가 이 파일을 내려받은 웹사이트는 다음과 같습니다.

 

lpdaac.usgs.gov/products/mod13q1v006/

 

그리고 파일은 제가 아래에 링크로 첨부해놓겠습니다. 링크를 눌러서 받으시면 됩니다.

 

파일 받기

 

이 파일의 이름을 보면 'MOD13Q1'로 시작하는데, 이는 MODIS 위성 관측 자료로서 16일 간격으로 제공되는 250m 공간해상도의 데이터가 저장되어 있으며, 내부적으로 NDVI와 EVI 등의 요소들을 포함하고 있습니다. 수록된 데이터의 배열 크기는 4800x4800입니다. 우선 이 파일에 수록된 내용을 미리 전반적으로 조회하는 것이 필요할 수도 있는데, 이를 위하여 HDF_BROWSER라는 함수를 사용하면 편리합니다. 다음과 같이 파일 이름을 변수에 저장하고 HDF_BROWSER 함수를 바로 사용해봅시다.

 

file = 'data/MOD13Q1.A2009049.h28v04.006.2015187063326.hdf'

ok = HDF_BROWSER(file)

 

그러면 다음 그림과 같이 팝업 GUI가 하나 뜨게 되며, 사용자는 이 안에서 간단한 마우스 조작에 의하여 파일 내 정보에 대한 전반적인 개요를 살펴볼 수 있습니다.

 

 

여기서 상단을 보면 Swath, Point, Grid 중 하나를 선택할 수 있는 풀다운 메뉴가 보입니다. 이 풀다운 메뉴에서 HDF-EOS Summary 항목이 선택된 상태에서는 이 파일에 1개의 Grid 데이터셋이 저장되어있다는 것을 확인할 수 있습니다. 그리고 Grid 항목을 선택하면 해당 데이터셋 내에 존재하는 데이터 필드들의 목록을 보여 줍니다. 실제로 확인해보면 250m 16days NDVI, 250m 16days EVI, 250m 16days VI Quality 등의 필드들이 보입니다. 여기서 하나의 필드를 선택하고, 바로 우측 상단에 있는 Preview 버튼을 누르면 대략적인 데이터의 형태를 투박하게나마 시각적으로도 확인할 수도 있습니다. 인터페이스의 우측 하단을 보면 Preview의 형태를 image, surface, contour 중에서 선택할 수 있게 되어 있습니다. 그리고 자료가 공간적으로 큰 편이기 때문에 Fit to window 옵션을 체크해두는 것이 더 나을 것입니다. 어쨌든 이와 같이 예제 파일 내에 한 개의 Grid 데이터셋이 존재하며 이 안에 여러 개의 데이터 필드들이 있다는 것도 확인하고 각 데이터 필드의 이름까지 확인하는 정도면 HDF_BROWSER 함수를 사용한 목적은 충분히 달성한 셈입니다.

 

하지만 본격적인 처리를 위해서는 위와 같이 팝업 GUI를 통하여 인터액티브하게 확인하기보다는 작업흐름을 프로그램으로 처리하는 것이 궁극적으로는 필요하기 때문에, HDF_BROWSER 함수의 사용은 그냥 사전 참조용이었다고 보는 것이 좋습니다. 이제 실제 프로그래밍을 위해서는 파일 내에 수록된 데이터에 대하여 조회하고 그 내용을 가져오는 것이 필요한데 이를 위하여 다음과 같이 EOS_QUERY 함수를 사용하게 됩니다.

 

result = EOS_QUERY(file, info)

HELP, info

 

여기서 좌변에 result라는 변수가 있는데, EOS_ 계열의 함수들 대부분이 이와 같이 좌변으로 명령 실행의 성공여부를 1 또는 0으로 되돌려줍니다. 사실 정말로 필요한 정보는 위에서 info라는 변수로 받게 되는데, 여기서 info는 실제로는 구조체(Structure)의 형태가 됩니다. HELP 명령에 의하여 출력된 내용을 보면 다음과 같습니다.

 

** Structure <65537788>, 6 tags, length=64, data length=60, refs=1:

   GRID_NAMES      STRING    'MODIS_Grid_16DAY_250m_500m_VI'

   NUM_GRIDS       LONG                 1

   NUM_POINTS      LONG                 0

   NUM_SWATHS      LONG                 0

   POINT_NAMES     STRING    ''

   SWATH_NAMES     STRING    ''

 

이 내용을 보면 GRID 데이터셋이 한 개 있고 그 이름이 'MODIS_Grid_16DAY_250m_500m_VI'라는 것을 확인할 수 있습니다. 이는 우리가 앞서 HDF_BROWSER를 통하여 확인했던 내용과 일치합니다. 일반적으로는 HDF-EOS 파일 내에 저장된 GRID 데이터셋의 갯수가 하나일 수도 있지만 다수가 될 수도 있습니다. 하나이든 다수이든간에 GRID 데이터셋의 이름들은 info 구조체의 GRID_NAMES라는 필드에 저장되는데 다음과 같이 그 필드에 대하여 구체적인 정보를 확인해볼 수도 있습니다.

 

HELP, info.GRID_NAMES

 

이 때 출력된 내용은 다음과 같습니다.

 

<Expression>    STRING    = 'MODIS_Grid_16DAY_250m_500m_VI'

 

만약 다수의 데이터셋들이 존재할 경우에는 info.GRID_NAMES는 문자값 배열이 될 것입니다. 물론 지금의 예제 파일에서는 데이터셋이 단 하나이기 때문에 위와 같이 그냥 단일 문자값으로 확인되는 상황이라고 보면 됩니다.

 

그러면 이제부터는 파일로부터 데이터를 추출하는 작업을 본격적으로 진행해봅시다. 가장 먼저 이 파일을 Open해야 합니다. 이것은 파일과의 교신을 위한 연결통로를 확보하는 것이라고 보면 됩니다. 이를 위해서는 다음과 같이 EOS_GD_OPEN 함수를 사용하면 됩니다.

 

hdf_fid = EOS_GD_OPEN(file, /READ)

 

여기서는 파일을 읽는 것이므로 /READ 키워드를 사용하였습니다. 이렇게 하면 좌변의 hdf_fid 변수를 통하여 연결통로에 대한 고유번호를 돌려받게 되며 이 번호를 후속 작업에서 계속 사용하게 됩니다. 그 다음 과정에서는 GRID 데이터셋에 대한 고유 번호를 가져오기 위하여 다음과 같이 EOS_GD_ATTACH 함수를 사용합니다.

 

grid_id = EOS_GD_ATTACH(hdf_fid, info.grid_names[0])

 

이렇게 하면 좌변의 grid_id 변수를 통하여 첫번째(그리고 여기서는 유일한) 데이터셋에 대한 고유번호를 돌려받게 되며, 이 번호를 후속 작업에서도 계속 사용하게 됩니다. 그 다음에는 이 데이터셋에 관한 몇가지 부가정보들을 추출하기 위하여 다음과 같이 EOS_GD_GRIDINFO 함수를 사용해봅시다.

 

result = EOS_GD_GRIDINFO(grid_id, xdim,ydim, upleft,lowright)

PRINT, xdim, ydim, upleft, lowright

 

여기서는 지정된 GRID 데이터셋에 관한 4종의 정보들을 개별 변수로 받아온 다음 출력하게 되는데, 실제로 출력된 내용은 다음과 같습니다.

 

        4800        4800

       11119505.       5559752.6

       12231456.       4447802.1

 

출력된 내용을 통하여 GRID 데이터셋의 배열 크기(4800x4800)와, 배열의 좌측상단 좌표, 우측하단 좌표를 확인할 수 있습니다. 그 다음 과정은 GRID 데이터셋 내에 구체적으로 어떤 데이터들이 수록되어 있는가를 확인하는 것인데요. 하나의 데이터셋이 거느리고 있는 실제 데이터 항목들 각각을 필드(Field)라고 합니다. 이러한 정보를 확인하기 위하여 다음과 같이 EOS_GD_INQFIELDS 함수를 사용합니다.

 

result = EOS_GD_INQFIELDS(grid_id, fieldlist, rank, numbertype)

HELP, fieldlist

 

여기서는 fieldlist라는 이름으로 받게 되는 결과를 주목할 필요가 있는데요. HELP에 의하여 출력된 내용은 다음과 같습니다.

 

FIELDLIST       STRING    = '250m 16 days NDVI,250m 16 days EVI,250m 16 days VI Quality,250m 16 days red refle'...

 

실제로 fieldlist 자체는 그냥 단일 문자값을 담은 변수입니다. 그런데 이 변수가 담고 있는 문자값의 내용을 자세히 보면, 데이터셋 내에 수록된 모든 필드들의 이름이 코마(,)로 연결된 거대한 단일 문자값임을 알 수 있습니다. 그래서 통상적으로는 이 거대한 문자값에 대하여 다음과 같이 STRSPLIT 함수를 추가적으로 사용함으로써, 필드별 이름 문자값들로 구성된 문자형 배열이 fieldlist가 되도록 바꿔버리는 것이 좋습니다.

 

fieldlist = STRSPLIT(fieldlist, ',', /EXTRACT)

HELP, fieldlist

FOR j = 0, N_ELEMENTS(fieldlist)-1 DO $

  PRINT, j, fieldlist[j], FORMAT='(I2, 3X, A)'

 

여기서 HELP에 의하여 출력된 내용은 다음과 같습니다.

 

FIELDLIST       STRING    = Array[12]

 

즉 원래의 거대한 문자값을 코마(,)로 분할하여 12개의 문자값들로 구성된 배열이 되었음을 알 수 있습니다. 그리고 반복형 구문의 PRINT에 의하여 출력된 내용은 다음과 같습니다.

 

 0   250m 16 days NDVI

 1   250m 16 days EVI

 2   250m 16 days VI Quality

 3   250m 16 days red reflectance

 4   250m 16 days NIR reflectance

 5   250m 16 days blue reflectance

 6   250m 16 days MIR reflectance

 7   250m 16 days view zenith angle

 8   250m 16 days sun zenith angle

 9   250m 16 days relative azimuth angle

10   250m 16 days composite day of the year

11   250m 16 days pixel reliability

 

이와 같은 방식으로 12개의 모든 필드들의 이름 및 인덱스 번호를 확인할 수도 있습니다. 이제 필드로부터 실제 데이터를 추출해봅시다. 여기서는 0번 NDVI 데이터를 추출하겠습니다. 이렇게 데이터셋으로부터 특정 필드의 데이터를 추출하기 위해서는 다음과 같이 EOS_GD_READFIELD 함수를 사용하면 됩니다.

 

result = EOS_GD_READFIELD(grid_id, fieldlist[0], ndvi)

HELP, ndvi

 

여기서는 맨 첫번째 필드를 ndvi라는 이름으로 추출하였는데, HELP에 의하여 출력된 내용을 보면 다음과 같습니다.

 

NDVI            INT       = Array[4800, 4800]

 

이와 같이 4800x4800의 구조를 갖는 정수형 배열로 추출되었습니다. 참고로 이렇게 정수형으로 수록되어 있는 이유는 용량의 절약을 위해서입니다. 데이터의 속성 자체로 보면 실수형이어야 맞습니다. 하지만 실수형 값은 개당 8바이트(일반 정밀도일 경우)인 반면 정수형 값은 개당 4바이트이기 때문에, 용량의 절약 목적으로 원래의 실수값에 10000을 곱하여 정수형으로 변환되어 저장된 것이라고 보면 됩니다. 따라서 이렇게 추출된 데이터를 실제로 사용하기 위해서는 0.0001을 곱하여 보정하는 것이 필요합니다. 다만 이러한 보정을 수행하기 전에 다음과 같이 배열 ndvi의 최소값 및 최대값을 먼저 출력해봅시다.

 

PRINT, MIN(ndvi), MAX(ndvi)

 

출력된 값들을 보면 다음과 같습니다.

 

   -3000    9986

 

즉 ndvi 배열의 경우를 보면 최소값이 -3000인데요. 실제로는 NDVI가 계산될 수 없는 부분에 대하여 이렇게 -3000이라는 값이 인위적으로 삽입되어 있습니다. 즉 -3000은 아무 의미가 없는 값입니다. 따라서 여기서는 이러한 값들을 걸러내어서 NaN값으로 대체하는 것이 필요하다고 봅니다. 그래서 이러한 처리를 위하여 다음과 같이 WHERE 함수를 사용하도록 하겠습니다.

 

bad = WHERE(ndvi EQ -3000)

ndvi = ndvi*0.0001

ndvi[bad] = !values.f_nan

HELP, ndvi

PRINT, MIN(ndvi), MAX(ndvi)

 

이와 같은 처리를 거친 후에 출력된 ndvi에 관한 정보들을 보면 다음과 같습니다.

 

NDVI            FLOAT     = Array[4800, 4800]

    -0.200000     0.998600

 

이와 같이 ndvi는 이제 정상 범위의 실수형 값들로 구성된 배열이 되었습니다. 그러면 이 정도의 처리를 거친 상태에서 이제 표출을 해봅시다. 이미지의 형태로 표출되도록 하기 위하여 다음과 같이 IMAGE, COLORBAR 등의 함수들을 사용해보았습니다. 과정은 다음과 같습니다.

 

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

i = IMAGE(ndvi, RGB_TABLE=34, AXIS_STYLE=2, MARGIN=0.1, /CURRENT)

cb = COLORBAR(TARGET=i, TITLE='NDVI', FONT_SIZE=11, $

  POSITION=[0.25, 0.95, 0.75, 0.98])

 

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

 

 

추가 작업을 조금만 더 해보자면, HDF-EOS 관련 기능들 중에서는 경위도 값을 받아 이에 해당하는 픽셀의 위치를 알려주는 기능도 있습니다. 이를 위하여 다음과 같이 EOS_GD_GETPIXELS 함수를 사용하면 됩니다.

 

lons = [142.17, 140.9]

lats = [43.5, 40.8]

status = EOS_GD_GETPIXELS(grid_id, N_ELEMENTS(lons), $

  lons, lats, rows, cols)

PRINT, 'ROWS :', rows

PRINT, 'COLS :', cols

 

여기서는 경위도 좌표가 (142.17, 43.5)인 지점과 (140.9, 40.8)인 지점 두 군데를 가정하여 해당 위치에 대응되는 픽셀 좌표를 산출하도록 하였습니다. 출력된 결과는 다음과 같습니다.

 

ROWS :        3119        4416

COLS :        1500        3197

 

이 결과에 의하면 경위도 좌표 (142.17, 43.5)가 픽셀 좌표 (1500, 3119)으로 대응되고, 경위도 좌표 (140.9, 40.8)이 픽셀 좌표 (3197, 4416)으로 환산된다는 것을 확인할 수 있습니다. 더 나아가서 이 좌표들에 대한 NDVI 값을 가져오는 것도 가능합니다. 방금 얻은 픽셀 좌표들을 앞서 우리가 이미 처리했던 ndvi 배열에 투입하여 해당 값들을 인덱싱하면 됩니다. 그 방법은 다음과 같습니다.

 

PRINT, ndvi[cols, rows]

 

출력된 결과는 다음과 같습니다.

 

     0.133600     0.185200

 

사실 이와 같이 특정 필드 데이터에 대하여 픽셀 좌표를 이용하여 해당 값을 가져오는데 작업을 위해서는, 다음과 같이 EOS_GD_GET_PIXELVALUES 함수를 사용하는 방법도 있습니다.

 

status = EOS_GD_GETPIXVALUES(grid_id, $

  N_ELEMENTS(rows), cols, rows, fieldlist[0], getpix)

PRINT, getpix

 

이 때 출력된 결과를 보면 다음과 같습니다.

 

    1336    1852

 

다만 여기서 유의할 것은, 이 값들은 필드로부터 그대로 가져온 것이기 때문에 0.0001을 곱해야 실제 관측값이 된다는 것입니다. 따라서 실제 관측값은 0.1336, 0.1852가 된다고 보면 됩니다. 어쨌든 이러한 방법들이 있다는 정도로만 이해하시면 될 것 같습니다.

 

이제 모든 작업들은 다 끝났습니다. 따라서 프로그램의 내용이 여기서 끝나도 큰 문제는 없습니다. 다만 우리가 HDF-EOS 파일을 읽는 작업을 위하여 앞서 EOS_GD_OPEN, EOS_GD_ATTACH 함수들을 사용하여 파일 및 데이터셋에 대한 연결통로를 오픈했었는데요. 이 통로들을 다음과 같이 EOS_GD_DETACH 및 EOS_GD_CLOSE 함수들을 사용하여 직접 닫아주는 것이 좋습니다.

 

ok = EOS_GD_DETACH(grid_id)

ok = EOS_GD_CLOSE(hdf_fid)

 

물론 이렇게 '닫는' 명령을 실행하지 않고 프로그램이 종료되도록 코딩해도 큰 문제는 없습니다. 어차피 프로그램이 종료되면, 프로그램 내에서 열어놨던 이런 연결통로들은 자동적으로 닫히기 때문입니다. 하지만 전체적인 작업이 명시적으로 마무리되도록 명문화하는 것이 아무래도 더 좋을 수도 있기 때문에, 이와 같은 방법도 알아두시는 것이 좋습니다.

 

내용이 좀 길긴 했지만 오늘은 HDF-EOS 파일을 읽어서 GRID 형식의 데이터를 추출하는 과정에 대한 예제를 소개해보았습니다. 오늘 소개한 내용은 예전에 (주)에스이랩에서 작성했던 HDF-EOS 읽기에 관한 다른 기술문서의 내용을 토대로 하여 약간의 수정 및 업데이트를 한 것입니다. 원래 문서의 링크는 아래와 같습니다.

 

idl-envi.co.kr/callcenter/tips/file/67_1406877911.pdf

 

그리고 앞서 이미 언급했듯이 HDF-EOS 파일 내에는 GRID 뿐만 아니라 POINT 또는 SWATH 형식의 데이터가 수록될 수도 있습니다. 오늘은 아무래도 GRID의 경우만 소개할 수 밖에 없었는데요. 일단 SWATH 형식의 데이터를 추출하는 과정에 대한 예제는 아래 링크에서 보실 수 있습니다. 따라서 이 내용도 아울러 참조하시면 좋을 것 같습니다.

 

idl.selab.re.kr/?p=1862