これまでSentinel-2のL2Aプロダクトという衛星画像データを中心に紹介してきました。このL2Aプロダクトというのは画像データ提供元のESAが大気補正などの補正処理をかけたデータです。
しかし、データ分析の内容によっては大気補正未処理のL1Cプロダクトを利用し、必要に応じて大気補正を行いたい場合があると思います。
今回はESAが公開しているデータ補正ツールSen2Corを紹介します。
Windows11(22H2), Sen2Cor 2.11
衛星データのプロダクトについて
人工衛星のデータは、検出器で検出された信号が格納された生データ(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に応じたものをダウンロードしてください。
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
で開始したローカル環境を終了し、元の環境に戻ります。
最後に
今回はSentinel-2のデータ補正ツールのSen2Corを紹介しました。
Sen2Corを利用すれば、L1Cの元データさえあればいつでもL2A相当のデータに補正することができ、さらにDEMを組み替えることで地形補正処理もカスタムすることができます。
処理自体はXMLファイルを修正し、後はバッチファイルを実行するだけなのでお手軽に利用できると思います。
参考
この記事を執筆するに当たり,以下の資料を参考にしました。
- Software User Manual: https://step.esa.int/thirdparties/sen2cor/2.11.0/docs/OMPC.TPZG.SUM.001 – i1r0 – Sen2Cor 2.11.00 Configuration and User Manual.pdf