Sentinel-2の衛星画像をコマンドライン（GDAL）でラスタタイルにする。（QGISなし）

はじめに

Sentinel-2の観測データをダウンロードして、トゥルーカラーの衛星画像タイルを作る手順をメモします。今回は、なるべく手間をかけたくないので、ほとんどのプロセスをコマンドラインで処理した手順です。
QGISは使いません。

Sentinel-2の可視光バンドの空間解像度は10メートルなので、ズームレベルは15までとします（最小ズームは9にしておく）。

（お断り）
複数シーンを使うときは、投影法を合わせてモザイク処理などが必要になりますが、今回は単一のシーンでの手順です。

環境

• Docker環境 (unvt/nanbanイメージ ubuntuベースです。)
  • GDAL 3.4.1, released 2021/12/27

手順

Step 1: データのダウンロード

Copernicus Open Access HubからSentinel-2データをダウンロードしてきます。ユーザー登録が必要です。

どこでもいいのですが、プロダクトタイプ S2MSI2A （処理レベル2A）のものをダウンロードします。ZIPファイルで1GB程度かと思います。ファイル名が長いですが、ZIPファイルがダウンロードできます。作業するフォルダにsrcというサブフォルダを作って入れておくことにします。

↓ダウンロードしてくるページ

↓ダウンロードしてきたデータ

初めての方は以下を参考にしてもよいと思います。

Step 2: Docker起動

gdal関係ツールが使えればなんでもいいのですが、bashコマンドも使いたかったのでDockerを使うことにします。unvt/nanban イメージを使います。

docker run -it --rm -v ${PWD}:/data unvt/nanban
cd /data

Step 3: ファイルを解凍して必要なファイルを取り出す

とりあえずZIPファイルを解凍します。コマンドラインで処理しますが、同じsrcフォルダに解凍することにします。

unzip src/*.zip -d src

解凍が終わると、.SAFEフォルダの中の、GRANULEフォルダの中の、プロダクトフォルダの中のIMG_DATAフォルダがありまして（下図）、、

その中の、R10mのフォルダの中に、必要なバンドのデータが入っています（下図）。今回は、トゥルーカラー合成なので、バンド２，３，４を使います。

01_bandというフォルダを作って、必要なファイルをコピーしてきます。

mkdir 01_band
cp src/*.SAFE/GRANULE/*/IMG_DATA/R10m/*B02_10m.jp2 01_band/band02.jp2
cp src/*.SAFE/GRANULE/*/IMG_DATA/R10m/*B03_10m.jp2 01_band/band03.jp2
cp src/*.SAFE/GRANULE/*/IMG_DATA/R10m/*B04_10m.jp2 01_band/band04.jp2

なお、ここまで出来たら、gdalinfoを使って、当該シーンの投影法を確認しておきます。
（今回は単一シーンなので、タイル化の前に投影変換をしないため、オリジナルの投影を確認しておきますが、複数のシーンを扱ったり、する場合はこのタイミングで投影変換をしてもいいと思います。）

gdalinfo 01_band/band02.jp2

下の方にある、UTM のEPSG番号を確認しておきましょう。ここでは、32654です。

ほかのバンドはいらないので、解凍したフォルダを削除しておきます。

rm -rf src/*.SAFE

Step 4: jp2をtiffに変換 (gdal_translate)

jp2は16bitのデータなので、RGB合成の前に計算をします。しかし gdal_calc.py はjp2を扱えないので、tiffに変換しておきます。

for f in 01_band/*.jp2; do gdal_translate ${f} 01_band/`basename ${f} .jp2`.tif; done;

これが終わると以下の感じになります。

Step 5: ストレッチ、および 16bit を　8bitに計算 (gdal_calc.py)

衛星画像の処理レベルは2Aなので、各バンドの値は反射率になっています。これは0～1の値をとりますがデータとしては0～1の反射率を0～10,000で表現しているようです。

私の経験上、反射率は0～0.35くらいを線形でストレッチ（0～250くらい）してあげればそこそこの画像になると思います。
ですので、gdal_calc.pyをつかって、ストレッチをするのと同時にピクセルの深度を16bitから8bitに落とします。

2_8bitというフォルダを作って計算結果を出力します。

mkdir 2_8bit
for f in 01_band/*.tif; do gdal_calc.py -A ${f} --calc="(A<=0)*0 + (A>0)*(A<=3500)*A/14 + (A>3500)*250" --type=Byte --outfile 2_8bit/`basename ${f}`; done;

（式の説明）
インプットラスタの0～3500をアウトプットラスタの0～250にするので、分布する値を14で割ればよいです。0より小さいものは0にして、3500より大きいものは250にします。

Step 6: RGB合成 (gdal_merge.py)

3つのバンドの画像ができたので、今度はRGB画像に合成します。3_mergeというフォルダを作って出力します。アウトプットファイルはrgb.tifという名前にしておきます。経験上、バンドの順は4-3-2（赤、緑、青）の順がよいみたいです。

mkdir 3_merge
gdal_merge.py -o 3_merge/rgb.tif -co PHOTOMETRIC=RGB -separate 2_8bit/band04.tif  2_8bit/band03.tif 2_8bit/band02.tif

Step 7: タイル化 (gdal2tiles.py)

あとは、ラスタタイルにします。ズームレベルは9～15にします。ソースファイルの投影法を指定するのを忘れないようにしましょう。ここでは、前に確認したEPSG:32654になります。

gdal2tiles.py -s "EPSG:32654" -z "9-15" --xyz 3_merge/rgb.tif

しばらく待つとタイルができます。入力したTIFFファイルの名前と同じ名前のフォルダができて、その中にラスタタイルが生成されています。このコマンドを使うと、フォルダの中にOpenLayersやLeafletなどのプレビューがついてくるので、そのままGitHubページなどにアップしても成果を閲覧できます。

なお、Sentinel-2 のシーンは大体100km四方ですが、ラスタタイルにして700MBくらいでした（ZL9-15）。1つのシーンのSentinel画像であれば、一つのGitHubレポジトリにアップできると思います。（中間ファイルは.gitignoreでよけておく必要あり。）

※ なお、mbutilでmbtilesにまとめても、データサイズはほとんど変わりませんでした。（むしろちょっと大きくなった。）

まとめ

最初のダウンロードのところは違いますが、後はコマンドラインの処理でラスタタイルを作りました。
地域、季節によっては、反射率0-0.35のストレッチが微妙かもしれないので、調整が必要になるかもしれません。（大気補正をしてある地表反射率なので、大気上端反射率の時ほど調整はいらないと思います。）

謝辞

コペルニクスプロジェクトに感謝します。新鮮な10メートル解像度の衛星画像がオープンに使えるのは素晴らしいです。また、GDALツールにも感謝します。便利です。

参考文献