はじめに

昨日早起きしてLandsat-8（単シーン）のラスタタイルを作ったので、今日はSentinel-2の4シーンでラスタタイルを作ってみます。QGISを使って作業します。

2021年11月末～12月上旬のデータ
対象地域は関東地方（Sentinel 2B で4シーン）
ZLは6～13
True color （B4,B3,B2）

二度寝したいので、作業時間は1:40AM～3:40AMの2時間を目標にしたいですが、Qiitaでメモをするともう少しかかるかもしれません。

##　環境

QGIS version 3.16.6
Google Chrome (データのダウンロードのためだけ)
Windows PC

手順

Step1：シーンの選択、データダウンロード

昨日（こちらの記事）と同様にEarthExplorerからダウンロードします。USGSで扱っているSentinel-2についての説明はこちら→USGS EROS Archive - Sentinel-2。EarthExplorerから雲がなくてなるべく新しい画像を選びます。

Sentinel-2の場合、2Aと2Bでシーン範囲が少し違いますが、関東の周辺は2Bだと、だいたい4シーンでカバーされています。幸いに西側も東側もUTM54ゾーンに入っています。今回は簡単のため、2Bのシーンを中心に選び、雲が少ない日のものにしました。西側は2021年12月3日、東側は2021年11月28日の画像です。処理レベルはL1C（DEMを使ってジオレファレンス済、TOA反射率）とのことです。

西側	東側
T54SUF	T54SVF
T54SUE	T54SVE

ダウンロードしたzipファイルを解凍すると、画像はGRANULEというフォルダのIMG_DATAというフォルダに入っていました。 USGS EROS Archive - Sentinel-2のページでバンドの波長を確認して、バンド4,3,2をつかってTrue Colorの画像を作ることにします。

Step2：QGSIでの読み込み

4つのシーンでそれぞれ3つの画像を読み込みます（B4,B3,B2）。一つの画像110MBと少しなのでそんなに大きくないのですが、バンドごとにマージしてからRGBにするか、先に各シーンでRGBを作るかデータを見て考えました。4シーンのBand4のヒストグラム見比べてみましたが結論がでないので、とりあえずシーンごとでRGB画像を作ってから調整することにします。（どのシーンも0－2000の間にだいたい分布しているんですが、地理的条件も違うのでなんともいえません。）

Step3: RGB画像作成（シーンごと）

昨日Landsatでやったのと同様に、RasterメニューからVirtual Rasterを作りました。4つのシーンについてそれぞれトゥルーカラーのRGB合成画像を作ります。
（インプットファイルを別バンドにするというチェックを忘れたのでここはやり直しました。気をつけましょう。）

レイヤを選ぶときは読み込んでいるリストがあるので、その中から対象シーンのRGBを選びます。念のため、B4,B3,B2の順番にしています。（インプットレイヤで3レイヤしか指定していないのに4レイヤ選んだことになってしまうのは謎ですが実害はないのでそのままにしておきます。）

すると合成された画像が4つできます。RGBの配分がどうもおかしそうなのでこの後のステップで調整します。

Step4: 色の調整

4つの合成画像について色を調整していきます。まずレイヤのプロパティから各バンドのヒストグラムをチェックしました。すると、1）なぜか3つのバンドが、バンド1、バンド3、バンド4に入っていることと、2）値はだいたい2000以下に分布していることがわかります。（Landsatのときは、バンド2、バンド3、バンド4でした。なんでこうなるのかは謎ですが気にせず次に進みます。）

上記の観察を踏まえて、レイヤのプロパティのSymbologyから色の調整をします。RGBは合成画像中のBand1,Bnd3,Band4を割り当てます。元の衛星画像のバンドと合成画像中のバンドの関係をしっかり追わないと混乱しそうですね。そして、色のストレッチについて、4シーンあって面倒なので、まずはどのバンドもminが0、maxが2000でいってみます。（これが以外とOKでした。観測日が近かったからでしょうか。）

結果はこんな感じです。4シーンあって、縦でならんでいる2つは同日観測なので同じストレッチでも問題ないだろうと予想していたのですが、横の接合もそんなに悪くなさそうです。

真ん中にある橋が1ピクセルくらいずれているのですが、ここにシーンの境界（南北のライン）があります。実験としては十分に感じましたので、この色の調整で次のステップに進みます。

Step5: タイル作成

昨日、LandsatではRGBのVirtual画像をファイルに出力したのですがこれは必須ではなさそうでした。現在表示している4つのレイヤをそのままラスタタイルにします。Processing ToolboxからGenerate XYZ Tiles (Directory)を使います。

出力範囲は4つのシーンをあわせた範囲を考えます。Xmin,Xmax,Ymin,Ymaxと並んでいて、4つのシーンはそれぞれ、
SUE: 300000.0000,409800.0000,3890220.0000,4000020.0000 [EPSG:32654]
SUF: 300000.0000,409800.0000,3990240.0000,4100040.0000 [EPSG:32654]
SVE: 399960.0000,509760.0000,3890220.0000,4000020.0000 [EPSG:32654]
SVF: 399960.0000,509760.0000,3990240.0000,4100040.0000 [EPSG:32654]
という範囲なので、ラスタタイルは、**300000.0000,509760.0000,3890220.0000,4100040.0000 [EPSG:32654]**にします。出力のズームレベルはZL6～ZL13までにしました。

そして出力先フォルダを指定するとデータが出力されます。この範囲のラスタタイルで436MBでした。GitHubでもホストできそうなサイズでした。リーフレットでみるとこんな感じです。30mより10mの方がだいぶよく見えます。

(追記2022-01-22)
この試行では最大ズームレベルを13にしたのですが、元データの解像度が10mなので、ZL14くらいまでにした方が良いかもしれません。ニューヨークの近辺でZL14まで作ってみたのがこちらです。
https://ubukawa.github.io/sentinel2-ny20211019/ny20211019/ny20211019.html
ラスタタイルが256×256として、ZL13だとウェブ地図だと全世界は東西方向が2,097,152ピクセル、ZL14だと4,194,304ピクセルです。赤道の長さがだいたい4万キロとして、赤道上の1ピクセルがZL13だと幅19メートルくらい、ZL14だとはば9.5メートルくらいです。北緯40度くらいで緯線の長さが3万キロと少しと考えると、一ピクセルの幅も赤道上の3/4くらいになります。ZL14が良さそうですね。
（追記ここまで）

Step6: タイルのホスト

GitHubレポジトリにホストしました。こちらです。
https://github.com/ubukawa/sentinel2-kanto2021

タイルのURLは　https://ubukawa.github.io/sentinel2-kanto2021/kanto2021-sentinel2/{z}/{x}/{y}.png　
サンプルを見たい人は　https://ubukawa.github.io/sentinel2-kanto2021/kanto2021-leaflet.html

Step7: ベクトルタイルの重ねあわせ

apLibre gl jsにしろ、Mapbox gl jsにしろ、スタイルファイルにソースとレイヤを追加すればウェブ地図になります。

styleファイルの一部（sourcesの中）

    "r": {
      "type": "raster",
      "tiles": ["https://ubukawa.github.io/sentinel2-kanto2021/kanto2021-sentinel2/{z}/{x}/{y}.png"],
      "attribution": "Copernicus Sentinel data (2021, EUROPEAN COMMISSION) downloaded through U.S. Geological Survey EarthExplorer.",
      "tileSize": 256,
      "maxzoom": 13
    }

styleファイルの一部（layersの中）

    {
      "id": "sentinel2",
      "type": "raster",
      "source": "r"
    }

今回は、地理院地図Vector（仮称）の標準地図風のデザインに追加してみました。みられるでしょうか？
https://ubukawa.github.io/sentinel2-kanto2021/style.json
https://ubukawa.github.io/sentinel2-kanto2021/test-map.html

みられました。ZL13から建物のポリゴンがでてくるので、~~Sentinelは13くらいまででちょうどいいかもしれません~~（2022.1.22修正。解像度から考えると14くらいがいいと思います）。

注意：データの利用について

コペルニクスのデターもオープンライセンスになっているということです（例えばCOPERNICUS SENTINEL SATELLITE IMAGERY UNDER OPEN LICENCEの記事）。ダウンロードしたページでデータの利用条件を確認しました（USGS explanation on Sentinel 2 Terms and Conditions）。

データの出典として、**Copernicus Sentinel data (2021, EUROPEAN COMMISSION) downloaded through U.S. Geological Survey EarthExplorer.**と明記しておきます。

まとめ

昨日紹介した手法で、Landsatにつづき、Sentinel2でもラスタタイルを作ることが確認出来ました。
10メートルの解像度なので、よく見えます。最近のデータを使えるところも魅力的ですね。

（今日の作業、Qiitaも書いたので3時間近くになってしまいました。目標に届かず・・。まだ5時前なので1時間ちょっと二度寝します。）

謝辞

コペルニクスプロジェクトの関係者に感謝します。USGSのEarthExplorerにもいつもお世話になっています。
QGISも使いました。便利で助かります。

参考文献