QGIS: Burn Severity(산불피해 심각도) 분석 (5) - 임상도를 이용한 공간통계 산출

REMOTE SENSING

QGIS: Burn Severity(산불피해 심각도) 분석 (5) - 임상도를 이용한 공간통계 산출

유병혁 2022. 4. 9. 15:07

안녕하세요? 이번 글은 Burn Severity(산불피해 심각도) 분석 실습 시리즈 다섯번째 글로 임상도를 이용한 공간통계 산출과정을 정리해 보겠습니다. 이전 글은 아래 링크를 산출하시기 바랍니다.

QGIS: Burn Severity(산불피해 심각도) 분석 (4) - dNBR-임상도 교차 영역 산출

안녕하세요? 이번 글은 Burn Severity(산불피해 심각도) 분석 실습 시리즈 네번째 글로 dNBR과 임상도(1:5000) 교차 영역을 산출해 보겠습니다. 이전 글은 아래 링크를 산출하시기 바랍니다. QGIS: Burn Seve

blog.daum.net

공간통계는 Python 프로그래밍으로 처리하려고 하는데요, 이를 위해 스프레드시트 데이터를 먼저 추출해 보겠습니다. 앞서 실습을 통해 산출한 dNBR.gpkg에서 'dNBR_FOREST_STOCKMAP_5000 레이어의 속성 테이블을 확인해 봅니다. 실습에 필요한 속성 값들이 정리되어 있습니다. 해당 속성 테이블에 면적 값 필드를 추가해 보겠습니다.

편집 모드 전환 후,

필드 계산기 열기를 실행합니다.

아래와 같이 $area / 10000 표현식을 통해 헥타르 면적 값을 포함하는 'HA' 필드를 추가해 줍니다.

자, 추가되었죠?!

이 상태에서 레이어 이름을 오른쪽 마우스 클릭한 후 'Export > 객체를 다른 이름으로 저장'을 실행합니다.

포맷은 'MS Office Open XML spreadsheet [XLSX]'를 선택하고 파일을 저장해 줍니다.

자, 이제 Python 프로그래밍으로 이어가 볼까요?! 사용할 모듈은 아래와 같습니다.

# 모듈 호출
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style='darkgrid', font='KoPubDotum', font_scale=1.2)

해당 *.xlsx 파일을 읽어와 봅니다.

df = pd.read_excel('D:\\GEODATA\\dNBR_FOREST_STOCKMAP_5000.xlsx')
df.head()

임상도는 우리나라 국토의 산림이 어떻게 분포하고 있는가를 보여주는 대표적인 산림지도로 임상, 영급, 경급, 임종, 수종 등 다양한 속성정보를 포함(아래 링크 글 참고)하고 있습니다. 우리는 앞서 dNBR을 기반으로 산불피해 지역을 산출했고 이 영역에 해당하는 임상도만을 별도 추출했습니다. 해당 임상도에 대한 공간통계를 통해 산불피해 지역에 특성을 확인해볼 수 있겠습니다.

임상도(1:5000) QGIS 스타일 파일(*.qml) 공유

안녕하세요? 이번 글은 임상도(1:5000) QGIS 스타일 파일(*.qml)을 공유 드려 봅니다. 임상도를 사용하시는 분들은 공유해드린 파일을 통해 스타일링을 적용하실 수 있습니다. 임상도는 우리나라 국

blog.daum.net

1. 임상 공간통계

임상코드를 기준으로 산불피해 통계를 확인해 보겠습니다. 아래 그림과 같이 임상유형별 헥타르 면적(HA)과 각 유형별 면적비율(PCT)을 확인할 수 있습니다.

# 임상
df1 = df.groupby('FRTP_CD').sum()[['HA']]
df1['PCT'] = (df1['HA'] / df1['HA'].sum()) * 100
df1

차트화 해보니 뚜렷하게 임상별 공간통계가 확인됩니다. 침엽수림이 75.2%, 활엽수림이 20.2% 해당하는 것을 확인할 수 있습니다. 최근 산불 관련 기사에 따르면 국내 산림 활엽수림은 27%에 불과하다고 합니다. 국립산림과학원에 따르면 침엽수는 화염 유지 시간이 57.3초, 활엽수는 23.0초로 침엽수림에서 산불 확산 속도가 빠르기 때문에 활엽수로의 수종 갱신이 앞으로 필수조건이라고 합니다.

# 임상
fig, ax = plt.subplots(figsize=(16, 4))
ax = sns.barplot(x = "FRTP_CD", y = "HA", estimator=np.sum, palette='magma', ax=ax, ci=None, data=df)
ax.set(xlabel='임상', ylabel='면적 (ha)')
plt.xticks(range(5), ['무립목지/비산림', '침엽수림', '활엽수림', '혼효림', '죽림'])
for p in ax.patches:
    total = df1['HA'].sum()
    ax.annotate('{:.1f}%'.format(100*p.get_height()/total),
                   (p.get_x() + p.get_width() / 2., p.get_height()),
                   ha='center', va='center',
                   size=12,
                   xytext=(0, 5),
                   textcoords='offset points')
plt.show()