elapid(엘라피드): Maxent 모델 적합 및 적용 예제 소개

안녕하세요? 이번 글은 Python용 SDM 도구, elapid(엘라피드)를 사용하여 Maxent 모델의 적합 및 적용 예제(end-to-end example)를 정리해 보겠습니다. 공식 링크는 다음과 같습니다. "end-to-end"는 한 시스템이 입력에서부터 출력까지 모든 단계를 담당하는 것을 의미합니다. 이 예시는 Maxent 모델의 적합(학습)과 적용(예측) 과정을 보여줄 것입니다.

 

이 예제는 모델 학습의 가장 간단한 패턴을 보여주기 위한 것입니다. 전체 모델 학습과 평가에는 훈련/시험 데이터 나누기(train/test splits), 교차 검증(cross-validation), 피처 선택(feature selection)과 같은 요소들이 포함되어야 합니다. 이러한 내용은 여기서 다루지 않았습니다.

 

예시로 사용되는 elapid 샘플 데이터는 팔색조(Pitta nympha) 출현 기록이며, 기후 및 식생 데이터로 주석이 달려 있습니다.

elapid(엘라피드): Python용 종 분포 모델링 도구 소개

A simple Maxent model - elapid

GEODATA.zip

5.49MB

간단한 Maxent 모델

일단, 예제에서 사용되는 여러 라이브러리를 불러옵니다.

import elapid as ela

import os
import warnings
from pprint import pprint
import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio as rio
from rasterio import plot as rioplot
from sklearn import metrics
import matplotlib as mpl
import matplotlib.pyplot as plt

초기 설정으로 Matplotlib의 스타일을 "ggplot"으로 설정합니다. 이 스타일은 R의 ggplot2 패키지의 스타일을 따르는 것으로 보다 보기 좋은 그래프를 생성합니다. 경고 메시지는 무시("ignore")하도록 설정하고, elapid 라이브러리의 버전 정보를 출력합니다.

# 초기 설정
mpl.style.use('ggplot')
warnings.filterwarnings("ignore")

print(f"Notebook last run with elapid version {ela.__version__}")

Notebook last run with elapid version 1.0.1

현재 작업 디렉토리를 데이터가 있는 경로(예: 'D:/GEODATA')로 변경합니다. 이렇게 하면 코드에서 해당 디렉토리에 있는 파일들을 접근할 수 있습니다. Geopackage 파일(pitta_nympha_data.gpkg)을 읽어와 데이터를 presence라는 데이터프레임에 저장하고 "species", "year", "month", "geometry" 열만 선택하여 저장합니다. 데이터프레임의 처음 다섯 개의 행을 출력하여 데이터의 일부를 확인합니다.

# 출현 데이터
os.chdir('D:/GEODATA')
vector = "pitta_nympha_data.gpkg"

presence = gpd.read_file(vector)[["species", "year", "month", "geometry"]]
presence.head()

출현 데이터를 플롯으로 시각화해 봅니다.

# 출현 데이터 플롯
fig, ax = plt.subplots(figsize=(5,5))
sp = presence.plot(
    column='species',
    ax=ax,
    legend=True,
    legend_kwds={"bbox_to_anchor": (1, 0.5)}
)

래스터 데이터 파일들의 이름을 처리합니다. rasters에는 원래 파일명 그대로가 담긴 리스트가 저장되고, labels에는 래스터 데이터 파일의 이름만 추출하여 담긴 리스트가 저장됩니다.

# 래스터 데이터
raster_names = [
    "bio04.tif", # Temperature seasonality
    "bio05.tif", # Max temperature of warmest month
    "bio06.tif", # Min temperature of coldest month
    "bio12.tif", # Auunal precipitation
    "pct_tree_cover.tif", # Percent tree cover
]
rasters = [f"{raster}" for raster in raster_names]
labels = [raster.split(".")[0] for raster in raster_names]

print(f"rasters = {rasters}")
print(f"labels = {labels}")

rasters = ['bio04.tif', 'bio05.tif', 'bio06.tif', 'bio12.tif', 'pct_tree_cover.tif']
labels = ['bio04', 'bio05', 'bio06', 'bio12', 'pct_tree_cover']

"Temperature seasonality" 래스터 데이터와 "presence" 데이터를 함께 플롯으로 표시하여, 기온 계절변동성과 해당 위치에서 발견된 생물종의 정보를 동시에 시각적으로 확인해 봅니다.

# Temperature seasonality 플롯
fig, ax = plt.subplots(figsize=(6,6))
tss_raster = rasters[0]
with rio.open(tss_raster, 'r') as src:
    profile = src.profile.copy()
    
    # 낮은 값: 짙은 회색, 높은 값: 밝은 회색
    rioplot.show(src, ax=ax, cmap="gray")
    
    # 출현 데이터 중첩
    presence.to_crs(src.crs).plot(column='species', ax=ax, legend=True)

이번에는 배경 데이터를 생성하고 포인트 플롯으로 시각화한 뒤, 데이터의 기술 통계량(describe)을 요약하여 보여줍니다. rasters[0]는 "rasters" 리스트의 첫 번째 요소인 "Temperature seasonality" 래스터 데이터를 가리킵니다. count=10_000는 10,000개의 배경 데이터를 생성하도록 지정한 것입니다. 배경 데이터는 출현 데이터가 아닌 랜덤한 위치에서 추출되며, 모델 학습에 사용될 환경적인 배경을 나타냅니다.

배경 데이터의 기술 통계량에서 count는 배경 데이터의 개수, unique는 고유한 값의 개수, top은 최빈값(most common value), freq는 최빈값의 빈도수를 나타냅니다.

background = ela.sample_raster(rasters[0], count=10_000)

background.plot(markersize=0.5)
background.describe()

stack_geodataframes() 함수를 사용하여 출현 데이터와 배경 데이터를 병합합니다.  두 데이터프레임을 병합할 때, 출현 데이터인지 배경 데이터인지를 나타내는 클래스 레이블(class label)을 추가합니다.

다음으로 annotate() 함수를 사용하여, 병합된 데이터프레임 merged에 있는 지점 위치에서 rasters에 포함된 각 래스터 데이터들의 값을 읽어와 공변량 데이터로 추가합니다. drop_na=True는 결측치가 있는 행을 제거하여 데이터의 일관성을 유지합니다. 

# 데이터셋 병합 후 각 지점 위치에서 공변량(covariates)을 읽어오기
merged = ela.stack_geodataframes(presence, background, add_class_label=True)
annotated = ela.annotate(merged, rasters, drop_na=True, quiet=True)

데이터프레임 x에는 종속변수와 지리 데이터를 제외한 공변량 데이터가 포함되고, 시리즈 y에는 클래스 레이블(출현 데이터인지 배경 데이터인지를 나타내는 정보)이 포함됩니다. 이 데이터를 활용하여 머신러닝 모델을 학습시키고 예측할 수 있습니다.

# x/y 데이터 분할
x = annotated.drop(columns=['class', 'geometry'])
y = annotated['class']

Maxent 모델을 생성하고 데이터를 사용하여 해당 모델을 학습시키는 과정입니다. transform은 모델의 변환 방식을 설정하는 옵션이고, beta_multiplier는 정규화(regularization)를 조절하는 파라미터입니다.

훈련 데이터 x와 해당 데이터의 레이블인 종속변수 y를 사용하여, Maxent 모델이 해당 데이터에 대해 가중치를 조정하고 학습하는 과정을 수행합니다. 학습이 완료되면 모델은 주어진 데이터에 대해 적절한 예측을 할 수 있도록 준비됩니다.

# 모델 학습
model = ela.MaxentModel(transform='cloglog', beta_multiplier=2.0)
model.fit(x, y)

학습된 Maxent 모델인 model을 사용하여 입력 데이터 x에 대한 예측값 ypred를 계산합니다. 이 예측값은 주어진 입력 데이터에 대한 모델의 출력으로서, 이진 분류 문제에서는 0과 1 사이의 확률값으로 표현됩니다. roc_auc_score 함수는 AUC(Area Under the Receiver Operating Characteristic Curve) 값을 계산하는 함수입니다. AUC는 Receiver Operating Characteristic(ROC) 곡선 아래 영역을 나타내는 지표로 0과 1 사이의 값을 가지며 높을수록 모델의 예측 성능이 좋다는 것을 의미합니다. 

# 학습 성능 평가
ypred = model.predict(x)
auc = metrics.roc_auc_score(y, ypred)
print(f"Training AUC score: {auc:0.3f}")

Training AUC score: 0.893

model 변수에 저장된 Maxent 모델을 'maxent-model.ela'라는 파일 이름으로 디스크에 저장합니다. 이렇게 저장된 파일은 나중에 필요할 때 다시 불러와서 모델을 재사용하거나 예측에 활용할 수 있습니다.

# 적합된 모델을 디스크에 저장
ela.save_object(model, 'maxent-model.ela')

Maxent 모델의 부분 의존도 그래프(PDP: Partial Dependence Plot)를 생성해 봅니다. 부분 의존도 그래프는 모델의 예측 결과가 공변량(독립변수)의 범위에 따라 어떻게 변하는지를 시각적으로 보여주는 도구입니다. PDP는 특정 변수의 영향력을 직관적으로 이해하는 데에 도움을 줍니다.

# 부분 의존도 그래프(PDP: Partial Dependence Plot)
# 모델의 예측이 공변량의 범위에 따라 어떻게 변하는지 보여준다.
fig, ax = model.partial_dependence_plot(x, labels=labels, dpi=100)

Maxent 모델의 예측 결과를 GeoTIFF 파일 형식으로 디스크에 저장하고, 저장된 결과를 메모리에 다시 불러옵니다.

# 모델의 예측결과를 디스크에 저장
output_raster = 'maxent-predictions.tif'
ela.apply_model_to_rasters(model, rasters, output_raster, quiet=True)

# 메모리로 불러오기
with rio.open(output_raster, 'r') as src:
    pred = src.read(1, masked=True)

적합성 예측(suitability predictions) 결과를 시각화하여, 픽셀별로 적합성 확률을 색상으로 표현한 그래프를 나타냅니다. 여기까지 elapid(엘라피드)를 사용하여 Maxent 모델의 적합 및 적용 예제(end-to-end example)를 정리해 봤습니다.

# 적합성 예측(suitability predictions) 그래프로 표현
fig, ax = plt.subplots(1, 1, figsize=(9, 6), dpi=100)
plot = ax.imshow(pred, vmin=0, vmax=1, cmap='GnBu')
ax.set_title('$Pitta\ nympha$ suitability')
ax.set_xticks([])
ax.set_yticks([])
cbar = plt.colorbar(plot, ax=ax, label="relative occurrence probability", pad=0.04)
plt.tight_layout()

1911_elapid(엘라피드)_Maxent 모델 적합 및 적용 예제 소개.ipynb

0.48MB