Google_Earth_Engine_Sentinel-1_알고리즘_소개.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyORn5ANvjS3kA5UwluA29xw",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "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/Google_Earth_Engine_Sentinel-1_%EC%95%8C%EA%B3%A0%EB%A6%AC%EC%A6%98_%EC%86%8C%EA%B0%9C.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Google Earth Engine: Sentinel-1 알고리즘 소개\n",
        "\n",
        "이번 학습은 Google Earth Engine Guides에 소개된 [**Sentinel-1 Algorithms**](https://developers.google.com/earth-engine/guides/sentinel1)을 따라가 봅니다. 해당 전체 글을 번역하고 일부 내용을 보완했습니다. GEE에서 Sentinel-1 데이터를 다룰 때 한번 읽어보시면 좋겠습니다.\n",
        "\n",
        "## Sentinel-1 데이터\n",
        "\n",
        "Sentinel-1은 유럽연합(EU: European Union)이 자금을 지원하고 유럽우주국(ESA: European Space Agency)이 코페르니쿠스 프로그램(Copernicus Programme) 내에서 수행하는 우주 미션입니다.\n",
        "\n",
        "Sentinel-1은 다양한 편파와 해상도에서 C-밴드 합성 개구 레이더(SAR: Synthetic Aperture Radar) 이미지를 수집합니다.\n",
        "\n",
        "  * SAR 시스템은 지상이나 다른 대상에 전자파를 쏘고, 그 반사된 신호를 수신하는 방식으로 작동합니다. 이 반사된 신호를 분석하여 지구의 표면, 지형, 식생과 같은 여러 특성을 세밀하게 이미징할 수 있습니다. SAR 기술은 날씨 조건에 영향을 받지 않으며, 낮과 밤 모두 작동할 수 있기 때문에 위성 이미징과 지구 관측에 매우 유용하게 사용됩니다.\n",
        "\n",
        "  * 편파(Polarization)는 전자파, 특히 레이더에서 사용되는 전자파의 전기장 방향이 어떻게 진동하는지를 나타내는 특성입니다. SAR 같은 시스템에서 편파는 타겟의 물리적 및 지형적 특성을 파악하는 데 중요한 역할을 합니다. SAR 이미지는 다양한 편파(예: 수직 송신/수직 수신(VV), 수직 송신/수평 수신(VH) 등)를 사용하여 물체의 물리적 구조와 조직을 더 잘 이해할 수 있게 도와줍니다.\n",
        "  \n",
        "  * C-밴드의 주파수 범위는 약 4GHz에서 8GHz 사이이며, 이는 파장으로 환산하면 대략 3.75cm에서 7.5cm 사이입니다. 광학위성의 NIR 밴드 파장은 일반적으로 약 0.7µm에서 1.4µm 사이입니다. 이를 통해 볼 때, C-밴드의 파장은 NIR 밴드의 파장보다 대략 1000배에서 10000배 정도 긴 것을 알 수 있습니다.\n",
        "\n",
        "레이더 데이터는 보정된 정사영상 이미지를 얻기 위해 여러 전문 알고리즘이 필요합니다. 따라서 Earth Engine에서도 Sentinel-1 데이터의 전처리가 적용됩니다.\n",
        "\n",
        "Sentinel-1 데이터는 다양한 기기 구성, 해상도, 대역 조합으로 상승 궤도와 하강 궤도 모두에서 수집됩니다. 이러한 이질성 때문에 처리를 시작하기 전에 데이터를 동질의 부분집합으로 필터링하는 것이 보통 필요합니다."
      ],
      "metadata": {
        "id": "hA4GHFQr7juF"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Sentinel-1 데이터 전처리에 앞서 ipyleaflet을 0.18 버전으로 설치해 주겠습니다. 현재 Colab에 설치된 ipyleaflet 0.19 버전에 [오류](https://foss4g.tistory.com/2016)가 있어 다운그레이드하는 조치입니다."
      ],
      "metadata": {
        "id": "IzEYOClfyzle"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install -q ipyleaflet==0.18\n",
        "\n",
        "# 세션 재시작\n",
        "import os\n",
        "os.kill(os.getpid(), 9)"
      ],
      "metadata": {
        "id": "Z4GBOi5Q6co5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# ipyleaflet 버전 확인\n",
        "import ipyleaflet\n",
        "print(ipyleaflet.__version__)"
      ],
      "metadata": {
        "id": "ReCsS-QO-PS3"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Earth Engine을 인증 및 초기화합니다."
      ],
      "metadata": {
        "id": "fV3lhg5gzyn1"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "VxKLQbad6IQF"
      },
      "outputs": [],
      "source": [
        "import ee\n",
        "import geemap\n",
        "import geopandas as gpd\n",
        "import requests\n",
        "\n",
        "# Earth Engine 인증\n",
        "ee.Authenticate()\n",
        "\n",
        "# Earth Engine 초기화\n",
        "ee.Initialize(project='ee-foss4g')"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 메타데이터와 필터링\n",
        "\n",
        "Sentinel-1 데이터의 동질 부분집합을 만들기 위해서는 보통 메타데이터 속성을 사용하여 컬렉션을 필터링할 필요가 있습니다. 필터링에 사용되는 일반적인 메타데이터 필드에는 다음과 같은 속성들이 포함됩니다:\n",
        "\n",
        "1. `transmitterReceiverPolarisation`: ['VV'], ['HH'], ['VV', 'VH'], 또는 ['HH', 'HV']\n",
        "2. `instrumentMode`: 'IW' (Interferometric Wide Swath), 'EW' (Extra Wide Swath) 또는 'SM' (Strip Map).\n",
        "3. `orbitProperties_pass`: 'ASCENDING' 또는 'DESCENDING'\n",
        "4. `resolution_meters`: 10, 25 또는 40\n",
        "5. `resolution`: 'M' (medium) 또는 'H' (high).\n",
        "\n",
        "Sentinel-1 위성의 세 가지 주요 SAR 모드는 각각 다른 관측 폭과 해상도를 가집니다:\n",
        "  * IW 모드: 250km의 관측 폭을 가지며, 해상도는 5m x 20m입니다. 이 모드는 주로 육지 관측에 사용되며, 지형 관측에 사용되는 TOPSAR(Terrain Observation with Progressive Scans SAR) 기술을 활용하여 더욱 균일한 이미지 품질을 제공합니다​.\n",
        "  * EW 모드: 410km의 매우 넓은 관측 폭을 가지며, 해상도는 20m x 40m입니다. 주로 극지방과 해양 지역의 감시에 사용되며, 빙하 및 해상 오염 감시에 적합합니다​.\n",
        "  * SM 모드: 80km의 관측 폭과 5m x 5m의 가장 높은 해상도를 제공합니다. 이 모드는 매우 상세한 지도 제작이 필요한 경우나 특정 상황에서 요청될 때 사용됩니다​.\n",
        "\n",
        "Google Earth Engine에서 Sentinel-1 데이터는 원래의 해상도를 재가공하여 사용하기 쉬운 형태(10, 25, 40미터)로 제공합니다.\n",
        "\n",
        "상승 궤도(ASCENDING)는 위성이 남극에서 북극 방향으로 이동하는 궤도(주로 지구의 야간 시간대에 데이터를 수집)를 말하고, 하강 궤도(DESCENDING)는 위성이 북극에서 남극 방향으로 이동하는 궤도(주로 지구의 주간 시간대에 데이터를 수집)합니다.\n",
        "\n",
        "IW와 SM 모드는 'H' 해상도 등급에 속하고, EW 모드는 'M' 해상도 등급에 속한다고 볼 수 있습니다.\n",
        "\n",
        "아래 코드들은 `transmitterReceiverPolarisation`, `instrumentMode` 및 `orbitProperties_pass` 속성에 따라 Sentinel-1 컬렉션을 필터링하고, 여러 관측 조합에 대한 복합 이미지를 계산한 후 이러한 특성이 데이터에 어떻게 영향을 미치는지 보여주기 위해 지도에 표시합니다."
      ],
      "metadata": {
        "id": "rqQ5aS7AAFX8"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "[Sentinel-1 SAR GRD](https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S1_GRD) 이미지 컬렉션을 불러오고, 2023년 6월부터 9월까지의 관측 데이터로 필터링합니다."
      ],
      "metadata": {
        "id": "2P1wfZVs9jH_"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Sentinel-1 이미지 컬렉션을 불러오고, 2023년 6-9월까지의 관측 데이터로 필터링\n",
        "sentinel_1 = ee.ImageCollection('COPERNICUS/S1_GRD').filterDate(\n",
        "    '2023-06-01', '2023-10-01'\n",
        ")"
      ],
      "metadata": {
        "id": "GnlPEgACFhTg"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "`VV`와 `VH` 이중 편파로 촬영된 이미지, IW 모드로 수집된 이미지로 필터링합니다."
      ],
      "metadata": {
        "id": "5MKP2QVi96RN"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# 메타데이터 속성에 따라 Sentinel-1 컬렉션 필터링\n",
        "vv_vh_iw = (\n",
        "    sentinel_1\n",
        "    # VV와 VH 이중 편파로 촬영된 이미지로 필터링\n",
        "    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))\n",
        "    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))\n",
        "    # IW 모드로 수집된 이미지로 필터링\n",
        "    .filter(ee.Filter.eq('instrumentMode', 'IW'))\n",
        ")"
      ],
      "metadata": {
        "id": "M-Cp57CnH4qD"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "`orbitProperties_pass` 속성에 따라 상승 궤도와 하강 궤도 이미지를 별도의 컬렉션으로 분리합니다. 상승 궤도는 주로 야간에, 하강 궤도는 주로 주간에 데이터를 수집합니다. SAR는 라이다나 광학 위성과 달리 자체적으로 전자파를 발사하고 반사된 신호를 수신하기 때문에 낮과 밤의 구분이 일반적인 광학 위성 이미지만큼 큰 차이를 만들지는 않습니다.\n",
        "\n",
        "그러나 상승 궤도와 하강 궤도 데이터를 별도로 분석하는 것은 지형의 광원과 관측 각도 차이로 인한 지표면의 다른 물리적 특성을 정밀하게 이해하고 해석하기 위해서 중요할 수 있습니다."
      ],
      "metadata": {
        "id": "ia0_Plgp-PcW"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# 상승 궤도와 하강 궤도 이미지를 별도의 컬렉션으로 분리\n",
        "vv_vh_iw_asc = vv_vh_iw.filter(\n",
        "    ee.Filter.eq('orbitProperties_pass', 'ASCENDING')\n",
        ")\n",
        "vv_vh_iw_desc = vv_vh_iw.filter(\n",
        "    ee.Filter.eq('orbitProperties_pass', 'DESCENDING')\n",
        ")"
      ],
      "metadata": {
        "id": "eVwzn4PaIEOd"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "상승 궤도의 VH 평균, 하강 궤도의 VH 평균, 상승 및 하강 궤도 이미지 컬렉션을 결합한 VH와 VH 평균을 계산해 보겠습니다."
      ],
      "metadata": {
        "id": "Ic9gApRaB6Sp"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# 상승 궤도 VH 평균\n",
        "vh_iw_asc_mean = vv_vh_iw_asc.select('VH').mean()\n",
        "# 하강 궤도 VH 평균\n",
        "vh_iw_desc_mean = vv_vh_iw_desc.select('VH').mean()\n",
        "# 상승 및 하강 궤도 이미지 컬렉션을 결합한 VV 평균.\n",
        "vv_iw_asc_desc_mean = vv_vh_iw_asc.merge(vv_vh_iw_desc).select('VV').mean()\n",
        "# 상승 및 하강 궤도 이미지 컬렉션을 결합한 VH 평균.\n",
        "vh_iw_asc_desc_mean = vv_vh_iw_asc.merge(vv_vh_iw_desc).select('VH').mean()"
      ],
      "metadata": {
        "id": "4otdxEaBIXyM"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "특정 지역 대상으로 해당 레이어들을 표출해 보겠습니다. 진주시 센서스경계를 다운로드 받아 `EPSG:4326`으로 좌표계 변환 후 중심 좌표를 추출해 봅니다.\n",
        "  > `\"wb\"`는 'write binary'의 줄임말로, 파일을 쓰기 위해 바이너리 모드를 사용하겠다는 의미입니다. 파일은 크게 텍스트 모드와 바이너리 모드로 두 가지 방식으로 다룰 수 있으며, 바이너리 모드는 파일을 파이트 형태로, 즉 파일의 원래 형식과 정확한 내용을 변경 없이 저장하는 방식을 말합니다."
      ],
      "metadata": {
        "id": "8bYQccGSRYm0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# 진주시 센서스경계 다운로드\n",
        "url = \"https://github.com/osgeokr/GEE-PAM-Book/raw/main/JINJU.gpkg\"\n",
        "response = requests.get(url)\n",
        "with open(\"JINJU.gpkg\", \"wb\") as file:\n",
        "    file.write(response.content)\n",
        "\n",
        "# 진주시 센서스경계 읽기\n",
        "gdf = gpd.read_file(\"JINJU.gpkg\")\n",
        "\n",
        "# 좌표계를 WGS84 (EPSG:4326)로 변환\n",
        "gdf_wgs84 = gdf.to_crs(epsg=4326)\n",
        "gdf_wgs84.head()"
      ],
      "metadata": {
        "id": "cOhQqrVkXJ0V"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# gdf의 중심 좌표를 계산\n",
        "centroid = gdf_wgs84.geometry.unary_union.centroid\n",
        "longitude, latitude = centroid.x, centroid.y\n",
        "print(\"Longitude:\", longitude, \"Latitude:\", latitude)"
      ],
      "metadata": {
        "id": "mrq-VvQt_W4C"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "이제 Sentinel-1 이미지 레이어를 표출해 보겠습니다."
      ],
      "metadata": {
        "id": "qQycyaMPCf8p"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from ipyleaflet import TileLayer\n",
        "\n",
        "# Vworld 영상지도 객체\n",
        "Vworld_Satellite = TileLayer(\n",
        "    url='https://xdworld.vworld.kr/2d/Satellite/service/{z}/{x}/{y}.jpeg',\n",
        "    name='Vworld Satellite',\n",
        "    attribution='Vworld',\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",
        "# 다양한 관측의 시간 평균을 표시 및 비교\n",
        "m = geemap.Map(width=\"800px\", height=\"500px\")\n",
        "m.add_layer(Vworld_Satellite)\n",
        "m.add_layer(vv_iw_asc_desc_mean, {'min': -12, 'max': -4}, 'vv_asc_desc_mean')\n",
        "m.add_layer(vh_iw_asc_desc_mean, {'min': -18, 'max': -10}, 'vh_asc_desc_mean')\n",
        "m.add_layer(vh_iw_asc_mean, {'min': -18, 'max': -10}, 'vh_asc_mean')\n",
        "m.add_layer(vh_iw_desc_mean, {'min': -18, 'max': -10}, 'vh_desc_mean')\n",
        "m.add_layer(Vworld_Hybrid)\n",
        "\n",
        "m.set_center(longitude, latitude, 15)  # 경상남도 진주시\n",
        "m"
      ],
      "metadata": {
        "id": "krihpMwUUSw-"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Sentinel-1 전처리\n",
        "\n",
        "Earth Engine의 `COPERNICUS/S1_GRD` Sentinel-1 `ImageCollection`에 포함된 이미지는 Level-1 Ground Range Detected(GRD), 레벨-1 지상 범위 감지 형식으로 처리되어 후방산란 계수(σ°)를 데시벨(dB) 단위로 나타냅니다.\n",
        "  * GRD는 데이터가 지상의 실제 거리(지상 범위)를 반영하도록 조정된 것을 의미하며, 후방산란 계수는 지면 단위 면적당 레이더 파동이 얼마나 반사되는지를 나타내는 지표입니다.\n",
        "  \n",
        "후방산란 계수는 매우 다양할 수 있어, 이를 데시벨(dB) 단위로 변환하여 사용합니다. 데시벨로의 변환 공식은 10*log10σ°입니다.\n",
        "  * 데시벨 값이 0보다 작으면 지형이 SAR 센서로부터 멀어지는 방향으로 파동을 산란시키는 것을 의미하며, 0보다 크면 센서 쪽으로 파동을 산란시키는 것을 의미합니다. 이러한 산란 현상은 지형의 형태(기하학적 특성)와 재질(전자기적 특성)에 따라 달라집니다.\n",
        "\n",
        "  * 데시벨(dB) 단위는 신호의 크기나 강도를 로그 스케일로 변환하여, 매우 넓은 값의 범위를 쉽게 비교하고 분석할 수 있게 해줍니다. 이는 특히 SAR 데이터에서 지형의 미세한 차이를 효과적으로 측정하고 비교하는 데 유용합니다.\n",
        "\n",
        "간단히 말해, 레이더가 지형에서 얼마나 많은 신호를 받아내는지를 정량화한 것으로, 지형의 세부 사항과 그 특성을 이해하는 데 중요한 정보를 제공합니다.\n",
        "\n",
        "Earth Engine은 각 픽셀에서 후방산란 계수를 도출하기 위해 Sentinel-1 Toolbox에서 구현된 다음과 같은 전처리 단계를 사용합니다:\n",
        "\n",
        "1. **Apply orbit file**(궤도 파일 적용): 복원된 궤도 파일(또는 복원된 파일이 사용 불가능한 경우 정밀 궤도 파일)로 궤도 메타데이터를 업데이트합니다.\n",
        "\n",
        "1. **GRD border noise removal**(GRD 테두리 잡음 제거): 장면 가장자리의 저강도 잡음과 유효하지 않은 데이터를 제거합니다. (2018년 1월 12일 기준)\n",
        "\n",
        "1. **Thermal noise removal**(열 잡음 제거): 다중-스워스 획득 모드에서 서브-스워스 간 불연속성을 줄이기 위해 서브-스워스의 가산 잡음을 제거합니다. (이 작업은 2015년 7월 이전에 생산된 이미지에는 적용할 수 없습니다)\n",
        "\n",
        "1. **Application of radiometric calibration values**(방사 보정 값 적용): GRD 메타데이터의 센서 보정 매개변수를 사용하여 후방산란 강도를 계산합니다.\n",
        "\n",
        "1. **Terrain correction**(지형 보정) 또는 Orthorectification(정사 보정): 지형을 고려하지 않는 지상 범위 기하학에서 SRTM 30미터 DEM 또는 고위도(60° 이상 또는 -60° 이하)의 경우 ASTER DEM을 사용하여 σ°로 데이터를 변환합니다."
      ],
      "metadata": {
        "id": "ihjc49EOPPW7"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "##데이터셋 주의사항\n",
        "\n",
        "- 산지경사에서 발생하는 오류 때문에 Radiometric Terrain Flattening(RTF)는 적용되지 않습니다. RTF는 SAR 이미지의 레이더 반사율 데이터를 지형의 영향으로부터 보정하는 기술입니다.\n",
        "\n",
        "- 단위 없는 후방산란 계수는 위에서 설명한 대로 dB로 변환됩니다.\n",
        "\n",
        "- Sentinel-1 SLC 데이터는 현재 Earth Engine에서 취급할 수 없습니다. 이는 Earth Engine이 피라미드화하는 동안 위상 정보의 손실 없이 복소수(complex number) 값을 평균화할 수 없기 때문에 복소수 값을 포함하는 이미지를 지원하지 않기 때문입니다. 복소수 데이터는 각 픽셀에서 수신된 레이더 신호의 진폭(amplitude)과 위상(phase) 정보를 모두 포함합니다.\n",
        "\n",
        "- GRD SM 자원은 S1 Toolbox의 테두리 잡음 제거 작업에서 computeNoiseScalingFactor() 함수가 SM 모드를 지원하지 않기 때문에 취급되지 않습니다."
      ],
      "metadata": {
        "id": "XAjvgPM9XZ_9"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "IsFGNW6EPTsC"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}