3장_위성기반_식생변화_분석.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/osgeokr/GEE-PAM-Book/blob/main/3%EC%9E%A5_%EC%9C%84%EC%84%B1%EA%B8%B0%EB%B0%98_%EC%8B%9D%EC%83%9D%EB%B3%80%ED%99%94_%EB%B6%84%EC%84%9D.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c98405ab-ccc0-4f02-8154-2149a5e4e63e",
      "metadata": {
        "id": "c98405ab-ccc0-4f02-8154-2149a5e4e63e"
      },
      "source": [
        "# 3장. 위성기반 식생변화 분석\n",
        "\n",
        "## 3.1. 설악산국립공원 식생지수 분석"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "7a0e8ae3-bb8a-4d44-9c90-49d7aa2f3fcf",
      "metadata": {
        "id": "7a0e8ae3-bb8a-4d44-9c90-49d7aa2f3fcf"
      },
      "outputs": [],
      "source": [
        "import ee\n",
        "import geemap\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "from ipyleaflet import TileLayer\n",
        "from google.colab import files\n",
        "\n",
        "# Earth Engine 인증\n",
        "ee.Authenticate()\n",
        "\n",
        "# Earth Engine 초기화\n",
        "ee.Initialize(project='ee-foss4g')"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2f4dad2c-01ce-48ff-ad50-70f22ce724e2",
      "metadata": {
        "id": "2f4dad2c-01ce-48ff-ad50-70f22ce724e2"
      },
      "source": [
        "### 3.1.1. 설악산국립공원 경계 추출"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e2b0ce1e-8d42-48b4-a070-ade5330fc507",
      "metadata": {
        "id": "e2b0ce1e-8d42-48b4-a070-ade5330fc507"
      },
      "outputs": [],
      "source": [
        "# 설악산국립공원 경계 추출\n",
        "wdpa = ee.FeatureCollection(\"WCMC/WDPA/current/polygons\")\n",
        "seoraksan = wdpa.filter(ee.Filter.eq('WDPAID', 768))\n",
        "\n",
        "# 선택된 보호지역 이름 확인\n",
        "wdpa_name = seoraksan.first().get('NAME').getInfo()\n",
        "print(\"Name:\", wdpa_name)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e0c364d3-0a35-44e9-82b5-20e1370eab8f",
      "metadata": {
        "id": "e0c364d3-0a35-44e9-82b5-20e1370eab8f"
      },
      "outputs": [],
      "source": [
        "# Vworld 배경지도 객체\n",
        "vworld_base = TileLayer(\n",
        "    url='https://xdworld.vworld.kr/2d/Base/service/{z}/{x}/{y}.png',\n",
        "    name='Vworld Base', attribution='Vworld',\n",
        ")\n",
        "\n",
        "# 설악산국립공원 경계\n",
        "m = geemap.Map(width=\"800px\", height=\"400px\")\n",
        "m.add_layer(vworld_base)\n",
        "m.addLayer(seoraksan, {'color': 'green'}, wdpa_name)\n",
        "m.centerObject(seoraksan, 11)\n",
        "m"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import requests\n",
        "import json\n",
        "\n",
        "# GitHub 파일 다운로드\n",
        "url = 'https://github.com/osgeokr/GEE-PAM-Book/raw/main/SEORAKSAN.geojson'\n",
        "response = requests.get(url)\n",
        "geojson = response.json()\n",
        "\n",
        "# 파일 저장\n",
        "with open('SEORAKSAN.geojson', 'wb') as file:\n",
        "    file.write(response.content)\n",
        "files.download('SEORAKSAN.geojson')\n",
        "\n",
        "# GeoJSON 데이터를 ee.FeatureCollection으로 변환\n",
        "seoraksan = ee.FeatureCollection(geojson)\n",
        "print(seoraksan.getInfo())"
      ],
      "metadata": {
        "id": "eEC9TkEteQVB"
      },
      "id": "eEC9TkEteQVB",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "id": "5adc2e2d-b995-462b-85d1-2079671efd1f",
      "metadata": {
        "id": "5adc2e2d-b995-462b-85d1-2079671efd1f"
      },
      "source": [
        "### 3.1.2. Sentinel-2 이미지 선택 및 필터링"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "fa78b221-0bc4-4623-b204-d8cd057844bb",
      "metadata": {
        "id": "fa78b221-0bc4-4623-b204-d8cd057844bb"
      },
      "outputs": [],
      "source": [
        "def mask_s2_clouds(image):\n",
        "    # QA(Quality Assurance) 밴드 사용, S2에서 구름 마스킹\n",
        "    qa = image.select('QA60')\n",
        "\n",
        "    # 비트 10은 구름(clouds), 11은 성층운(cirrus)\n",
        "    cloud_bit_mask = 1 << 10\n",
        "    cirrus_bit_mask = 1 << 11\n",
        "\n",
        "    # 구름과 성층운이 0이면 맑은 상태로 간주함.\n",
        "    mask = (\n",
        "        qa.bitwiseAnd(cloud_bit_mask)\n",
        "        .eq(0)\n",
        "        .And(qa.bitwiseAnd(cirrus_bit_mask).eq(0))\n",
        "    )\n",
        "\n",
        "    return image.updateMask(mask).divide(10000) # 스케일링"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "1baeb159-9f41-464d-ab59-df83b3d944c8",
      "metadata": {
        "id": "1baeb159-9f41-464d-ab59-df83b3d944c8"
      },
      "outputs": [],
      "source": [
        "# Sentinel-2 이미지 선택 및 필터링\n",
        "s2_images = (\n",
        "    ee.ImageCollection(\"COPERNICUS/S2_SR_HARMONIZED\")\n",
        "    .filterDate(\"2024-01-01\", \"2024-01-31\")\n",
        "    .filterBounds(seoraksan)\n",
        "    # 구름이 5% 미만인 이미지 필터링\n",
        "    .filter(ee.Filter.lt(\"CLOUDY_PIXEL_PERCENTAGE\", 5))\n",
        "    .map(mask_s2_clouds)\n",
        ")\n",
        "\n",
        "# 이미지 컬렉션의 이미지 개수 확인\n",
        "image_count = s2_images.size()\n",
        "\n",
        "# 이미지 개수 출력\n",
        "print(\"Image count:\", image_count.getInfo())"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# 중간값 이미지 계산\n",
        "s2_image = s2_images.median()"
      ],
      "metadata": {
        "id": "CTW82bUGH2rU"
      },
      "id": "CTW82bUGH2rU",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "aab62fbb-1f52-422c-a1dd-7c82aae485af",
      "metadata": {
        "id": "aab62fbb-1f52-422c-a1dd-7c82aae485af"
      },
      "outputs": [],
      "source": [
        "# TCC\n",
        "vis_params = {\n",
        "    'min': 0.0,\n",
        "    'max': 0.3,\n",
        "    'bands': ['B4', 'B3', 'B2'],\n",
        "}\n",
        "\n",
        "m = geemap.Map(width=\"800px\", height=\"400px\")\n",
        "m.add_layer(s2_image, vis_params, 'TCC')\n",
        "m.centerObject(seoraksan, 11)\n",
        "m"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#FCC\n",
        "vis_params = {\n",
        "    'min': 0.0,\n",
        "    'max': 0.3,\n",
        "    'bands': ['B8', 'B4', 'B3'],\n",
        "}\n",
        "\n",
        "m = geemap.Map(width=\"800px\", height=\"400px\")\n",
        "m.add_layer(s2_image, vis_params, 'FCC')\n",
        "m.centerObject(seoraksan, 11)\n",
        "m"
      ],
      "metadata": {
        "id": "34AZOwfWHpkD"
      },
      "id": "34AZOwfWHpkD",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# 밴드 이름 목록\n",
        "band_names = s2_image.bandNames()\n",
        "print(band_names.getInfo())"
      ],
      "metadata": {
        "id": "P68XnQ79gLvm"
      },
      "id": "P68XnQ79gLvm",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# B3, B4, B8 밴드 선택\n",
        "fcc_image = s2_image.select(['B3', 'B4', 'B8']).clip(seoraksan.geometry())"
      ],
      "metadata": {
        "id": "5G6gUT17j9aV"
      },
      "id": "5G6gUT17j9aV",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import time\n",
        "\n",
        "# 이미지를 Google Drive에 내보내기\n",
        "task = ee.batch.Export.image.toDrive(\n",
        "    image=fcc_image,\n",
        "    description='S2_FCC_202401',\n",
        "    folder='export',\n",
        "    scale=30,  # 이미지의 해상도\n",
        "    region=seoraksan.geometry(),  # 내보낼 영역\n",
        "    maxPixels=1e13,\n",
        "    fileFormat='GeoTIFF'\n",
        ")\n",
        "\n",
        "# 내보내기 작업 시작\n",
        "task.start()\n",
        "\n",
        "# 내보내기 작업 상태 확인\n",
        "print('Export task started. Checking status...')\n",
        "\n",
        "while True:\n",
        "    status = task.status()\n",
        "    state = status['state']\n",
        "    print('Polling for task (id: {}). Status: {}'.format(task.id, state))\n",
        "    if state in ['COMPLETED', 'FAILED']:\n",
        "        break\n",
        "    time.sleep(30)  # 30초 간격으로 상태 확인\n",
        "\n",
        "# 완료 후 상태 출력\n",
        "print('Task completed. Final status:')\n",
        "print(status)"
      ],
      "metadata": {
        "id": "LNcG865ghnKq"
      },
      "id": "LNcG865ghnKq",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "id": "314da205-feb7-44cb-b6c8-ef249eb3dae4",
      "metadata": {
        "id": "314da205-feb7-44cb-b6c8-ef249eb3dae4"
      },
      "source": [
        "### 3.1.3. 국립공원 경계 내에서 NDVI 계산 및 통계치 산출"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# NDVI 계산: (NIR - Red) / (NIR + Red)\n",
        "ndvi = s2_image.normalizedDifference(['B8', 'B4'])"
      ],
      "metadata": {
        "id": "Kd3J5FYTJhQU"
      },
      "id": "Kd3J5FYTJhQU",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b692709c-6e45-479c-84c2-830e408a2c6f",
      "metadata": {
        "id": "b692709c-6e45-479c-84c2-830e408a2c6f"
      },
      "outputs": [],
      "source": [
        "# NDVI 색상 팔레트 정의\n",
        "ndvi_palette = [\n",
        "    'FE8374',  # 낮은 NDVI - 갈색\n",
        "    'FED976',  # 낮은-중간 NDVI - 밝은 녹색\n",
        "    'CAE23C',  # 중간 NDVI - 녹색\n",
        "    '98B718',  # 중간-높은 NDVI - 진한 녹색\n",
        "    '059033',  # 높은 NDVI - 매우 진한 녹색\n",
        "]\n",
        "\n",
        "# Vworld 하이브리드지도 객체\n",
        "vworld_hybrid = TileLayer(\n",
        "    url='https://xdworld.vworld.kr/2d/Hybrid/service/{z}/{x}/{y}.png',\n",
        "    name='Vworld Hybrid',\n",
        "    attribution='Vworld',\n",
        ")\n",
        "\n",
        "m = geemap.Map(width=\"800px\", height=\"400px\")\n",
        "m.add_layer(ndvi, {'min': 0, 'max': 0.5, 'palette': ndvi_palette}, 'NDVI')\n",
        "m.add_layer(vworld_hybrid)\n",
        "m.centerObject(seoraksan, 11)\n",
        "m"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b7f10d4c-5d54-49d9-a2fc-0cb74c60e8d3",
      "metadata": {
        "id": "b7f10d4c-5d54-49d9-a2fc-0cb74c60e8d3"
      },
      "outputs": [],
      "source": [
        "# NDVI 통계치 계산 (최소값, 평균, 중간값, 최대값, 표준편차)\n",
        "stats = ndvi.reduceRegion(\n",
        "    reducer=ee.Reducer.min()\n",
        "    .combine(reducer2=ee.Reducer.mean(), sharedInputs=True)\n",
        "    .combine(reducer2=ee.Reducer.median(), sharedInputs=True)\n",
        "    .combine(reducer2=ee.Reducer.max(), sharedInputs=True)\n",
        "    .combine(reducer2=ee.Reducer.stdDev(), sharedInputs=True),\n",
        "    geometry=seoraksan.geometry(),\n",
        "    scale=10,\n",
        "    maxPixels=1e9,\n",
        ")\n",
        "\n",
        "# 통계치 결과를 DataFrame으로 변환\n",
        "df_stats = pd.DataFrame(\n",
        "    [stats.getInfo()],\n",
        "    columns=[\"nd_min\", \"nd_mean\", \"nd_median\", \"nd_max\", \"nd_stdDev\"],\n",
        "    index=[\"Seoraksan\"],\n",
        ")\n",
        "df_stats.columns = [\"Min\", \"Mean\", \"Median\", \"Max\", \"StdDev\"]\n",
        "df_stats"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "3908fe65-c2a1-4042-b73c-407490b5624e",
      "metadata": {
        "id": "3908fe65-c2a1-4042-b73c-407490b5624e"
      },
      "outputs": [],
      "source": [
        "# DataFrame을 CSV 파일로 저장하기\n",
        "df_stats.to_csv('df_stats.csv', index=True)\n",
        "files.download('df_stats.csv')"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 3.2. 설악산생물권보전지역 산림손실 분석\n",
        "\n",
        "### 3.2.1. 설악산 생물권보전지역 경계 추출"
      ],
      "metadata": {
        "id": "J27Ipx9NPi54"
      },
      "id": "J27Ipx9NPi54"
    },
    {
      "cell_type": "code",
      "source": [
        "# 설악산 생물권보전지역 경계 추출\n",
        "wdpa = ee.FeatureCollection(\"WCMC/WDPA/current/polygons\")\n",
        "seorak = wdpa.filter(ee.Filter.eq('WDPAID', 5085))\n",
        "\n",
        "# 선택된 보호지역의 이름 확인\n",
        "wdpa_name = seorak.first().get('NAME').getInfo()\n",
        "print(\"Name:\", wdpa_name)"
      ],
      "metadata": {
        "id": "5CrfD7GMO53I"
      },
      "id": "5CrfD7GMO53I",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# 1: 외곽선의 픽셀 값, 2: 외곽선의 두께(픽셀 단위)\n",
        "wdpa_raster = ee.Image().paint(seorak, 1, 2)\n",
        "\n",
        "# 시각화 파라미터 설정\n",
        "wdpa_vis = {\n",
        "    'palette': ['#FFFF00'], # Yellow 색상 적용\n",
        "    'opacity': 1\n",
        "}"
      ],
      "metadata": {
        "id": "rhMBXWhbPzgt"
      },
      "id": "rhMBXWhbPzgt",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 3.2.2. 산림변화 데이터셋 사용"
      ],
      "metadata": {
        "id": "jdm6PFZPTDV-"
      },
      "id": "jdm6PFZPTDV-"
    },
    {
      "cell_type": "code",
      "source": [
        "# Hansen Global Forest Change v1.10 (2000-2022) 데이터셋\n",
        "gfc = ee.Image('UMD/hansen/global_forest_change_2022_v1_10')\n",
        "\n",
        "# 수관피복 시각화 파라미터\n",
        "tree_cover_vis = {\n",
        "    'bands': ['treecover2000'],\n",
        "    'min': 0,\n",
        "    'max': 100,\n",
        "    'palette': ['black', 'green']\n",
        "}\n",
        "\n",
        "# 산림손실 시각화 파라미터\n",
        "tree_loss_vis = {\n",
        "    'bands': ['lossyear'],\n",
        "    'min': 0,\n",
        "    'max': 22,\n",
        "    'palette': ['yellow', 'red']\n",
        "}"
      ],
      "metadata": {
        "id": "AzZt6IICS-DX"
      },
      "id": "AzZt6IICS-DX",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# gfc 레이어 추가\n",
        "m = geemap.Map(width=\"800px\", height=\"400px\")\n",
        "m.addLayer(gfc, tree_cover_vis, 'tree cover')\n",
        "m.addLayer(gfc, tree_loss_vis, 'tree loss year')\n",
        "m.addLayer(wdpa_raster, wdpa_vis, wdpa_name)\n",
        "m.add_layer(vworld_hybrid)\n",
        "m.centerObject(seorak, 10) # 지도의 중심 설정\n",
        "m # 지도 객체 출력"
      ],
      "metadata": {
        "id": "DKwE2fWsTp7w"
      },
      "id": "DKwE2fWsTp7w",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 3.2.3. 산림손실 면적 계산 및 시각화"
      ],
      "metadata": {
        "id": "4FtIZlBHUAlh"
      },
      "id": "4FtIZlBHUAlh"
    },
    {
      "cell_type": "code",
      "source": [
        "# 산림손실 이미지\n",
        "loss_image = gfc.select(['loss'])\n",
        "area_image = loss_image.multiply(ee.Image.pixelArea())\n",
        "\n",
        "# 설악산 BR 내 산림손실 면적 합산\n",
        "stats = area_image.reduceRegion(\n",
        "    reducer=ee.Reducer.sum(),\n",
        "    geometry=seorak.geometry(),\n",
        "    scale=30,\n",
        "    maxPixels=1e9,\n",
        ")\n",
        "\n",
        "# 산림손실 면적 출력\n",
        "loss_area = stats.get('loss').getInfo()\n",
        "print(f\"산림손실 면적: {loss_area:.2f}m²\")"
      ],
      "metadata": {
        "id": "KLCFfraMT9cj"
      },
      "id": "KLCFfraMT9cj",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# 산림손실연도 이미지\n",
        "lossyear_image = gfc.select(['lossyear'])\n",
        "\n",
        "# 2001-2022년 산림손실 면적 계산\n",
        "years = range(1, 23)\n",
        "loss_areas = []\n",
        "\n",
        "for year in years:\n",
        "    # 해당 연도의 산림손실 픽셀 선택\n",
        "    yearly_loss = lossyear_image.eq(year)\n",
        "\n",
        "    # 산림손실 면적 계산\n",
        "    area_image = yearly_loss.multiply(ee.Image.pixelArea())\n",
        "    stats = area_image.reduceRegion(\n",
        "        reducer=ee.Reducer.sum(),\n",
        "        geometry=seorak.geometry(),\n",
        "        scale=30,\n",
        "        maxPixels=1e9,\n",
        "    )\n",
        "\n",
        "    # 연도별 손실 면적을 리스트에 추가 (단위: 제곱미터)\n",
        "    loss_area = stats.getInfo()['lossyear']\n",
        "    loss_areas.append(loss_area)\n",
        "\n",
        "# 연도와 손실 면적을 데이터프레임으로 변환\n",
        "df = pd.DataFrame({\n",
        "    'Year': range(2001, 2023),\n",
        "    'Loss_Area_m2': loss_areas\n",
        "})\n",
        "\n",
        "df['Loss_Area_ha'] = df['Loss_Area_m2'] / 10000\n",
        "df.head(1)"
      ],
      "metadata": {
        "id": "vGeZJXmEUXyq"
      },
      "id": "vGeZJXmEUXyq",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# 연도별 산림손실 면적 막대 그래프 생성\n",
        "fig, ax = plt.subplots(figsize=(12, 6))\n",
        "colors = plt.cm.Reds(df['Loss_Area_ha'] / df['Loss_Area_ha'].max())\n",
        "bars = ax.bar(df['Year'], df['Loss_Area_ha'], color=colors)\n",
        "\n",
        "# 막대 위에 헥타르 단위로 면적 표시\n",
        "for bar, area_ha in zip(bars, df['Loss_Area_ha']):\n",
        "    yval = bar.get_height()\n",
        "    ax.text(bar.get_x() + bar.get_width()/2, yval, f'{area_ha:.0f}ha', va='bottom', ha='center', color='black')\n",
        "\n",
        "# 레이블 설정\n",
        "ax.set_xlabel('Year')\n",
        "ax.set_ylabel('Forest Loss Area (ha)')\n",
        "ax.set_title('Annual Forest Loss Area')\n",
        "\n",
        "# 차트를 이미지로 저장\n",
        "plt.savefig('forest_loss_area_colored.png', dpi=100)\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "up00If18U1Ua"
      },
      "id": "up00If18U1Ua",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "Fw7QR2ekLx-A"
      },
      "id": "Fw7QR2ekLx-A",
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.12.1"
    },
    "colab": {
      "provenance": [],
      "include_colab_link": true
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}