【衛星データ】Sen2Corによる衛星データのデータ補正について【データ補正】

これまでSentinel-2のL2Aプロダクトという衛星画像データを中心に紹介してきました。このL2Aプロダクトというのは画像データ提供元のESAが大気補正などの補正処理をかけたデータです。

しかし、データ分析の内容によっては大気補正未処理のL1Cプロダクトを利用し、必要に応じて大気補正を行いたい場合があると思います。

今回はESAが公開しているデータ補正ツールSen2Corを紹介します。

動作検証済み環境

Windows11(22H2), Sen2Cor 2.11

Google Earth Engineを用いた衛星データ解析入門 ~タイムラプスを作成し地球の環境変動を観察~

Pythonを使った衛星データ解析の初心者向け技術書を販売中！

¥2,500 → ¥800 今なら約70%OFF!!

Google Earth Engineを使えば簡単に
自宅にいながら地球環境の変動を観察できます

技術書ページへ

衛星データのプロダクトについて

人工衛星のデータは、検出器で検出された信号が格納された生データ（Rawデータ）から様々な補正処理が施され、その処理レベルに応じたプロダクトが提供されます。一般的にはL1やL2といったレベルのプロダクトが利用者には提供されます。

以下、Sentinel-2の例を紹介します。

L0 (Level 0): 生データ L0データは、Sentinel-2の検出器で検出された信号が格納された最も基本的なデータです。これは、光子が検出器で受け取られた際に生成される電荷の情報を含んでいます。L0データは、画像としてはまだ解釈できない生の信号データであり、これを基に次の処理ステップに進みます。通常L0のデータは人工衛星を管理している機関や会社が保有し、あまり公開されることがありません。
L1 (Level 1): 放射量補正データ L1データは、L0データに対して放射量補正が行われたデータです。L1データはさらに以下の3つのサブレベルに分類されます。
- L1A: Level-1A – 人工衛星がスキャンした帯状のL0データの解凍、圧縮処理が行われています。
- L1B: Level-1B – ラジオメトリック補正が実施されたデータです。この処理により、L1Aの輝度値が大気上端の放射輝度値として補正されます。また、幾何補正以外にも、欠損画素の補間や暗電流によるノイズの補正などが含まれています。これにより、より正確で解釈しやすい画像データが得られます。
- L1C: Level-1C – 輝度値を反射率に変換し、オルソ補正が実施されたデータです。デジタル標高モデル（DEM）を使用して画像データを地図投影することで位置情報が補正されます。また、反射率は大気上端の反射率（TOA: Top of Atmosphere）として提供されます。
L2 (Level 2): 大気補正データ L2A: Level-2A – L1Cに対して大気補正が行われたデータです。大気補正により、大気中のガス（水蒸気など）やエアロゾルによる影響が除去され、併せて地形高度補正も行い影などの影響を除去しています。反射率としては地表面反射率（BOA: Bottom of Atmosphere）が提供されます。

このうち、Sentinel-2ではL1CとL2Aが提供されています。

主なデータ補正処理について

衛星画像データには撮影時の衛星の位置や大気状態により誤差やゆがみが含まれています。そのため、画像データを扱う際には複数のデータ補正処理が行われます。

代表的なデータ補正は下記があります。

ラジオメトリック補正：センサーの応答特性を考慮して、画像データの放射測定値を補正する方法です。
オルソ補正：地形の影響を考慮して、画像データを地図上に正確に配置する方法です。
大気補正：大気の影響を考慮して、画像データの放射測定値を補正する方法です。

Sen2Corについて

Sen2Corは、ESAが提供するSentinel-2衛星データに対する大気補正ツールです。このツールは、主にSentinel-2のL1CデータをL2Aデータに変換するために使用されます。

Sen2Corの主な処理は以下の通りです：

大気補正: 大気中のガスやエアロゾルによる画像データへの影響を除去し、地表面の反射率を正確に計算します。
雲マスク: 雲の存在を検出し、その影響を考慮した画像解析が可能になります。
水蒸気補正: 画像データ中の水蒸気の影響を除去します。
地表高度補正: 地形による影響を考慮して、地表面の反射率を補正します。

また、ユーザ側ではSen2CorのL2A_GIPP.xmlというXMLファイルのパラメータを編集することでカスタマイズした補正を行うことができます。

大気モデルやエアロゾルモデルの選択
雲検出や雪氷検出の閾値
シーン分類の方法やクラスの定義
標高データや土地被覆データの指定

Sen2Corを使ったデータ補正方法

Sen2Corのダウンロード

Sen2CorはESAの以下のサイトから最新版がダウンロードできます。最新版の2.11はWindowsとLinuxに対応しているので使用するOSに応じたものをダウンロードしてください。

Sen2Cor v2.11

Windows版ではSen2Cor-02.11.00-win64.zipというZipファイルがダウンロードされるので解凍します。フォルダ構成は下記の通りです。L2A_Process.batというバッチファイルを実行することで補正処理が行われます。

Sen2Cor-02.11.00-win64
|_ LICENSE.txt
|_ L2A_Process.bat
|_ share
|_ Lib
|_ DLLs
|_ copyright
|_ bin

L1Cプロダクトのダウンロード

今回は衛星データはESAのCopernicus Open Access Hubからダウンロードします。アカウント登録が必要なため、下記サイトから必要事項を記入してアカウントを登録しましょう。

https://scihub.copernicus.eu/dhus/#/self-registration

アカウント登録が完了したら、ブラウザ上でもダウンロード出来ますが、今回は作業効率化のため以下のコードを実行してダウンロードを行います。コード内では比較用にL2Aのダウンロードも行っています。

import os
import zipfile
from datetime import date

import geopandas as gpd
from sentinelsat import SentinelAPI, geojson_to_wkt, read_geojson
from shapely.geometry import Polygon

# 座標リストを使ってポリゴンを作成し、GeoJSONファイルに保存する関数
def create_polygon(coordinates):
    polygon = Polygon(coordinates)
    gdf = gpd.GeoDataFrame(geometry=[polygon], crs="EPSG:4326")
    gdf.to_file("search_area.geojson", driver='GeoJSON')
    return read_geojson('search_area.geojson')

# Sentinel APIを使って衛星画像を検索する関数
def search_sentinel_images(api, footprint, start_date, end_date, platformname, processinglevel, cloudcover):
    products = api.query(footprint,
                         date=(start_date, end_date),
                         platformname=platformname,
                         processinglevel=processinglevel,
                         cloudcoverpercentage=cloudcover)
    return products

# 検索結果から画像をダウンロードする関数
def download_best_image(api, products, download_folder):
    if len(products) == 0:
        print("No products found within the specified parameters.")
        return None

    best_product_id = sorted(products, key=lambda p: products[p]['cloudcoverpercentage'])[0]
    downloaded_product = api.download(best_product_id, directory_path=download_folder)
    return downloaded_product

# ダウンロードしたzipファイルを解凍する関数
def unzip_file(file_path, output_folder):
    with zipfile.ZipFile(file_path, 'r') as zip_ref:
        zip_ref.extractall(output_folder)

# メイン処理(衛星画像のダウンロード)
if __name__ == '__main__':
    user = 'your_user_name' # OpenAccessHubのユーザ名
    password = 'your_password' # OpenAccessHubのパスワード
    api = SentinelAPI(user, password, '<https://scihub.copernicus.eu/dhus>')

    # 目的地のポリゴンの座標リストを定義
    coordinates = [[130.43061, 30.39225],
                   [130.43061, 30.28647],
                   [130.61729, 30.28647],
                   [130.61729, 30.39225]]

    # 座標リストを使ってポリゴンを作成し、そのfootprintを取得
    footprint = geojson_to_wkt(create_polygon(coordinates))

    start_date = '20230307'
    end_date = '20230308'
    download_folder = 'your_download_folder'
    platformname = 'Sentinel-2'
    cloudcover = (0, 30)

    # L1C衛星画像を検索し、ダウンロードして解凍
    products_l1c = search_sentinel_images(api, footprint, start_date, end_date, platformname, 'Level-1C', cloudcover)
    downloaded_product_l1c = download_best_image(api, products_l1c, download_folder)
    if downloaded_product_l1c is not None:
        unzip_file(downloaded_product_l1c['path'], download_folder)

    # L2A衛星画像を検索し、ダウンロードして解凍
    products_l2a = search_sentinel_images(api, footprint, start_date, end_date, platformname, 'Level-2A', cloudcover)
    downloaded_product_l2a = download_best_image(api, products_l2a, download_folder)
    if downloaded_product_l2a is not None:
        unzip_file(downloaded_product_l2a['path'], download_folder)

コードを実行した場合、S2A_MSIL1C_hogehoge.SAFE、S2A_MSIL2A_hogehoge.SAFEというような名前のフォルダが作成されていれば成功です。

XMLファイルの修正

Sen2Corは補正を行う際にL2A_GIPP.xmlというXMLファイルに記載されたパラメータを参照します。L2A_GIPP.xmlはC:/Users/<ユーザ名>/Documents/sen2cor/2.11/cfgに配置する必要があります。元ファイルはSen2Cor-02.11.00-win64/Lib/site-packages/sen2cor/cfgにあるのでコピーして配置しましょう。

また、地形補正はデフォルトのパラメータでは未設定なので、あらかじめ設定する必要があります。<DEM_Directory>と<DEM_Reference>を以下のように修正しましょう。

<DEM_Directory>**dem/CopernicusDEM90_DGED**</DEM_Directory>
<!-- should be either a directory in the sen2cor home folder or 'NONE'. If NONE, no DEM will be used -->
<DEM_Reference>**<https://prism-dem-open.copernicus.eu/pd-desk-open-access/prismDownload/COP-DEM_GLO-90-DGED__2021_1/**></DEM_Reference>
<!-- DEM_Reference><http://srtm.csi.cgiar.org/wp-content/uploads/files/srtm_5x5/TIFF/></DEM_Reference
     disable / enable the upper two rows if you want to use an SRTM DEM
     The SRTM DEM will then be downloaded from this reference, if no local DEM is available
     if you use Planet DEM you can optionally add the local path instead,
     which then will be inserted in the datastrip metadata -->

<DEM_Directory>では地形補正用のDEMデータを配置するディレクトリを入力します。基本的にはC:/Users/<ユーザ名>/Documents/sen2cor/2.11/dem/CopernicusDEM90_DGEDのような構成にします。次に<DEM_Reference>にDEMデータを取得するリファレンスを入力します。L1CからL2Aへの補正には通常GLO-90というモデルが利用されているようです。そのため、取得用のURLhttps://prism-dem-open.copernicus.eu/pd-desk-open-access/prismDownload/COP-DEM_**GLO-90**-DGED__2021_1/を入力します。他にはGLO-30というモデルも用意されているようなので、COP-DEM_**GLO-30**-DGED__2021_1に修正すればよさそうです。

Sen2Corを実行する

次にPowerShellを開き、以下のコマンドでSen2Cor-02.11.00-win64にディレクトリを移動します。

$ cd <保存したフォルダ名>/Sen2Cor-02.11.00-win64

あとは以下のコマンドを実行するだけです。引数には先ほどダウンロードしたS2A_MSIL1C_hogehoge.SAFEを指定します。（ $マークは無視してください）

$ ./L2A_Process.bat <S2A_MSIL1Cから始まる取得したL1Cデータが格納されたフォルダ>

実行結果

途中でエラーメッセージが出ずに、上記の様なProgress[%]: 100.00 : Application terminated successfully. というメッセージがPowerShell上で表示されていれば成功です！S2A_MSIL2A_hogehoge_実行日.SAFEという名前のフォルダが作成されていれば成功です！

では、実際に処理が行われているかを画像データで確認します。左からL1C、L2A、Sen2Corで補正したL1Cになります。L1Cと比較すると、画像データの全体的な青白さが抑えられ、さらに山による陰影も補正されているのがわかります。

コードの解説

上記で実行したバッチファイルL2A_Process.batについて解説をしていきます。

L2A_Process.batを開くと中は下記のように記述されています。

@echo off

setlocal

set CURRENT_SCRIPT_DIR=%~dp0

set SYSPATH=%SystemRoot%\\system32;%SystemRoot%;%SystemRoot%\\System32\\Wbem;
set PATH=%CURRENT_SCRIPT_DIR%\\bin;%SYSPATH%

set GDAL_DATA=%CURRENT_SCRIPT_DIR%\\share\\data
set GEOTIFF_CSV=%CURRENT_SCRIPT_DIR%\\share\\epsg_csv
:: unset some variables that allows to have an isolate python environment

set PYTHONPATH=
set PYTHONHOME=
set PYTHONEXECUTABLE=
set PYTHONUSERBASE=

set GDAL_DRIVER_PATH=disable

set SEN2COR_HOME=%USERPROFILE%\\Documents\\sen2cor\\2.11
set SEN2COR_BIN=%CURRENT_SCRIPT_DIR%\\Lib\\site-packages\\sen2cor

::echo "setting SEN2COR_HOME is %SEN2COR_HOME%"

if not exist "%SEN2COR_HOME%" ( 
mkdir %SEN2COR_HOME%
)
if not exist "%SEN2COR_HOME%\\cfg" mkdir %SEN2COR_HOME%\\cfg

if not exist "%SEN2COR_HOME%\\cfg\\L2A_GIPP.xml" (
 xcopy /F "%SEN2COR_BIN%\\cfg\\L2A_GIPP.xml" "%SEN2COR_HOME%\\cfg\\"
)

%CURRENT_SCRIPT_DIR%\\bin\\python.exe -s %SEN2COR_BIN%\\L2A_Process.py %*

endlocal

このスクリプトの処理内容は環境変数の設定とL2A_Process.pyの実行を行い、L1Cデータの補正処理を行っています。

@echo off：コマンドプロンプトで表示されるコマンドの出力をオフにします。
setlocal：環境変数の変更をローカル化します。このバッチファイルの実行が終了すると元に戻ります。
スクリプトがあるディレクトリをCURRENT_SCRIPT_DIRに設定します。
システムパスをSYSPATHに設定し、その後PATH変数にbinディレクトリを追加します。
GDAL_DATAおよびGEOTIFF_CSV環境変数を設定します。これらは、地理空間データフォーマットをサポートするためのライブラリであるGDALに関連しています。
Python環境変数（PYTHONPATH、PYTHONHOME、PYTHONEXECUTABLE、PYTHONUSERBASE）をリセット（未設定）します。
GDAL_DRIVER_PATHを無効に設定します。
SEN2COR_HOMEとSEN2COR_BIN環境変数を設定します。これらは、Sen2Corのホームディレクトリとバイナリディレクトリを指します。
必要に応じて、Sen2Corのホームディレクトリと構成ディレクトリを作成します。
L2A_GIPP.xmlファイルが存在しない場合、SEN2COR_BINディレクトリからSEN2COR_HOMEディレクトリにコピーします。
L2A_Process.py Pythonスクリプトを実行して、Sentinel-2データを処理します。Sen2Corにあらかじめ格納されているPythonモジュールで、他にESA側で用意されたPythonモジュールの呼び出しを行うなどのメインの処理が管理されたモジュールです。ユーザ側では編集することは基本的にないと思われます。
endlocal：setlocalで開始したローカル環境を終了し、元の環境に戻ります。