Landsat 8画像のダウンロードとQGISでラスタタイル化

はじめに

早起きしたので、衛星画像をタイル化することを試してみました。オープンソースのものだけでどれだけできるか試してみました。最新の衛星画像でラスタタイルをつくって背景にするなんて素敵だと思います。

##　環境

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

##　手順
###　Step1: Landsat画像データのダウンロード（USGSのEarthExplorerから）
しばらく使っていなかったので、昔使っていたUSGSのEarthExplorerのアカウントがなくなってしまいました。新たに登録してデータをダウンロードします。

ユーザー登録後、EarthExplorerで、画像がほしい地区の指定、観測日時や雲の割合でのフィルター設定、その後にデータベースを選びます。今回は、1）東京近辺の、2）この1年分くらい、3）そしてLandsat-8のCollection2のlevel2の処理レベルのものを選びます。

Collecation 2 の Level 2の処理レベルは、地形補正等に加えて、輝度値を反射率に変換しているようです。しかも大気上端反射率（TOA）でなくて、地表の反射率（surfece reflectance）にしているということなので驚きです。これをオープンで利用させてくれるなんて素晴らしい。一番上の処理レベルであるU.S. Analysis Ready Dataは日本近辺では使えないのですが、数年前はL1GTの補正レベルのものをダウンロードして、DNから輝度、そして反射率にするプロセスが必要だったのでずいぶん進んだなぁと思います。（もちろん、レベル1GTの補正でも大変感謝していました。）

雲が少なくて、なるべく新しいものということで、パス・ロウが107・035で2021年12月28日の画像をダウンロードすることにしました。（Landsat画像のファイルの命名ルールやWRS2についてはここでは説明しませんが、知っておくと便利です。）
ダウンロードしたファイルはtar圧縮されているファイルで約850MBです。解凍すると以下の様なファイルがあります。

可視光のバンドに加えて、近赤外、中間赤外、品質情報などのバンドもあるのがわかります。データはたくさんありますが、今回の作業でこれから使うのは可視光のバンド2，3，4の3つ画像です。

###　Step2: QGISで画像をみて方針を決める

####　2-1.QGISで画像を開く
QGISで新規のプロジェクトを開いて、そこにファイル名がB2,B3,B4で終わる画像ファイルを追加します。Landsat8のバンドの波長はこちらにありますが、バンド2が青、バンド3が緑、バンド4が赤の画像です。

画像を追加すると、プロジェクトの投影法がEPSG4326から32654に変わったのがわかります（ウィンドウの右下）。これはこの地域のUTMの投影法ですので、ダウンロードしてきた画像がきちんとジオレファレンスされていることがわかります。（Natural Earthの陸地ポリゴンを重ねてみましたが、画像が正しい位置に出てきていました。）

画像を見てみると全体的に色が暗いです。また、値の分布をみると例えば4872～58147となっていて、8bitではないことがわかります。MTLファイルをみると、バンド2,3,4はみなUINT16（符号なしの16ビット整数値,1-65455）と確認出来ます。これについて調べてみると、地表面反射率にスケーリングファクターをかけた値が格納されているようです（詳細：Collecation 2 Level 2の説明）。

####　2-2.方針を決める
この3つの画像からラスタタイルを作るのには以下のことが必要です。

• 各バンドの画像をRGBに当てはめ、カラー合成した画像ファイルを作る。
• 各バンドの描画ストレッチを直して、見やすい色合いにする。そして8bit（×3band）の画像として出力する。
• 出来た画像をラスタタイルにする。

これらはどうやらQGISで全部できそうです。

###　Step 3：画像の合成と色の調整
####　3-1.RGB画像の合成
QGISのラスターメニューに、Build Virtual Rasterというものがあるので選びます（下図）。

実際のコマンドはgdalのようなので、この部分はコマンドラインでも出来そうですね。

バンド4,3,2の順にRGBを重ねたいので、念のためバンドの順番も気をつけます。実行するとVirtural Rasterが出来ます。

####　3-2.色調調整
出来たレイヤのプロパティからヒストグラムを確認します。

ヒストグラムから以下の二つがわかりました。二つ目のことについては、ファイル名で自動的にバンド2～4が割り当てられたのかとも考えましたが、バンドの順番はb4,b3,b2の順に並んでいました。理由はよくわかりません。

• 値がだいたい5,000から12,000に分布していること
• なぜか、band1か空で、band2,band3,band4に値が入っていること

上記の発見を踏まえて、レイヤのプロパティからシンボロジーを調節します。以下の図のとおり、RGBにそれぞれのバンドを割り当て、ストレッチのminとmaxを調節します。true color合成なので、RGBに衛星のRGBをそれぞれ割り当てます。ストレッチは、今回の場合、青から赤にかけてmaxの値を少し増やすことで自然な色合いになりました（青の方が散乱しやすいから暗くなりがちなのかもしれませんね。）。ENVIなどのリモセンソフトでは色のストレッチについていろいろな方法が利用できるのですが、QGISでは高度なストレッチは見つけられませんでした。最大と最小を指定してその間は線形の関数でストレッチをかけるんだと思います。

次に、作成したストレッチを画像として出力します。ラスタメニューのコンバージョンから以下の図の様にやったのですが、もともとの値（16bitのUINT）をもったままの画像になってしまいました。

そこで、レイヤを右クリックして、単純にsave raster layer asを選びます。一番上でraw data ではなくて、rendered imageを選ぶように気をつけましょう。すると画像が出力できます（ 8bit × 3band のカラー画像です）。

###　Step3: QGISでタイルの作成
Step 2から続いてQGISで作業します。QGISのProcessing Toolboxには、なんとラスタタイルを作成するツールがあるそうです。Generate XYZ Tiles (Directory)を選びます。

Extentは画像ファイルの範囲から決めればよいですし、maxとminズームは好きなようにするとよいです。Landsatの解像度が30meterということを考えれば最大ZLは12くらいでいいと思います。出力先のフォルダと、必要に応じてLeafletの地図ファイルを指定して実行を押します。すると、ZXYのベクトルタイルが作成されます。

なお、このGenerate XYZ Tilesは、入力ラスタを指定するわけでないので、どうやら現在描画されているQGISプロジェクトのシーンを出力するのではないかと思います。ですので、Step2でRGB合成画像を必ずしも出力しなくてもいいかもしれませんね。

Step 4: Web地図での表示

Step 4-1：ラスタタイルのアップロード

私の場合は、作成したラスタタイルをGitHubのレポジトリーにアップロードしました。
その後、GitHubページの設定をすることでラスタタイルのホスト完了です。
　　レポジトリはこちら　→　https://github.com/ubukawa/landsat-test01

Step 4-2：ウェブ地図の作成

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

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

    "r": {
      "type": "raster",
      "tiles": ["https://ubukawa.github.io/landsat-test01/test-image01/default/{z}/{x}/{y}.png"],
      "attribution": "Landsat-8 image courtesy of the U.S. Geological Survey",
      "tileSize": 256,
      "maxzoom": 12
    }

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

    {
      "id": "landsat",
      "type": "raster",
      "source": "r",
      "maxzoom": 12
    }

すると、こんな感じでみられます。

実際のファイルはこちら

• styleファイル　https://ubukawa.github.io/landsat-test01/test-image01/test-simple.json
• mapのhtml　https://ubukawa.github.io/landsat-test01/test-image01/test-map.html

注意：コピーライトについて

USGSのページにあるとおり、USGSからダウンロードしたLandsatの利用制限はありません。ただし、引用元を明記することが必要です。（私もウェブ地図に追加したときはattributionでLandsat-8 image courtesy of the U.S. Geological Surveyと書いています。）

まとめ

今回は、パブリックドメインの衛星画像（Landsat-8）とオープンソースのQGISをつかって、ラスタタイルを作成してみました。さらに、それをベクトルタイルに重ねてみました。ソースデータの下処理がかなり高度になっているので、ラスタタイル作成はQGISだけでもできることが確認できました。ただし、より高度な色調補正には専門のリモセンソフトでの対応が必要になるかもしれません。

また、将来に向けて、さらに深掘りするとするとネタは以下の通りです。Sentinel-2での実験は試したのを報告するかもしれませんが、あとはもしかしたら試してみるかもしれないくらいです。

今後のネタ

Sentinel-2での実験

じつはもう試しています。解像度が10mになるのでめちゃめちゃ面白かったです。Sentinel-2もUSGSとESAで協定があるようで、EarthExplorerからダウンロードできます（全てではないようです）。

画像のモザイク

時期やシーンが違う画像の間でモザイクするのにどのくらいの手間が必要なのか、QGISでの処理だけで大丈夫かは試してもよいかもしれません。Landsatの処理レベルがCollection 2 level 2になって、surface reflectanceになっているので、昔よりは手間が小さいかもしれません。

処理の自動化

対象が小さいと今回のように手動でいいのですが、多数の処理をしたい場合にはコマンドラインでの処理なども検討が必要かもしれません。一部の処理はGDALで出来るかもしれませんが、色調補正はどうしても画像を確認しながらという部分があるので自動化しにくいところがあるかもしれません。
（Surface Reflectanceとscaling factorが上手くいっていればもしかしたら自動化できるかも）

雲の除去

モザイクとも関係しますが、雲の除去をしたい人もいるかもしれませんね。私の実験では使う画像がいつであってよかったので、単純にシーン全体が綺麗に見えている日のものを選びました。雲除去のことは特に気にしていません。LandsatのQAのデータを使うとか、複数シーンをコンポジットするとか、いくつか方法はあるのではないかと思います。（データサイズが大きいから大変ですね。）

謝辞

今回の実験では、U.S.Geological Surveyからダウンロードした　Landasat-8画像を利用しました。こんなにいいデータを自由に利用させてくれるUSGSに感謝します。
また、QGISの機能もすばらしいと思います。開発された皆様に感謝します。

参考

• USGS EarthExplorer
• Landsatの処理レベル概要
• Level2 Science Product
• Level2 Science ProductのScaling factorなど
• Landsat band