【ArcGIS Pro】Landsat8を用いた地表面温度解析

Last updated at 2024-06-05Posted at 2024-06-04

はじめに

本記事では，広域における熱環境を分析する前段階として，ArcGIS Proを用いてLandsat8から地表面温度を推定する方法について備忘録的に解説していく。
ArcGISによる地表面温度の算出について参考になるサイトには，
「ArcGIS で Landsat 画像を地表面温度画像に変換しよう！（ArcGISブログ）」

QGISによる地表面温度の算出について参考になるサイトには，
「地表面温度プロダクトを表示する方法（G-Portal）」

「第2節　QGISで地表温度を計算しよう（岩田　敏彰ホームページ）」

「東京五輪マラソンは本当に暑い？衛星データで過去と比較してみた（宙畑）」

などがある。
本記事はそれらの二番煎じにはなるが，ArcGIS Pro初心者でも進められるように丁寧に解説していく。地表面温度の算出ということで，本記事では「あついぞ！熊谷！」でおなじみ
埼玉県熊谷市を対象に，Landsat8の衛星画像から地表面温度ラスタを作成していく。

本記事を参考にした際には，引用文献として明示してください。
本記事内で作成した解説画像の改変・再配布は認めません。

1. Landsat8の衛星画像の取得

1.1 Landsatとは？

そもそもLandsatとは何なのか？
Landsatは，アメリカ合衆国の地質調査所（USGS：United States Geological Survey）とアメリカ航空宇宙局（NASA：National Aeronautics and Space Administration）が管理する地球観測衛星である。毎日地球を周回しながら，中分解能マルチスペクトル画像及び，熱赤外画像の取得・保存を行っている。Landsatシリーズは，1972年に1号機が打ち上げられ，2024年現在までで9号機まで打ち上げられている。

出典：Landsatの衛星概念図（Credit: NASA）

1.2 Landsat8のバンドと解像度

Landsat8は，2013年2月11日に打ち上げられ，2013年5月30日からの画像を集めている。地上700 kmを周回しており，回帰日数（地点Aを観測してから次に地点Aを観測するまでの期間）は16日である。Landsat8の簡単な衛星諸元を以下の表に示す。

解像度（分解能）とは，1つのピクセルが表す実際の地上の広がりである。解像度が小さいほど，細かい地物を識別することができる。Landsat8の分解能は Band Number によって異なるが，おおよそ15～30 mに収まる。そのため，道路・田畑など比較的大きな地物は識別できるが，建物・人・車など小さいものは識別できない。例えば，可視線のBand2（Blue）では，30 m × 30 m を1ピクセルとして表現されている。
各ピクセルには，デジタル値（DN：Digital Number）が格納されている。

出典：光の波長と波長による観測しやすい観測対象を表すグラフ（Credit：JAXA）
（https://sorabatake.jp/364/ ）

本記事のような表面温度の推定を行う際には，解像度30 m のBand 10（TIRS）を用いることが多い。しかしながら，範囲を狭めた局所的な熱的環境は把握できない。研究対象や用途に応じて，適切な衛星を選択する必要がある。

1.3 EarthExplorerからの衛星画像取得

Landsat8のデータは，EarthExplorer・GloVis・LandsatLook Viewerから無料でダウンロードすることができる。本記事では，EarthExplorer から衛星画像を取得する。以下の①～④は，EarthExplorer(https://earthexplorer.usgs.gov/ )から衛星画像を取得する手順である。

① EarthExplorer から衛星画像をダウンロードをする際には，アカウントの登録が必要である。新たなユーザーアカウントを作成し，ログインする。

② 検索範囲を設定する。地図上で数点を左クリックしてエリアを選ぶか，座標を直接入力する。今回は熊谷市を含むエリアにする。

③ 画面左下にあるData RangetタブとCloud Coverタブから，日付範囲と雲の被覆率を設定する。本記事では，2024年1月1日～2024年6月4日までに撮影された衛星画像の中で，雲の被覆率が0%～10 %の画像のみを選択してく。

タブ	項目	入力値	選択 or 入力
Data Ranget	Search from	01/01/2024（mm / dd / yyyy 形式）	入力
Data Ranget	Search to	06/03/2024（mm / dd / yyyy 形式）	入力
Cloud Cover	Search to	0 % - 10 %	選択

④ 画面左上にあるData Setsタブから，
「Landsat ➡ Landsat Collection 2 Level-1 ➡ Landsat 8-9 OLI/TIRS C2 L1」を探し，
チェックマークを入れる。そして，「Results」を押して条件に沿った画像を検索する。衛星プロダクトの種類についての説明は割愛する。

今回の条件だと7枚の画像が合致した。検索結果から，画像の取得日を確認することができる。一番上の画像は2024年5月18日であった。7つのアイコンから様々な差操作を行える。

左から1つ目のアイコン：画像の範囲の可視化
左から2つ目のアイコン：画像の可視化
左から4つ目のアイコン：メタデータの表示
左から5つ目のアイコン：ダウンロード方式の選択

まず，左から4つ目のアイコンを選択しメタデータの表示を行う。メタデータを見ると，画像に関する様々な情報を確認することができる。

画像の取得開始時間を示す「Start TIme」は世界標準時（GMT：Greenwich Mean Time）で表記されているため，注意が必要である。
日本時間で考える際には +9時間 して考えてほしい。

メタデータの確認が終わったら，左から5つ目のアイコンからダウンロード方式を選択する。Product Optionsを開き次の2つのファイルをダウンロードする。
「LC09_L1TP_107035_20240518_20240518_02_T1_B10.TIF」「LC09_L1TP_107035_20240518_20240518_02_T1_MTL.txt」
前者は，先ほど説明した解像度30 m のBand 10（TIRS）である。
後者は，メタデータが記載されたテキストファイルであり，後の解析で必要となってくる。

7-Zipを用いたLandsat8画像の展開について

もし，全てのデータ（1.14GB）をダウンロードした場合は，アーカイブファイル（展開前のファイル）の展開が必要となる。ダウンロードしたファイル名を確認すると「LC09_L1TP_107035_20240518_20240518_02_T1.tar」となっている。「.tar」とは，複数のファイルがまとまったアーカイブファイルの拡張子を表している。本記事では，Windowsに標準搭載の解凍ソフトではなく，圧縮・解凍ソフトの「7-Zip」を用いる。

公式HP（https://7-zip.opensource.jp/ ）から
「7-Zip 24.05 (2024-05-14) Windows x64用 (64-bit)」をダウンロードし，ファイルを開いてインストールを行う。そして，展開したいファイルを7-Zip File Manager のアイコンの上に持っていき，アイコンの上で離す。展開先のフォルダを指定する。

2. ArcGIS Proを用いた地表面温度の算出

2.1 メタデータ情報の確認

ダウンロードしたファイルのうち，テキストファイル（本記事通りに行うと「LC09_L1TP_107035_20240518_20240518_02_T1_MTL.txt」）を開く。このファイルには，メタデータの情報が格納されている。このうちBand 10に関する以下のパラメータを探し，数値をメモしておく。

RADIANCE_MULT_BAND_10
RADIANCE_ADD_BAND
K1_CONSTANT_BAND_10
K2_CONSTANT_BAND_10

この値は，固定ではなく常に変わる可能性があるため，
複数枚の場合は，画像 & Band 毎に確認する必要がある。

2.2 地表面温度の算出方法

Landsat8 の Band 10 から地表面温度を算出するための式は，次の式(1)～(3)である。

2.3 カスタム関数の作成

ArcGIS Proを開き1.3節で取得したBand 10の「LC09_L1TP_107035_20240518_20240518_02_T1_B10.TIF」をドラッグ&ドロップで表示させる。名前が長いため，ラスタ名を右クリックしてプロパティを開き，「Band10」という名前に変更する。

2.2節の式をもとに，ラスター関数で計算してくのだが，画像の枚数が増えるほどラスター関数を使う回数（処理の回数）が増えて大変である。そこで，式(1)～(3)を一括して計算できるようなオリジナルの関数を作成し，処理手順を簡略化していく。
手順は①～⑨まである。

① 画像タブを開き，関数エディターを選択するラスター関数テンプレートウィンドウの上部のにあるアイコンのうち，「自動レイアウト」「ラスター変数の追加」「名前を付けて保存」を使っていく。

② 「ラスター変数の追加」を左クリックすると緑色のラスタアイコンが追加される。追加したラスタアイコンを右クリックし，名前の変更から「Band 10」に変更する。
③ 右側のラスター関数ウィンドウから，「システム」➡「数学関数」➡「バンド演算」を探す。バンド演算のアイコンを持ち，ドラッグ&ドロップで左側のラスター関数テンプレートウィンドウに追加する。

④ 追加したバンド演算をダブルクリックし，プロパティを開く。

タブ	項目	入力値	選択 or 入力
一般	名前	TOAの計算	入力
パラメーター	方法	ユーザー定義	選択
パラメーター	バンドインデック	*0.00038 B1 + 0.10000**	入力
変数	Rasterのパブリック	☑チェックを付ける	選択

バンドインデックの入力値にあるB1は，1つ目のBandのDN値を表す。
マルチバンドなら熱赤外バンドの数値に変更する必要がある。

上記に設定しOKを押す。
そして緑色のラスタアイコンからTOAの計算アイコンまで矢印を伸ばす。

⑤ ③と同様に，右側のラスター関数ウィンドウからバンド演算をドラッグ&ドロップで左側のラスター関数テンプレートウィンドウに追加する。④と同様に，追加したバンド演算をダブルクリックしプロパティを開く。

タブ	項目	入力値	選択 or 入力
一般	名前	Kの計算①	入力
パラメーター	方法	ユーザー定義	選択
パラメーター	バンドインデック	(799.0284 / B1)+ 1	入力

ここのB1は，仮にマルチバンドであっても変更しない。

上記に設定しOKを押す。
そしてTOAの計算アイコンからKの計算①アイコンまで矢印を伸ばす。

⑥ ③と同様に，右側のラスター関数ウィンドウからLnアイコンをドラッグ&ドロップで左側のラスター関数テンプレートウィンドウに追加する。そして，Kの計算①アイコンからLnアイコンまで矢印を伸ばす。

⑦ ③と同様に，右側のラスター関数ウィンドウからバンド演算をドラッグ&ドロップで左側のラスター関数テンプレートウィンドウに追加する。④と同様に，追加したバンド演算をダブルクリックしプロパティを開く。

タブ	項目	入力値	選択 or 入力
一般	名前	Kの計算②	入力
パラメーター	方法	ユーザー定義	選択
パラメーター	バンドインデック	1329.2405 / B1	入力

ここのB1は，仮にマルチバンドであっても変更しない。

上記に設定しOKを押す。
そしてLnアイコンからKの計算②アイコンまで矢印を伸ばす。

タブ	項目	入力値	選択 or 入力
一般	名前	Tの計算	入力
パラメーター	方法	ユーザー定義	選択
パラメーター	バンドインデック	B1 - 273.15	入力

ここのB1は，仮にマルチバンドであっても変更しない。

上記に設定しOKを押す。
そしてKの計算②アイコンからTの計算アイコンまで矢印を伸ばす。

⑨ 全てのアイコンを選択し，「自動レイアウト」できれいに整列させる。「名前を付けて保存」を押し，名前を付けて保存ウィンドウの名前を「地表面温度」とする。
これでカスタム関数が完成する。

2.4 カスタム関数を使った地表面温度の算出

2.3節でカスタム関数の保存が完了すると，右側のラスター関数ウィンドウから，
「カスタム」➡「Custom1」➡ 「地表面温度」関数が表示される。
パラメータータブのBand 10 を2.3節で追加し名前を変更した「Band10」に設定し実行。

3. 結果と可視化

作成した「地表面温度_Band10」のシンボルを変更し，埼玉県熊谷駅付近まで拡大してみる。熊谷駅からラグビーロードを中心に高温域が広がっているのがわかる。本記事で使用したLandsat 8 のBand 10 の解像度は30 mであるため局所的な解析にはあまり向いてはいない。

esriジャパンの「全国市町村界データ（https://www.esrij.com/products/japan-shp/ ）」を重ね，レイアウトに出力すると関東近郊の地表面温度を色分けした画像が完成する。

まとめ

以上でLandsat8を用いた地表面温度の算出に関する解説を終わりにする。本記事では紹介しなかったが，Landsat 8 画像の取得と地表面温度ラスタへの変換をプログラムを用いて行う方法もある。次の機会に解説していきたい。

参考文献

[1] ArcGISブログ，ArcGIS で Landsat 画像を地表面温度画像に変換しよう！．
https://blog.esrij.com/2014/08/26/arcgis-landsat-890b/ ．（2024年6月5日確認）

[2] G-Portal，地表面温度プロダクトを表示する方法．
https://gportal.jaxa.jp/gpr/assets/mng_upload/COMMON/upload/Imaging_procedure_for_Level-2_LST_product_of_SHIKISAI(GCOM-C)_using_QGIS_jp.pdf ．（2024年6月5日確認）

[3] 宇宙技術汎用化トレーナー　岩田　敏彰のホームページ，第4章　地表温度を観察してみよう　第2節　QGISで地表温度を計算しよう．
https://www5.hp-ez.com/hp/tottyiwata/page22/11 ．（2024年6月5日確認）

[4] 宙畑，東京五輪マラソンは本当に暑い？衛星データで過去と比較してみた．
https://sorabatake.jp/737/ ．（2024年6月5日確認）

[5] 宙畑，光の波長って何？なぜ人工衛星は人間の目に見えないものが見えるのか．
https://sorabatake.jp/364/ ．（2024年6月5日確認）

[6] USGS，EarthExplorer．
https://earthexplorer.usgs.gov/ ．（2024年6月5日確認）

[7] 7-zip，7-zip．
https://7-zip.opensource.jp/ ．（2024年6月5日確認）

[8] esriジャパン，全国市町村界データ．
https://www.esrij.com/products/japan-shp/ ．（2024年6月5日確認）

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up