QGIS에서 gdal_translate를 이용하여 위성 데이터 포맷 변환하기

QGIS에서 gdal_translate를 이용하여 위성 데이터 포맷을 변환해 보겠습니다.

gdal_translate | http://www.gdal.org/gdal_translate.html

이 기능은 '래스터 > 변환 > 변환 (포맷 변경)'을 통해 이용하실 수 있는데요,

실행 창은 아래와 같습니다.

gdal_translate의 개요는 아래와 같은데요, 이것을 참고하여 실행 창의   버튼을 클릭하고 내용 편집을 하면 되겠습니다.

실습 위성데이터는 PROBA-V(프로바브이)입니다. 아래 웹사이트에서 받아보실 수 있습니다.

VITO Earth Observation | http://www.vito-eodata.be/PDF/portal/Application.html#Home

PROBA-V 위성영상 다운로드와 WMTS레이어 정보는 아래 2개 링크를 참고하시기 바랍니다.

프로바브이(PROBA-V) 위성영상 다운로드 방법 소개 | http://blog.daum.net/geoscience/1004

QGIS에서 WMTS 레이어로 프로바브이(PROBA-V) 위성영상 열어보기 | http://blog.daum.net/geoscience/1005

여기서는 S10 TOC NDVI - 1km[C1] 데이터를 HDF-5 (Hierarchical Data Format 5) 포맷으로 다운로드 받았습니다.

참고로, S10은 10일주기 합성, TOC는 대기보정이 적용된 Top-Of-Canopy, 1km는 공간해상도, C1은 보정 수준을 뜻합니다.

HDF 포맷은 HDF Group(https://www.hdfgroup.org/)에 의해 관리되고 있는데요,
공식 홈페이지: https://support.hdfgroup.org/HDF5/

QGIS는 HDF5와 HDF4(기존 포맷)을 모두 공식 지원하고 있습니다.

자, 먼저 HDF5 포맷 데이터를 QGIS에서 열어볼까요?!

1. HDF5를 GeoTIFF로 포맷 변환하기:

gdal_translate
-of GTiff
D:/PROBAV_Sample.hdf5 D:/PROBA_V_Sample.tif

현재 QGIS는 HDF5 Geotie point definition과의 호환성 문제로 PROBA-V의 CRS를 지정해도 영상 좌표계로 열리게 됩니다. 

이에 대해 PROBA-V 설명서는 타 GIS 프로그램 사용을 권장하고 있는데요, 이 부분을 gdal_translate로 해결해 보겠습니다.

PROBA-V Product User Manual (61p 참조) | http://proba-v.vgt.vito.be/sites/proba-v.vgt.vito.be/files/Product_User_Manual.pdf

2. 대상 좌표계를 EPSG:4326으로 정의하기:

gdal_translate
-of GTiff
-a_srs EPSG:4326
D:/PROBAV_Sample.hdf5 D:/PROBA_V_Sample.tif

3. 대상 좌표계(EPSG:4326)를 기준으로 레이어 범위 정의하기:

gdal_translate
-of GTiff
-a_srs EPSG:4326
-a_ullr 119.9955357142857082 45.0044642857142847 129.9955357142852392 35.0044642857147608
D:/PROBAV_Sample.hdf5 D:/PROBAV_Sample.tif

위와 같이 정의하면, QGIS의 HDF5 Geotie point definition 호환 문제를 해결할 수 있겠죠?!

이번에는 데이터 유형을 실수로 변경한 후, DN값을 NDVI 값으로 변환해 보도록 하겠습니다.

4. 데이터 유형을 Float32 - 32 bit 실수로 변경하기:

gdal_translate
-ot Float32
-of GTiff
-a_srs EPSG:4326
-a_ullr 119.9955357142857082 45.0044642857142847 129.9955357142852392 35.0044642857147608
D:/PROBAV_Sample.hdf5 D:/PROBAV_Sample.tif

5. DN 값을 NDVI 값으로 변환하기:

gdal_translate
-ot Float32
-of GTiff
-scale 0 250 -0.08 0.92
-a_srs EPSG:4326
-a_ullr 119.9955357142857082 45.0044642857142847 129.9955357142852392 35.0044642857147608
D:/PROBAV_Sample.hdf5 D:/PROBAV_Sample.tif

DN 값을 NDVI 값으로 변환하는 식은 아래와 같습니다(PROBA-V Product User Manual : 53p 참조).

위 식에서 PV는 Physical Values, DN은 Digital Count Number, offset과 scale은 각각 20, 250입니다.

계산해 보면 DN 값이 [0, 250] 범위에서 NDVI 값은 [-0.08, 0.92]범위를 갖습니다.

6. No data 설정하기:

gdal_translate
-ot Float32
-of GTiff
-scale 0 250 -0.08 0.92
-a_srs EPSG:4326
-a_ullr 119.9955357142857082 45.0044642857142847 129.9955357142852392 35.0044642857147608
-a_nodata 0.939999997615814
D:/PROBAV_Sample.hdf5 D:/PROBAV_Sample.tif

PROBA-V DN 값에서 255는 No Data에 해당됩니다. 255의 PV 값인 0.939999997615814를 No Data로 정의합니다(아래 화면은 No Data 예시).

7. QGIS 파이썬 콘솔에서 gdal_translate 반복 실행하기:

import os
srcPath = "D:/GEOINFO/PROBA-V_1km/X30Y03_HDF5"
dstPath = "D:/GEOINFO/PROBA-V_1km/X30Y03_GeoTIFF"
for srcName in os.listdir(srcPath):
    dstName = srcName.replace(".hdf5", ".tif")
    srcDataset = os.path.join(srcPath, srcName)
    dstDataset = os.path.join(dstPath, dstName)
    os.system("gdal_translate
-ot Float32
-of GTiff
-scale 0 250 -0.08 0.92
-a_srs EPSG:4326
-a_ullr 119.9955357142857082 45.0044642857142847 129.9955357142852392 35.0044642857147608
-a_nodata 0.939999997615814 " + srcDataset + " " + dstDataset)

시계열 위성영상 처리에 유용하겠죠?! 아래는 SPOT Vegetation을 처리한 예시입니다.

Southeast Asia 영역을 반복 변환한 후, 위 PROBA-V 레이어 영역으로 잘라냈습니다.

import os
srcPath = "D:/GEOINFO/SPOT_VGT/SE_Asia_NDV_HDF"
dstPath = "D:/GEOINFO/SPOT_VGT/SE_Asia_NDV_GeoTIFF"
for srcName in os.listdir(srcPath):
    dstName = srcName.replace(".HDF", ".tif")
    srcDataset = os.path.join(srcPath, srcName)
    dstDataset = os.path.join(dstPath, dstName)
    os.system("gdal_translate
-ot Float32
-of GTiff
-scale 0 255 -0.1 0.92
-a_srs EPSG:4326
-a_ullr 67.9955357143000043 55.0044642857000028 147.0044640328999890 4.9955358743000033 "
+ srcDataset + " " + dstDataset)

import os
srcPath = "D:/GEOINFO/SPOT_VGT/SE_Asia_NDV_GeoTIFF"
dstPath = "D:/GEOINFO/SPOT_VGT/X30Y03_GeoTIFF"
for srcName in os.listdir(srcPath):
    dstName = srcName.replace("_SE-Asia_", "_X30Y03_")
    srcDataset = os.path.join(srcPath, srcName)
    dstDataset = os.path.join(dstPath, dstName)
    os.system("gdal_translate
-projwin 119.9955357142857082 45.0044642857142847 129.9955357142852392 35.0044642857147608
-projwin_srs EPSG:4326 " + srcDataset + " " + dstDataset)