QGIS: Python 종 분포 모델링(SDM) (3) - 종 적합성(species suitability) 매핑

안녕하세요? QGIS와 Python을 이용한 종 분포 모델링(SDM: Species Distribution Modeling)을 시리즈 글로 정리해 보겠습니다. 이번에는 세 번째 실습으로 종 적합성(species suitability) 매핑을 주된 내용으로 정리해 보겠습니다.

*이 글은 국립공원 비대면 학습동아리 '파크랩(ParkLab)'의 학습용 자료입니다. 스터디에 참여하는 조경학, 컴퓨터공학 등 다양한 전공자분들의 참여를 통해 콘텐츠는 수시로 수정, 보완될 수 있습니다. 유튜브 영상과 프레젠테이션 자료는 아래 링크를 참고하시면 됩니다.

QGIS: Python 종 분포 모델링(SDM)

이전 실습 글은 다음 링크를 참조하시면 됩니다.

QGIS: Python 종 분포 모델링(SDM) (1) - Python 설치 및 환경 설정

QGIS: Python 종 분포 모델링(SDM) (2) - 실습 데이터 전처리

조슈아 트리 연구지역으로 자른 생물 기후 특징 ASCII Grid 포맷을 INPUT 폴더에 복사합니다. 아래 코드를 풀어보면, glob.glob은 제시한 조건에 맞는 파일 이름을 리스트형으로 반환하고, sorted는 새로 정렬된 리스트를 만들어 줍니다. 그리고 for in 반복문을 통해 각각의 파일은 shutil.copy로 INPUT 폴더에 복사됩니다.

# 조슈아 트리 연구지역으로 자른 생물 기후 특징
for f in sorted(glob.glob('DATA/BIOCLIM/bclim*.asc')):
    shutil.copy(f,'INPUT/')

해당 리스트를 raster_features(래스터 특징)라는 변수로 정의한 후 출력한 결과입니다. 앞서 우리는 train_vec에다가 종 발생 여부와 19가지 환경변수 값을 담았었습니다. 한 가지 유의할 부분은, 래스터 특징과 훈련 데이터의 레이어 순서는 서로 일치해야 한다는 것입니다. 예를 들면, raster_features는 bclim1, bclim2... 순이고 train_vec는 bclim1, bclim10... 순으로 불일치하면 안됩니다.

# 래스터 특징
raster_features = sorted(glob.glob(
    'INPUT/bclim*.asc'))
print(raster_features)

래스터 특징 수를 확인해봅니다. 기온과 강수량 통계로 된 19개 생물 기후 특징이니 19가 나와야겠죠?!

# 래스터 특징 수 확인
print(len(raster_features), '개 래스터 특징')

자, 이번에는 필요한 pyimute 모듈을 불러와 봅니다. scikit-learn과 rasterio를 이용한 지리공간 예측용 Python 모듈입니다. 각각의 공식 웹페이지는 다음을 참고하시기 바랍니다.

• pyimpute(지리공간 예측용 Python 모듈): https://github.com/perrygeo/pyimpute
• scikit-learn(Python용 자유 소프트웨어 ML 라이브러리): https://scikit-learn.org/stable/
• rasterio(지리공간 래스터 데이터 입출력 모듈): https://rasterio.readthedocs.io/en/latest/ 

# pyimpute 모듈
from pyimpute import load_targets

아래는 pyimpute 홈페이지에서 제공하는 이미지로 '머신러닝 감독 분류(Machine Learning Supervised Classification)'를 표현한 것입니다.  훈련 데이터(training data)로 알려진 관측치(observation)는, 반응변수(response variables)와 설명변수(explanatory variables)로 구성되어 있습니다.

• 반응 변수(response variables, Y): 우리가 예측(predict)하려고 하는 것
• 설명 변수(explanatory variables, Xs): 반응의 공간적 패턴(spatial patterns)을 설명하는 변수

그리고 대상 데이터(target data)는 (래스터 데이터셋으로 대표되는) 설명 변수만으로 구성되며, 사용 가능한 반응 변수는 없습니다. 우리의 목표는 그 반응 변수의 래스터 표면을 예측(predict)하는 것으로, 해당 값은 이산(분류) 또는 연속(회귀)일 수 있습니다.

이미지 출처: https://github.com/perrygeo/pyimpute

자, 그럼 train_vec 변수에서 훈련 데이터를 로드해 보겠습니다. 실습에서 우리는 조슈아 트리 7,200건을 다루고 있습니다. 따라서 설명 변수 Xs의 유형은 7200*19 행렬, 반응 변수 Y의 유형은 7200*1 열벡터가 될 것입니다. 아래 소스에서 pyimpute 모듈의 load_targets는 래스터 특징 변수를 scikit-learn 데이터 구조로 로드하는 역할을 수행해 줍니다.

train_xs, train_y = train_vec.iloc[:,1:20].values, train_vec.iloc[:,0].values # 훈련 데이터 로드
target_xs, raster_info = load_targets(raster_features) # scikit-learn용 데이터 구조로 래스터 특징 로드
train_xs.shape, train_y.shape # 관측치 크기와 일치하는지 shape 확인

*Python 언어 중 iloc 함수가 익숙치 않으신 분들을 위해 예제로 작성해 봤습니다. iloc는 pandas에서 제공하는 행과 열 선택 함수로, 아래 코드를 참고하시기 바랍니다. 

train_vec.iloc[1:5] # 1~5번째(전까지) 행 선택

train_vec.iloc[:, 1:5] # 1~5번째(전까지) 열 선택

자, 이번에는 ML 분류기를 로드해 봅니다. 여기서는 4종의 ML 분류기를 사용하며, 각각의 링크는 아래와 같습니다.

• RandomForest: https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html
• ExtraTrees: https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.ExtraTreesClassifier.html
• XGB: https://xgboost.readthedocs.io/en/stable/python/python_api.html
• LGBM: https://lightgbm.readthedocs.io/en/latest/Python-Intro.html

# ML 분류기 로드
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import ExtraTreesClassifier 
from xgboost import XGBClassifier
from lightgbm import LGBMClassifier

4종 분류기를 하나의 딕셔너리(dictionary)로 정의합니다. 앞서 사용한 리스트(list)는 인덱스 번호 순으로 값을 저장하는 반면, 딕셔너리는 키(key)와 값(value)을 한 쌍으로 저장합니다. 여기서는 이름, (모델)을 ML 분류기 딕셔너리 CLASS_MAP으로 정의했습니다.

# ML 분류기 딕셔너리: 이름, (모델)
CLASS_MAP = {
    'RF': (RandomForestClassifier()),
    'ET': (ExtraTreesClassifier()),
    'XGB': (XGBClassifier()),
    'LGBM': (LGBMClassifier())
    }

자, 다음으로 pyimpute, sklearn을 로드합니다.

# pyimpute, sklearn 로드
from pyimpute import impute
from sklearn import model_selection

이제 아래 소스코드를 통해 모델 피팅(model fitting)과 공간 예측(spatial range prediction)을 진행합니다. 단계별로 보면 일단 for in 반복문은 4종 ML 분류기를 순차적으로 처리하게 됩니다.

 

다음으로 교차 검증(cross-validation) 단계에서는 관측치를 훈련(training)과 시험(test) 목적의 두 개 모집단으로 분리하여 모델의 정확도를 계산합니다. *이 글은 소스코드 중심으로 설명합니다. 전체 데이터셋을 훈련, 검증, 시험 3종으로 분류하는 개념은 다음 링크의 이미지를 참고하시기 바랍니다: https://scikit-learn.org/stable/modules/cross_validation.html

 

아래 코드에서 사용한 'k-fold 교차 검증'은 관측치를 k개의 동일한 크기 하위 샘플로 무작위 분할한 후,  k-1개는 훈련용으로, 1개는 시험용으로 사용합니다. 모델은 k-1개 훈련 데이터를 통해 학습되고, 그 결과 모델은 나머지 1개 부분으로 검증됩니다. 이렇게 하면 모든 관측치가 훈련과 검증에 사용되면서도, 각 관측치가 정확히 검증에 1번 사용된다는 장점이 있습니다.

이미지 출처: https://en.wikipedia.org/wiki/Cross-validation_(statistics)#/media/File:K-fold_cross_validation_EN.svg

그런 다음 결과는 평균화 집계하여 단일 추정치를 생성합니다. 실습 코드에서는 각 모델마다 5-fold 교차 검증 정확도를 평균과 표준편차로 출력하였습니다. 훈련 데이터로 학습한 모델로 반응 변수의 래스터 표면을 예측하고, 그 결과를 OUTPUT 폴더에 저장합니다.

3.1. Cross-validation: evaluating estimator performance

# 모델 피팅 및 공간 예측(model fitting and spatial range prediction)
for name, (model) in CLASS_MAP.items():
    # 교차 검증(cross validation)
    k = 5 # k-fold
    kf = model_selection.KFold(n_splits=k)
    accuracy_scores = model_selection.cross_val_score(model, train_xs, train_y, cv=kf, scoring='accuracy')
    print(name + " %d-fold 교차 검증 정확도: %0.2f (+/- %0.2f)"
          % (k, accuracy_scores.mean() * 100, accuracy_scores.std() * 200))
    # 공간 예측(spatial prediction)
    model.fit(train_xs, train_y)
    os.mkdir('OUTPUT/' + name + '-IMAGES')
    impute(target_xs, model, raster_info, outdir='OUTPUT/' + name + '-IMAGES',
           class_prob=True, certainty=True)

RF 5-fold 교차 검증 정확도: 90.22 (+/- 7.68)
ET 5-fold 교차 검증 정확도: 89.74 (+/- 7.77)
XGB 5-fold 교차 검증 정확도: 90.08 (+/- 7.12)
LGBM 5-fold 교차 검증 정확도: 90.26 (+/- 7.52)

아래와 같이 4종 ML 분류기로 공간예측 결과가 매핑되었고 certainty.tif, probability_0.0.tif, probability_1.0.tif, response.tif 4개 파일이 각각의 모델 디렉터리에 생성된 것을 확인할 수 있습니다!

이중 4종 모델의 probability_1.0.tif 파일을 평균화한 후, 가시화 해보겠습니다!

import rasterio
distr_rf = rasterio.open("OUTPUT/RF-IMAGES/probability_1.0.tif").read(1)
distr_et = rasterio.open("OUTPUT/ET-IMAGES/probability_1.0.tif").read(1)
distr_xgb =  rasterio.open("OUTPUT/XGB-IMAGES/probability_1.0.tif").read(1)
distr_lgbm =  rasterio.open("OUTPUT/LGBM-IMAGES/probability_1.0.tif").read(1)
distr_averaged = (distr_rf + distr_et + distr_xgb + distr_lgbm)/4

# 종 적합성 매핑의 평균값 가시화
from pylab import plt
def plotit(x, title, cmap="Blues"):
    plt.imshow(x, cmap=cmap, interpolation='nearest')
    plt.colorbar()
    plt.title(title, fontweight = 'bold')
plotit(distr_averaged, "Joshua Tree Range, averaged", cmap="Greens")

조슈아 트리 국립공원 지역으로 확대해 본 결과는 다음과 같습니다!

# 조슈아 트리 국립공원의 종 적합성 가시화
plotit(distr_averaged[100:150, 100:150], "Joshua Tree National Park Suitability", cmap="Greens")

여기까지 QGIS와 Python을 이용한 종 분포 모델링(SDM: Species Distribution Modeling)을 시리즈 글로 정리해 봤습니다. UC 버클리 데이터 과학자이신 다니엘 퍼먼(Daniel Furman) 님이 좋은 글을 공유해 주신 덕분에 가능한 실습이었습니다. 참고로, 다니엘 퍼먼 님이 개발하신 PySDMs 파이썬용 종 분포 모델링(SDM) 패키지를 소개 드려봅니다. 'pip install PySDMs'로도 내려받을 수 있습니다.

GitHub - daniel-furman/PySDMs: An object-oriented Python class for semi-auto ML geo-classification (running on PyCaret). Compare