Google Earth Engine을 활용한 설악산국립공원의 식생 분석

안녕하세요? 이번 코드 실습은 Google Earth Engine(GEE)에서 설악산국립공원의 식생 분석을 진행해 보겠습니다. 분석 단계는 다음과 같습니다:

• GEE의 Python API를 사용하여 설악산국립공원의 경계를 추출합니다.
• 2024년 1월 동안 구름 없는 Sentinel-2 이미지를 선택합니다.
• NDVI 계산을 수행하고, 국립공원 경계 내에서 NDVI 통계치를 계산해 봅니다.

 

 

그럼 실습을 시작해볼까요?! 먼저 GEE Python API를 설치하고 인증해야 합니다.

import ee
import geemap
import pandas as pd

# Earth Engine 인증
ee.Authenticate()

# Earth Engine 초기화
ee.Initialize(project='my-project')

 

설악산국립공원 경계 추출하기

세계 보호지역 데이터베이스(WDPA: World Database on Protected Areas)는 매월 업데이트되는 보호지역에 대한 가장 최신이며 완전한 정보의 원천입니다.

 

이는 유엔 환경 프로그램의 세계 보전 모니터링 센터(UNEP-WCMC: United Nations Environment Programme's World Conservation Monitoring Centre)가 IUCN 및 세계보호지역위원회(WCPA: World Commission on Protected Areas)의 지원을 받아 관리하고 있습니다.

 

GEE에서 WDPA 폴리곤은 FeatureCollection 데이터 구조로 접근 가능하며, 설악산국립공원의 WDPA ID는 768입니다.

# 설악산국립공원 경계 추출
wdpa = ee.FeatureCollection("WCMC/WDPA/current/polygons")
seoraksan = wdpa.filter(ee.Filter.eq('WDPAID', 768))

# 선택된 보호지역의 이름 확인
wdpa_name = seoraksan.first().get('NAME').getInfo()
print("Name:", wdpa_name)

Name: Seoraksan

from ipyleaflet import TileLayer

# Vworld 배경지도 객체
vworld_base = TileLayer(
    url='https://xdworld.vworld.kr/2d/Base/service/{z}/{x}/{y}.png',
    name='Vworld Base',
    attribution='Vworld',
)

# 설악산국립공원 경계 가시화
m = geemap.Map(width="800px", height="500px")
m.add_layer(vworld_base)
m.addLayer(seoraksan, {'color': 'green'}, wdpa_name) # 레이어 추가
m.centerObject(seoraksan, 11) # 지도의 중심 설정
m # 지도 객체 출력

 

2024년 1월 동안 구름 없는 Sentinel-2 이미지 선택하기

Sentinel-2 (S2)는 전 세계적으로 5일마다 재방문하는 광역(wide-swath), 고해상도(high-resolution), 다중 스펙트럼 이미징 임무(multispectral imaging mission)입니다. S2 다중 분광 기기(MSI: MultiSpectral Instrument)는 13개의 분광 밴드를 샘플링하는데, 가시광선 및 NIR을 10미터, red edge와 SWIR을 20미터, 그리고 대기 밴드를 60미터 공간 해상도로 측정합니다. 이는 식생, 토양, 수체의 상태 및 변화를 평가하는 데 적합한 데이터를 제공합니다.

S2 데이터는 Level-1C와 Level-2A로 구분되며, GEE에서는 이를 ImageCollection 데이터 구조로 제공합니다. Level-1C는 대기 상부 반사율(Top-of-Atmosphere Reflectance), Level-2A는 지표면 반사율(Surface Reflectance)에 해당하며, 두 데이터는 모두 정사보정(Orthorectified)되어 있어, 지리적 위치가 정확하게 조정되어 있습니다. Level-1C는 대기를 통과한 후의 반사율 데이터인데 반해 Level-2A는 대기의 영향을 보정하여 얻은, 대기보정(atmospherically corrected)된 지표면의 반사율 데이터입니다. 여기서는 Level-2A를 사용하겠습니다.

 

아래 코드는 Sentinel-2 위성 이미지에서 구름과 성층운을 마스킹하여 맑은 상태의 이미지만을 추출하는 역할을 합니다. 구체적으로, 'QA60' 품질 보증(Quality Assurance) 밴드를 선택하여 구름과 성층운을 식별합니다. 비트 10과 11을 사용하여 각각 구름과 성층운을 나타내며, 이 비트들이 0으로 설정되어 있을 때만 맑은 상태로 간주합니다. 마지막으로, 이 마스크를 사용하여 구름과 성층운이 없는 영역만을 남기고, 이미지의 스케일을 조정하기 위해 결과를 10000으로 나눕니다.

def mask_s2_clouds(image):
    # QA(Quality Assurance) 밴드 사용, S2에서 구름 마스킹
    qa = image.select('QA60')
    
    # 비트 10은 구름(clouds), 11은 성층운(cirrus)
    cloud_bit_mask = 1 << 10
    cirrus_bit_mask = 1 << 11
    
    # 구름과 성층운이 0이면 맑은 상태로 간주함.
    mask = (
        qa.bitwiseAnd(cloud_bit_mask)
        .eq(0)
        .And(qa.bitwiseAnd(cirrus_bit_mask).eq(0))
    )
    
    return image.updateMask(mask).divide(10000) # 스케일링

 

Sentinel-2 위성 이미지를 선택하고 필터링하는 과정을 수행합니다. 구체적으로, 2024년 1월 1일부터 1월 31일까지의 기간에 촬영된 이미지 중에서, 지정된 seoraksan 지역을 포함하는 이미지를 필터링합니다. 그 다음, 구름이 차지하는 픽셀의 비율이 5% 미만인 이미지만을 필터링하여 선택합니다. 마지막으로, mask_s2_clouds 함수를 사용하여 각 이미지에서 구름과 성층운을 마스킹하여 맑은 상태의 이미지만을 남깁니다. 이 과정을 통해 구름이 거의 없는 맑은 상태의 이미지들만을 추출하여 분석이나 가시화에 사용할 수 있는 이미지 컬렉션을 생성합니다.

# Sentinel-2 이미지 선택 및 필터링
s2_images = (
    ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
    .filterDate("2024-01-01", "2024-01-31")
    .filterBounds(seoraksan)
    # 구름이 5% 미만인 이미지 필터링
    .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 5))
    .map(mask_s2_clouds)
)

 

필터링된 이미지의 총 개수는 4장입니다.

# 이미지 컬렉션의 이미지 개수 확인
image_count = s2_images.size()

# 이미지 개수 출력
print("Image count:", image_count.getInfo())

Image count: 4

 

필터링된 이미지 컬렉션 내의 모든 이미지들로부터 밴드별 중간값을 계산하여 새로운 단일 이미지를 생성합니다. 이는 구름이 적은 여러 날짜의 이미지들을 통합하여 더 깨끗한 대표 이미지를 얻기 위해 사용됩니다. 이미지는 RGB 컬러로 표시합니다. `B4`, `B3`, `B2`는 각각 적색, 녹색, 청색 색상에 해당하는 Sentinel-2의 밴드입니다.

# 중간값 이미지 계산
s2_image = s2_images.median()

visualization = {
    'min': 0.0,
    'max': 0.3,
    'bands': ['B4', 'B3', 'B2'],
}

m = geemap.Map(width="800px", height="500px")
m.add_layer(s2_image, visualization, 'RGB')
m.centerObject(seoraksan, 11) # 지도의 중심 설정
m # 지도 객체 출력

 

국립공원 경계 내에서 NDVI 계산 및 통계치 산출하기

프리즘을 통해 볼 수 있듯이, 태양광 스펙트럼은 많은 다른 파장으로 구성되어 있습니다. 태양광이 물체에 비추어질 때, 특정 파장은 흡수되고 다른 파장은 반사됩니다. 식물 잎의 색소인 클로로필(chlorophyll)은 광합성에 사용되는 가시광선을 강하게 흡수합니다. 반면에, 잎의 세포 구조는 근적외선을 강하게 반사합니다. 나무가 클로로필과 클로로필을 포함하는 잎을 더 많이 가질수록, 이러한 파장의 빛은 더 많이 영향을 받습니다. 과학자들은 식물이 빛과 상호작용하는 이 지식을 활용하여 지구 표면 전역의 식물이 흡수하고 반사하는 적색과 근적외선의 파장을 측정하기 위해 위성 센서를 설계함으로써 지구의 풍경을 통틀어 녹색 식생 밀도를 매핑합니다.

식물이 반사하는 적색 빛의 반사율을 근적외선 빛의 반사율에서 빼고, 그 차이를 적색과 근적외선 빛의 반사율의 합으로 나누면 과학자들이 정규식생지수(NDVI: Normalized Difference Vegetation Index)라고 부르는 값을 얻을 수 있습니다.

 

아래 코드는 Sentinel-2 위성 이미지를 사용하여 NDVI를 계산하고, 계산된 NDVI 값을 색상 팔레트에 따라 시각화하여 지도에 표시하는 과정을 수행합니다. NDVI는 Sentinel-2 이미지의 B8 밴드(NIR)와 B4 밴드(Red)를 사용하여 계산합니다.

# NDVI 계산: (NIR - Red) / (NIR + Red)
ndvi = s2_image.normalizedDifference(['B8', 'B4'])

# NDVI 색상 팔레트 정의
ndvi_palette = [
    'FE8374',  # 낮은 NDVI - 갈색
    'FED976',  # 낮은-중간 NDVI - 밝은 녹색
    'CAE23C',  # 중간 NDVI - 녹색
    '98B718',  # 중간-높은 NDVI - 진한 녹색
    '059033',  # 높은 NDVI - 매우 진한 녹색
]

# Vworld 하이브리드지도 객체
vworld_hybrid = TileLayer(
    url='https://xdworld.vworld.kr/2d/Hybrid/service/{z}/{x}/{y}.png',
    name='Vworld Hybrid',
    attribution='Vworld',
)

m = geemap.Map(width="800px", height="500px")
m.add_layer(ndvi, {'min': 0, 'max': 0.5, 'palette': ndvi_palette}, 'NDVI')
m.add_layer(vworld_hybrid)
m.centerObject(seoraksan, 11) # 지도의 중심 설정
m # 지도 객체 출력

 

마지막으로, Sentinel-2 이미지를 기반으로 계산된 NDVI에 대한 다양한 통계치를 계산하고, 이를 데이터프레임으로 변환하여 CSV 파일로 저장하는 과정을 수행합니다. 

`reduceRegion` 메서드를 사용하여 지정된 지역(`seoraksan.geometry()`)에서 NDVI의 최소값, 평균값, 중간값, 최대값, 표준편차를 계산합니다. 이는 `ee.Reducer` 객체를 사용하여 여러 통계치를 결합함으로써 한 번의 연산으로 여러 통계치를 얻습니다. `scale` 파라미터는 해상도를 10m로 설정하며, `maxPixels` 파라미터는 처리할 최대 픽셀 수(1e9는 10억 개)를 지정합니다.

# NDVI 통계치 계산 (최소값, 평균, 중간값, 최대값, 표준편차)
stats = ndvi.reduceRegion(
    reducer=ee.Reducer.min()
    .combine(reducer2=ee.Reducer.mean(), sharedInputs=True)
    .combine(reducer2=ee.Reducer.median(), sharedInputs=True)
    .combine(reducer2=ee.Reducer.max(), sharedInputs=True)
    .combine(reducer2=ee.Reducer.stdDev(), sharedInputs=True),
    geometry=seoraksan.geometry(),
    scale=10,
    maxPixels=1e9,
)

# 통계치 결과를 DataFrame으로 변환
df_stats = pd.DataFrame(
    [stats.getInfo()],
    columns=["nd_min", "nd_mean", "nd_median", "nd_max", "nd_stdDev"],
    index=["Seoraksan"],
)
df_stats.columns = ["Min", "Mean", "Median", "Max", "StdDev"]

# DataFrame을 CSV 파일로 저장하기
df_stats.to_csv('df_stats.csv', index=True)

# NDVI 통계치 출력
print(df_stats)