Earth Engine & Geemap: 종 분포 모델링(SDM) 구현 (2) 예측 변수 추가

안녕하세요? 이번 글은 시리즈 글의 일부로 Google Earth Engine과 geemap을 이용한 종 분포 모델링(SDM: Species Distribution Modeling) 구현 방법을 소개해 보겠습니다. 이 글의 내용은 스미소니언 보전생물연구소 연구진 분들이 제공한 JavaScript 소스코드를 Python으로 변환하여 수정, 보완한 것입니다.

 

Crego, R. D., Stabach, J. A., & Connette, G. (2022). Implementation of species distribution models in Google Earth Engine. Diversity and Distributions, 28, 904–916. https://doi.org/10.1111/ddi.13491

 

사례 연구 1: 출현-전용 데이터를 이용한 팔색조(Pitta nympha) 서식지 적합성 및 예측 분포 모델링

2. 관심 영역 정의

다음 단계는 관심 영역(AOI: Area Of Interest)의 범위를 정의하는 것입니다. 여기서는 전체 출현 데이터를 포함하는 경계 상자(bounding box) 주위에 100km 버퍼를 정의했습니다. 버퍼 거리는 미터(m) 단위입니다.

# AOI 정의
AOI = Data.geometry().bounds().buffer(distance=50000, maxError=1000)

# 연구지역 경계를 지도에 추가
outline = ee.Image().byte().paint(featureCollection=AOI, color=1, width=3)
Map.addLayer(outline, {'palette': 'FF0000'}, "Study Area")
Map.centerObject(Data.geometry(), 6)
Map

3. 예측 변수 추가

WorldClim BIO 변수 V1: BIO 변수는 "WORLDCLIM/V1/BIO" 이미지를 로드합니다. 이는 WorldClim 데이터셋에서 가져온 1km 해상도의 총 19개 기후 변수로 구성된 이미지입니다. BIO 변수는 기온, 강수량 등과 같은 다양한 기후 변수 정보를 가지고 있습니다.

NASA SRTM Digital Elevation: Terrain 변수는 "USGS/SRTMGL1_003" 이미지를 통해 NASA의 SRTM(Shuttle Radar Topography Mission) 디지털 고도 데이터를 로드합니다. 이 데이터는 지구 표면의 고도 정보를 제공하며 30m 해상도로 제공됩니다.

 

MOD44B.006 Terra MODIS Vegetation Continuous Fields Yearly Global (Percent Tree Cover): MODIS 변수는 "MODIS/006/MOD44B" 이미지 컬렉션을 로드합니다. 이 컬렉션에는 2003년부터 2020년까지의 Terra 위성 MODIS 데이터를 기반으로 한 수목 피복 확률 정보가 포함되어 있습니다. filterDate를 사용하여 2003년 1월 1일부터 2020년 12월 31일 사이의 데이터를 필터링하고, 'Percent_Tree_Cover' 밴드를 선택하여 해당 기간 동안의 수목 피복 확률 데이터를 250m급 해상도로 가져옵니다. 그런 다음 median()을 사용하여 해당 기간 동안 중간값을 계산하여 MedianPTC 변수에 저장합니다.

# WorldClim BIO Variables V1(1㎞): 총 19개 기후 변수
BIO = ee.Image("WORLDCLIM/V1/BIO")

# NASA SRTM Digital Elevation 30m: 디지털 고도 데이터
Terrain = ee.Algorithms.Terrain(ee.Image("USGS/SRTMGL1_003"))

# MOD44B.006 Terra MODIS Vegetation Continuous Fields Yearly Global 250m
# 수목 피복 확률(percent tree cover): 2003~2020년 사이 중간값
MODIS = ee.ImageCollection("MODIS/006/MOD44B")
MedianPTC = MODIS.filterDate('2003-01-01', '2020-12-31').select(['Percent_Tree_Cover']).median()

predictors 변수는 BIO, Terrain, 그리고 MedianPTC 변수들을 각각의 밴드로 가지는 단일 다중-밴드 이미지로 결합합니다. Terrain 변수에서 'elevation' 밴드를 선택하여 0보다 큰 값으로 해양 픽셀을 마스킹하고 관심 영역 내에 해당하는 이미지만 클리핑합니다.

# 밴드들을 단일 다중-밴드 이미지로 결합
predictors = BIO.addBands(Terrain).addBands(MedianPTC)

# 예측 변수 이미지에서 해양 픽셀 마스킹
watermask =  Terrain.select('elevation').gt(0) # 수역 마스크 생성
predictors = predictors.updateMask(watermask).clip(AOI) # 해양 픽셀 마스킹 & 관심 영역 클리핑

예측 변수 이미지에서 임의로 5,000개 지점을 선택하고, 해당 지점에서 다중-밴드 예측 이미지의 공변량 값을 추출합니다.

# 5000개 임의 지점에서 다중-밴드 예측 이미지 간 공변량 값 추출
DataCor = predictors.sample(scale=GrainSize, numPixels=5000, geometries=True) # 5000개 임의 지점 생성
PixelVals = predictors.sampleRegions(collection=DataCor, scale=GrainSize, tileScale=16) # 공변량 값 추출

Earth Engine에서 추출한 지점별 공변량 값 데이터를 판다스 데이터프레임 형식으로 변환하고 처음 5개 행을 확인합니다.

# 공변량 값 ee를 df로 변환
PixelVals_df = geemap.ee_to_pandas(PixelVals)
PixelVals_df.head()

PixelVals_df의 컬럼 이름들을 추출하여 columns 변수에 저장합니다.

columns = PixelVals_df.columns
columns

Index(['elevation', 'bio08', 'bio19', 'bio09', 'slope', 'bio10', 'hillshade',
       'aspect', 'bio06', 'bio17', 'bio07', 'bio18', 'bio04', 'bio15', 'bio05',
       'bio16', 'bio02', 'bio13', 'bio03', 'bio14', 'bio11', 'PTC', 'bio01',
       'bio12'],
      dtype='object')

스피어만(Spearman) 상관계수를 사용하여 주어진 예측 변수들 간의 상관관계를 계산하고, 이를 시각화한 히트맵을 출력해 봅니다.

def plot_correlation_heatmap(dataframe, h_size=10):
    # 스피어만 상관계수 계산
    correlation_matrix = dataframe.corr(method="spearman")

    # 히트맵 그리기
    plt.figure(figsize=(h_size, h_size-2))
    plt.imshow(correlation_matrix, cmap='coolwarm', interpolation='nearest')

    # 히트맵에 수치 표시
    for i in range(correlation_matrix.shape[0]):
        for j in range(correlation_matrix.shape[1]):
            plt.text(j, i, f"{correlation_matrix.iloc[i, j]:.2f}",
                     ha='center', va='center', color='white', fontsize=8)  # fontsize 조정

    columns = dataframe.columns.tolist()
    plt.xticks(range(len(columns)), columns, rotation=90)
    plt.yticks(range(len(columns)), columns)
    plt.title("Variables correlation matrix")
    plt.colorbar(label="Spearman Correlation")
    # plt.savefig('correlation_heatmap_plot.png')
    plt.show()

# 변수 상관행렬 히트맵
plot_correlation_heatmap(PixelVals_df)

상관관계가 높은 변수들을 함께 모델에 포함하면 다중공선성(Multicollinearity) 문제가 발생할 수 있습니다. 다중공선성은 모델 내에서 독립 변수들 간에 강한 선형 관계가 있을 때 발생하는 현상으로, 이로 인해 모델의 계수 추정치가 불안정해지고 해석이 어려워집니다.

 

모델은 각 변수가 종속 변수에 어떤 영향을 미치는지 정확하게 분리하기 어려워지고, 이는 계수(가중치) 추정의 불안정성을 야기하며 변수들의 계수가 원래의 의미와는 다른 값으로 추정될 수 있습니다. 이러한 불안정성은 모델의 신뢰도를 떨어뜨리며, 새로운 데이터에 대한 예측 또는 해석을 어렵게 만들 수 있습니다.

 

VIF(Variance Inflation Factor)는 다중공선성을 평가하고 변수를 선택하는데 사용되는 통계적 지표입니다. 각 독립 변수의 다른 독립 변수에 대한 선형 관계의 정도를 나타내며, 높은 VIF 값은 다중공선성의 증거일 수 있습니다. 일반적으로 VIF 값은 5또는 10보다 크면 해당 변수가 다른 변수들과 강한 상관관계를 가지며 모델의 안정성과 해석력을 저해할 가능성이 높다는 것을 나타냅니다. 여기서는 10보다 작은 값을 기준으로 변수를 선택했습니다.

def filter_variables_by_vif(dataframe, threshold=10):
    original_columns = dataframe.columns.tolist()
    remaining_columns = original_columns[:]
    
    while True:
        vif_data = dataframe[remaining_columns]
        vif_values = [variance_inflation_factor(vif_data.values, i) for i in range(vif_data.shape[1])]
        
        max_vif_index = vif_values.index(max(vif_values))
        max_vif = max(vif_values)
        
        if max_vif < threshold:
            break
        
        print(f"Removing '{remaining_columns[max_vif_index]}' with VIF {max_vif:.2f}")
        
        del remaining_columns[max_vif_index]
    
    filtered_data = dataframe[remaining_columns]
    return filtered_data

# VIF(Variance Inflation Factor, 분산 팽창 요인)로 필터링 수행
filtered_PixelVals_df = filter_variables_by_vif(PixelVals_df)
bands = filtered_PixelVals_df.columns.tolist()
predictors = predictors.select(bands)
print('Bands:', bands)

Removing 'bio06' with VIF inf
Removing 'bio04' with VIF 213143.75
Removing 'bio10' with VIF 68437.04
Removing 'bio17' with VIF 48555.55
Removing 'bio05' with VIF 32374.85
Removing 'bio07' with VIF 22348.86
Removing 'bio01' with VIF 9725.41
Removing 'bio16' with VIF 3785.62
Removing 'bio03' with VIF 2515.02
Removing 'bio18' with VIF 1791.83
Removing 'bio08' with VIF 1383.23
Removing 'bio12' with VIF 618.77
Removing 'bio19' with VIF 488.40
Removing 'bio15' with VIF 340.60
Removing 'hillshade' with VIF 120.15
Removing 'bio13' with VIF 49.94
Removing 'bio11' with VIF 35.68
Removing 'bio02' with VIF 15.55
Bands: ['elevation', 'bio09', 'slope', 'aspect', 'bio14', 'PTC']

# 변수 상관행렬 히트맵
plot_correlation_heatmap(filtered_PixelVals_df, h_size=6)

6개 예측 변수 데이터(['elevation', 'bio09', 'slope', 'aspect', 'bio14', 'PTC'])에 대해 최소/최대값, 팔레트를 지정하고 레이어와 컬러바를 추가하여 시각화한 결과입니다.

cm.plot_colormaps(width=8.0, height=0.2)

cm.plot_colormap('terrain', width=8.0, height=0.2, orientation='horizontal')

# elevation
Map = geemap.Map()
vis_params = {'bands':['elevation'], 'min': 0, 'max': 1800, 'palette': cm.palettes.terrain}
Map.addLayer(predictors, vis_params, 'elevation')
Map.add_colorbar(vis_params, label="Elevation (m)", orientation="vertical", layer_name="elevation")
Map.centerObject(Data.geometry(), 7)
Map

# bio09 최소값 & 최대값
min_val = predictors.select("bio09").multiply(0.1).reduceRegion(reducer=ee.Reducer.min(), scale=1000).getInfo()
max_val = predictors.select("bio09").multiply(0.1).reduceRegion(reducer=ee.Reducer.max(), scale=1000).getInfo()
print("Minimum Value:", min_val)
print("Maximum Value:", max_val)

Minimum Value: {'bio09': -10.9}
Maximum Value: {'bio09': 13.5}

# bio09
Map = geemap.Map()
vis_params = {'min': -11, 'max': 14, 'palette': cm.palettes.hot}
Map.addLayer(predictors.select("bio09").multiply(0.1), vis_params, 'bio09')
Map.add_colorbar(vis_params, label="Mean temperature of driest quarter (℃)", orientation="vertical", layer_name="bio09")
Map.centerObject(Data.geometry(), 7)
Map

# slope
Map = geemap.Map()
vis_params = {'bands':['slope'], 'min': 0, 'max': 25, 'palette': cm.palettes.RdYlGn_r}
Map.addLayer(predictors, vis_params, 'slope')
Map.add_colorbar(vis_params, label="Slope", orientation="vertical", layer_name="slope")
Map.centerObject(Data.geometry(), 7)
Map

# aspect
Map = geemap.Map()
vis_params = {'bands':['aspect'], 'min': 0, 'max': 360, 'palette': cm.palettes.rainbow}
Map.addLayer(predictors, vis_params, 'aspect')
Map.add_colorbar(vis_params, label="Aspect", orientation="vertical", layer_name="aspect")
Map.centerObject(Data.geometry(), 7)
Map

# bio14 최소값 & 최대값
# min_val = predictors.select("bio14").reduceRegion(reducer=ee.Reducer.min(), scale=1000).getInfo()
# max_val = predictors.select("bio14").reduceRegion(reducer=ee.Reducer.max(), scale=1000).getInfo()
# print("Minimum Value:", min_val)
# print("Maximum Value:", max_val)

# bio14
Map = geemap.Map()
vis_params = {'bands':['bio14'], 'min': 0, 'max': 90, 'palette': cm.palettes.Blues}
Map.addLayer(predictors, vis_params, 'bio14')
Map.add_colorbar(vis_params, label="Precipitation of driest month (mm)", orientation="vertical", layer_name="bio14")
Map.centerObject(Data.geometry(), 7)
Map

# PTC
Map = geemap.Map()
vis_params = {'bands':['PTC'], 'min': 1, 'max': 100, 'palette': ['bbe029', '0a9501', '074b03']}
Map.addLayer(predictors, vis_params, 'PTC')
Map.add_colorbar(vis_params, label="Percent Tree Cover (%)", orientation="vertical", layer_name="PTC")
Map.centerObject(Data.geometry(), 7)
Map

1917_ee_geemap_종 분포 모델링(SDM) 구현 (2) 예측 변수 추가.ipynb

0.88MB