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](https://scihub.copernicus.eu/dhus/#/home)からSentinel-2データをダウンロードしてきます。ユーザー登録が必要です。

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

↓ダウンロードしてくるページ
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/bcf16aaa-c821-ccb2-9bd9-324f183fb017.png)

↓ダウンロードしてきたデータ
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/2bdf13a9-6b59-8989-674f-dae067cd300b.png)


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

https://qiita.com/T-ubu/items/e40e49b3321dfc720b78

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

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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/d4fe332a-1c65-c5f7-584d-0dd0420fc4a6.png)

## Step 3: ファイルを解凍して必要なファイルを取り出す
とりあえずZIPファイルを解凍します。コマンドラインで処理しますが、同じsrcフォルダに解凍することにします。
```bash
unzip src/*.zip -d src
```

解凍が終わると、.SAFEフォルダの中の、GRANULEフォルダの中の、プロダクトフォルダの中のIMG_DATAフォルダがありまして（下図）、、
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/7d174baa-4703-731d-3b11-d0b7c471433a.png)
その中の、R10mのフォルダの中に、必要なバンドのデータが入っています（下図）。今回は、トゥルーカラー合成なので、バンド２，３，４を使います。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/a064b51b-92d1-45fc-5c84-762073caa313.png)

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

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

```bash
gdalinfo 01_band/band02.jp2
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/46d1b26b-3a85-38c6-1792-c5238178f488.png)
下の方にある、UTM のEPSG番号を確認しておきましょう。ここでは、32654です。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/0618de8f-ca58-0b27-8d7a-30054432d0ef.png)

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

```bash
rm -rf src/*.SAFE
```

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

```bash
for f in 01_band/*.jp2; do gdal_translate ${f} 01_band/`basename ${f} .jp2`.tif; done;
```
これが終わると以下の感じになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/6e004c6a-96ae-b469-7d9c-09554b50e006.png)



## 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というフォルダを作って計算結果を出力します。

```bash
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にします。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/9af3cdb5-13a7-00e2-a140-368d07f159f8.png)

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

```bash
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
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/d7a504d6-fe8d-a904-ae31-87debe26fb2c.png)

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

```bash
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ツールにも感謝します。便利です。

# 参考文献
https://scihub.copernicus.eu/dhus/#/home

https://gdal.org/

https://qiita.com/T-ubu/items/e40e49b3321dfc720b78