Google Earth Engine: 백두대간보호지역의 NPP(순일차생산량) 계산하기

안녕하세요? 이번 글은 Google Earth Engine을 사용하여 백두대간보호지역의 NPP(순일차생산량)을 계산해 보겠습니다. 백두대간보호지역은 하나의 예시이고, 다른 지역에 대해서도 NPP를 계산할 수 있는 방법을 소개해 봅니다.

 

사용할 MOD17A3HGF V6.1 이미지 컬렉션은 500m 픽셀 해상도로 연간 총일차생산량(GPP: Gross Primary Productivity)과 순일차생산량(NPP: Net Primary Productivity)에 대한 정보를 제공합니다.

연간 NPP는 해당 연도의 8일 주기 PSN(Net Photosynthesis) 데이터(MOD17A2H)의 합에서 도출됩니다. PSN 값은 총일차생산량에서 유지 호흡(MR: Maintenance Respiration)을 뺀 값(GPP-MR)입니다.

MOD17A3HGF.061: Terra Net Primary Production Gap-Filled Yearly Global 500m  |  Earth Engine Data Catalog  |  Google for Deve

 

먼저 Google Earth Engine에서 백두대간보호지역 경계를 호출해 보겠습니다. WDPA ID를 가지고 해당 경계를 불러올 수 있습니다.

Protected Planet | Baekdudaegan Mountains Reserve

# 백두대간보호지역
wdpa = ee.FeatureCollection("WCMC/WDPA/current/polygons")
wdpa = wdpa.filter(ee.Filter.eq('WDPAID', 478389))

# 1: 외곽선의 픽셀 값, 2: 외곽선의 두께(픽셀 단위)
wdpa_raster = ee.Image().paint(wdpa, 1, 2)

wdpa_vis = {
    'palette': ['#FFFF00'], # 노란색 적용
    'opacity': 1
}

m = geemap.Map(width="800px", height="400px",
               basemap='Esri.WorldImagery')
m.addLayer(wdpa_raster, wdpa_vis, '백두대간보호지역')
m.centerObject(wdpa, 6)
m

 

2023년 NPP를 호출해 봅니다. 연간 NPP이므로 레이어는 1장입니다.

dataset = (
    ee.ImageCollection('MODIS/061/MOD17A3HGF')
    .filterDate("2023-01-01", "2023-12-31")
    .filterBounds(wdpa)
)

layer_count = dataset.size().getInfo()
print("Layer count:", layer_count)

 

축척계수, 좌표계, 해상도를 설정하고 백두대간보호지역과 함께 NPP를 가시화해 봅니다.

npp = dataset.select('Npp').first().multiply(0.0001)
crs = ee.Projection('EPSG:5179')
npp = npp.reproject(crs=crs, scale=500)

npp_vis = {
    'min': 0,
    'max': 1.9,
    'palette': ['bbe029', '0a9501', '074b03']
}

m = geemap.Map(width="800px", height="400px",
               basemap='Esri.WorldImagery')
m.addLayer(npp, npp_vis, 'NPP')
m.addLayer(wdpa_raster, wdpa_vis, '백두대간보호지역')
m.centerObject(wdpa, 6)
m

 

백두대간보호지역의 2023년 NPP를 Mg으로 출력해 보겠습니다. 해상도는 500m로 유지하고 단위변환을 위한 계산을 수행합니다. 값은 2,213,523.68 Mg으로 제시됩니다.

# NPP 픽셀 값 합계 계산
npp_sum = npp.reduceRegion(
    reducer=ee.Reducer.sum(),
    geometry=wdpa.geometry(),
    scale=500,
    maxPixels=1e13
).get('Npp').getInfo()

total_npp = npp_sum * 250000 * 0.001

# 결과 출력
print(f"{total_npp:.2f} Mg")

 

이번에는 2001년부터 2023년까지 백두대간보호지역의 연간 NPP를 전부 추출해 보겠습니다. 반복문을 통해 각 연도에 대한 NPP를 리스트로 저장하는 방식입니다.

# 연도 리스트 생성
years = list(range(2001, 2024))

# 각 연도에 대한 NPP 값 저장할 리스트 초기화
npp_values = []

# 각 연도에 대해 NPP 값을 계산
for year in years:
    dataset = (
        ee.ImageCollection('MODIS/061/MOD17A3HGF')
        .filterDate(f"{year}-01-01", f"{year}-12-31")
        .filterBounds(wdpa)
    )
    
    # 첫 번째 이미지를 선택하고 필요한 연산 수행
    npp = dataset.select('Npp').first().multiply(0.0001)
    npp = npp.reproject(crs=crs, scale=500)
    
    # NPP 픽셀 값 합계 계산
    npp_sum = npp.reduceRegion(
        reducer=ee.Reducer.sum(),
        geometry=wdpa.geometry(),
        scale=500,
        maxPixels=1e13
    ).get('Npp').getInfo()
    
    # 총 NPP 계산 및 리스트에 추가
    total_npp = npp_sum * 250000 * 0.001
    npp_values.append(total_npp)

 

각 연도의 NPP를 출력해 보겠습니다.

# 각 연도의 NPP 값 출력
for year, npp_value in zip(years, npp_values):
    print(f"{year}: {npp_value:.2f} Mg")

2001: 2095738.30 Mg
2002: 2241697.64 Mg
2003: 2145454.74 Mg
2004: 2180237.28 Mg
2005: 2066550.91 Mg
2006: 2129783.49 Mg
2007: 2158312.93 Mg
2008: 2165888.99 Mg
2009: 2214339.63 Mg
2010: 1946082.52 Mg
2011: 1952524.44 Mg
2012: 2163056.21 Mg
2013: 2062665.49 Mg
2014: 2178660.25 Mg
2015: 2132293.20 Mg
2016: 1941293.69 Mg
2017: 1953672.64 Mg
2018: 2120833.62 Mg
2019: 2171249.27 Mg
2020: 2272044.38 Mg
2021: 2191090.26 Mg
2022: 2105580.53 Mg
2023: 2213523.68 Mg

 

마지막 연도와 첫번째 연도의 NPP 변화율은 5.62%로 확인됩니다.

# 전체 데이터에서의 변화율 계산
total_change_rate = ((npp_values[-1] - npp_values[0]) / npp_values[0]) * 100
print(f"Total NPP Change Rate from {years[0]} to {years[-1]}: {total_change_rate:.2f}%")

 

NPP 변화를 그래프로 표현해본 결과는 다음과 같습니다.

# 컬러맵 적용을 위한 색상 생성
norm = plt.Normalize(min(npp_values), max(npp_values))
colors = cm.viridis(norm(npp_values))

# NPP 값을 백만 단위로 변환
npp_values_million = [value / 1e6 for value in npp_values]

# 그래프 그리기
plt.figure(figsize=(10, 3))
plt.bar(years, npp_values_million, color=colors)
plt.xlabel('Year')
plt.ylabel('Annual NPP (million Mg)')
plt.title('Annual NPP from 2001 to 2023')
plt.xticks(
    ticks=[2000, 2005, 2010, 2015, 2020],
    labels=["2000", "2005", "2010", "2015", "2020"],
)  # x축 5년 주기 설정
plt.grid(axis='y', linestyle='--', alpha=0.7)

# 컬러바 추가
sm = plt.cm.ScalarMappable(cmap='viridis', norm=norm)
sm.set_array([])
plt.colorbar(sm)

plt.show()