프로바브이(PROBA-V) 위성 시계열을 이용한 아고산 식생 연간 변화탐지

REMOTE SENSING

프로바브이(PROBA-V) 위성 시계열을 이용한 아고산 식생 연간 변화탐지

유병혁 2017. 1. 13. 00:22

PyQGIS scripting과 D3.js를 이용하여 PROBA-V 위성 시계열을 처리하는 과정을 정리해 보도록 하겠습니다.

프로바브이(PROBA-V) 위성은 유럽우주기구(ESA)에서 2013년 5월 6일에 발사하였으며, SPOT-VGT(1998년

3월 - 2014년 5월)과 Sentinel-3호 영상과 함께 지구 궤도를 돌며 식생을 측정하는 임무를 수행하고 있습니다.

아래 그림처럼 프로바브이 위성의 체적은 불과 0.05 ㎥(50리터)이며, 무게는 33kg입니다.

그러나 이 위성의 관측 폭(Field Of View)은 102˚, 주사 폭(swath width)은 2295 km 로써,

일일 지구영역의 90%를 촬영할 수 있으며 2일 단위로 전 지구 범위 데이터를 취득할 수 있습니다.

사진 출처: http://www.esa.int/Our_Activities/Observing_the_Earth/Proba-V/About_Proba-V

사진 출처: http://www.esa.int/For_Media/Photos/Highlights/Proba-V

이 위성영상을 다운로드 받는 방법은 아래 글을 통해 정리하였습니다.

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

저는 아래기사와 관련하여 지리산국립공원 아고산 일원의 위성 시계열을 분석해 보려고 합니다. 이 글에서는

2014년 1년간 10일 단위로 합성된 수관상부(TOC, Top of Canopy) 정규식생지수(NDVI, 해상도 300m) C1를 이용합니다.

구상나무의 비명…수령 70년 거목 생선 가시처럼 말라죽어

기사원문: http://news.joins.com/article/20510879

기사에 따르면 2014년은 이상기후가 발생했던 시기로 봄철 평균 기온이 평년보다 1.4도 높았었고,

그해 5월 기온은 1973년 이후 가장 높았다고 합니다. 그 달 강수량은 평년의 51.6%에 그쳤습니다.

일단 아래와 같이 2014년 위성 시계열을 다운로드 받고,

아래와 같이 국립공원 구역(Boundary zone)와 해발 1,200m 이상 아고산대 구역(subalpine zone) 폴리곤을 준비합니다.

프로바브이 영상은 8비트 DN 값으로 제공되기 때문에, 아래 표의 계산식으로 실제 정규식생지수(NDVI) 변환을 하도록 합니다.

https://earth.esa.int/web/guest/missions/esa-operational-eo-missions/proba-v/faqs#question_1

PyQGIS scripting을 아래와 같이 적용하면 각각 폴리곤에 구역 통계(zonal statistics)가 반복 계산됩니다.

구역 통계 계산은 아래 글을 참고하시기 바랍니다.

PyQGIS: QGIS 구역 통계(Zonal Statistics)를 반복적으로 실행하기

http://blog.daum.net/geoscience/933

#전체 레이어 갯수 확인
QgsMapLayerRegistry.instance().count()
#벡터 레이어 확인
vlayer = qgis.utils.iface.mapCanvas().layer(0)
print vlayer.name()
#래스터 레이어 확인
rfile = qgis.utils.iface.mapCanvas().layer(1).source()
print rfile
print rfile[92:100]
#구역 통계 반복 계산
from qgis.analysis import QgsZonalStatistics
i = 1
layercount = QgsMapLayerRegistry.instance().count()
while i < layercount:
    rfile = qgis.utils.iface.mapCanvas().layer(i).source()
    # 구역 통계 계산
    zonalstats = qgis.analysis.QgsZonalStatistics(vlayer, rfile,
    rfile[92:100] + "_", 1, QgsZonalStatistics.Mean|QgsZonalStatistics.StDev)
    zonalstats.calculateStatistics(None)
    i += 1