안녕하세요? 이번 글에서는 QGIS 구역 통계(Zonal Statistics)를 반복적으로 실행하는 방법을 정리해 보도록 하겠습니다.
일단, 아래와 같이 QGIS Oceancolor Data Downloader 플러그인을 통해 Aqua 위성에 탑재된 MODIS(Moderate
Resolution Imaging Spectroradiometer) 센서의 해수면온도(SST; Sea Surface Temperature) 데이터를 내려받았습니다.
2002년 7월부터 2016년 2월까지 총 164장의 월 단위 4km 해상도 데이터인데요,
저는 우리나라 서해에 위치한 태안해안 국립공원, 남해에 위치한 다도해해상과 한려해상 국립공원 경계 범위를 구역 통계의 기준으로 설정하였습니다.
표층 수온의 변화(특히 서해)로 인해 해상형 국립공원 어종의 변화도 뚜렷해지고 있는데요,
구역 통계를 통해 해상형 국립공원 별로 지난 14년간 어떠한 추세가 있는지 살펴보고자 합니다.
자, 상단 메뉴에서 '플러그인 > 파이썬 콘솔' 메뉴를 클릭합니다.
아래와 같이 '파이썬 콘솔' 창이 추가되었는데요,
일단, 현재 QGIS에 추가된 레이어 갯수를 확인해 보겠습니다. 총 165장이라고 나오네요?!
첫 번째 레이어는 해상형 국립공원 경계 범위인데요, 아래 코드에서 layer(0)이 최상위에 위치한 레이어를 지시합니다.
단, 해당 레이어가 최상위에 있더라도 현재 꺼져있는 상태(예: 모든 레이어 숨기기)라면 위 코드는 적용되지 않습니다.
이번에는 래스터 레이어 명을 한 번 읽어봤습니다.
구역 통계를 컬럼으로 저장할 때, 다른 컬럼과의 구분을 위해 접두사가 필요한데요, 아래 코드는
전체 파일명에서 접두사로 쓸 부분만 선별('2016032는 2016년 32번째 날의 의미임')하기 위한 과정입니다.
자, 원래의 구역 통계는 어떻게 계산되는지 한 번 살펴볼까요?!
아래와 같이 상단 메뉴에서 '래스터 > 구역 통계 > 구역 통계'를 클릭하면,
아래와 같이 구역 통계가 필행됩니다. 기능은 아주 간단하죠?! 저는 여기서는 3개 공원 범위 내에서 평균(Mean),
중앙값(Median), 표준편차(Standard Deviation), 최소값(Min), 최대값(Max)을 반복 계산해 보도록 하겠습니다.
다시 파이썬 콘솔로 돌아와 아래와 같이 코드를 입력해 볼 건데요, 아래 코드는 위 구역 통계 기능을 파이썬에서
그대로 구현한 것입니다. 아래와 같이 값이 '0'이라고 표시되면 해당 기능이 잘 완료되었다고 보시면 되겠습니다.
테이블을 열어본 결과입니다. 간단하죠?! 아래는 2016년 2월 자료인데요, 평균값에서 태안해안 표층수온은
4.85℃에 불과한 반면, 다도해하상과 한려해상의 표층수온은 각각 8.31℃, 9.62℃인 것을 확인하실 수 있습니다.
이제 아래와 같이 레이어 갯수만큼 구역 통계를 반복적으로 실행합니다.
끝으로 본 글에서 사용한 파이썬 소스 코드는 아래와 같습니다.
#전체 레이어 갯수 확인 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[28:35] #구역 통계 계산 from qgis.analysis import QgsZonalStatistics zonalstats = qgis.analysis.QgsZonalStatistics(vlayer, rfile, rfile[28:35] + "_", 1, QgsZonalStatistics.Mean|QgsZonalStatistics.Median |QgsZonalStatistics.StDev|QgsZonalStatistics.Min|QgsZonalStatistics.Max) zonalstats.calculateStatistics(None) #구역 통계 반복 계산 i = 1 layercount = QgsMapLayerRegistry.instance().count() while i < layercount: rfile = qgis.utils.iface.mapCanvas().layer(i).source() # 구역 통계 계산 zonalstats = qgis.analysis.QgsZonalStatistics(vlayer, rfile, rfile[28:35] + "_", 1, QgsZonalStatistics.Mean|QgsZonalStatistics.Median |QgsZonalStatistics.StDev|QgsZonalStatistics.Min|QgsZonalStatistics.Max) zonalstats.calculateStatistics(None) i += 1
다음 글에서는 이렇게 추출한 월별 해수면 온도 데이터를 R 프로그래밍을 통해 다뤄보도록 하겠습니다.