見出し画像

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

you

はじめに

本記事では，広域における熱環境を分析する前段階として，
ArcGIS Proを用いてLandsat8から地表面温度を推定する方法について
備忘録的に解説していく。
ArcGISによる地表面温度の算出について参考になるサイトには，
「ArcGIS で Landsat 画像を地表面温度画像に変換しよう！（ArcGISブログ）」などがあり，
QGISによる地表面温度の算出について参考になるサイトには，
「地表面温度プロダクトを表示する方法（G-Portal）」
「第2節　QGISで地表温度を計算しよう（岩田　敏彰ホームページ）」
「東京五輪マラソンは本当に暑い？ 衛星データで過去と比較してみた（宙畑）」などがある。

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


1. Landsat 8画像の取得

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 Landsat 8のバンドと解像度

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

Landsat8の衛星諸元
（出典）光の波長と波長による観測しやすい観測対象を表すグラフ（Credit：JAXA）
（光の波長って何？ なぜ人工衛星は人間の目に見えないものが見えるのか / 宙畑）

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

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

Landsat8のBand10の概念図

1.3 EarthExplorerからの衛星画像取得

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

(1) EarthExplorer（https://earthexplorer.usgs.gov/）にアクセスする。
　　サイトから衛星画像をダウンロードをする際には，アカウントの登録が
　　必要である。新たなユーザーアカウントを作成し，ログインする。

アカウントの作成とログイン

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

衛星画像の欲しいエリアの選択

(3) 画面左下にあるData Rangetタブの，
　　Search from を「01/01/2024」（mm / dd / yyyy 形式で記入すること）
　　Search to を「06/03/2024」（mm / dd / yyyy 形式で記入すること）
　　画面左下にあるCloud Cover タブの，
　　Cloud Cover Range を「0 % - 10 %」
　　に設定。ここでは，2024年1月1日～2024年6月4日までに撮影された
　　衛星画像の中で，雲の被覆率が0～10 % の画像のみを選択するという
　　作業を行っている。

日付（Data Rangetタブ）と雲の被覆率（Cloud Cover タブ）の設定

(4) 画面左上にある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」を開き，
①「LC09_L1TP_107035_20240518_20240518_02_T1_B10.TIF」と
②「LC09_L1TP_107035_20240518_20240518_02_T1_MTL.txt」を
ダウンロードする。
①は，先ほど説明した解像度30 m の Band 10（TIRS）である。
②は，メタデータが記載されたテキストファイルであり，後の解析で必要となってくる。

Band10とテキストファイルのダウンロード

※ 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 のアイコンの上に持っていき，アイコンの上で離す。展開先のフォルダを指定する。

7-zip を用いたアーカイブファイルの展開

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 毎に確認する必要がある。
今回は以下の数値であることを確認した。

Band 10 に関するパラメータ

2.2 地表面温度の算出方法

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

式(1)～(3)と変数

2.3 カスタム関数の作成

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

Band 10 の表示と名前の変更

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

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

ラスター関数エディター
ラスター関数テンプレートウィンドウとアイコンの説明

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

ラスター関数テンプレートウィンドウへのバンド演算の追加

④ 追加したバンド演算をダブルクリックし，プロパティを開く。
　 一般タブの名前を「TOAの計算」
　 パラメータータブの方法を「ユーザー定義」
　 パラメータータブのバンドインデックを「0.00038 * B1 + 0.10000」
　 （B1＝一つ目のBandのDN値，マルチバンドなら熱赤外バンド数に変更）
　 変数タブのRasterのパブリックに「☑チェック」
　 に設定しOKを押す。
　 そして緑色のラスタアイコンからTOAの計算アイコンまで矢印を伸ばす。

TOAの計算に関するパラメータの設定
緑色のラスタアイコンからTOAの計算アイコンまで矢印を伸ばす

⑤ ③と同様に，右側のラスター関数ウィンドウからバンド演算を
　 ドラッグ&ドロップで左側のラスター関数テンプレートウィンドウに
　 追加する。
　 ④と同様に，追加したバンド演算をダブルクリックしプロパティを開く。
　 一般タブの名前を「Kの計算①」
　 パラメータータブの方法を「ユーザー定義」
　 パラメータータブのバンドインデックを「(799.0284 / B1)+ 1」
　 （※ここのB1はマルチバンドであっても変更しない）
　 に設定しOKを押す。
　 そして，TOAの計算アイコンからKの計算①アイコンまで矢印を伸ばす。

Kの計算①に関するパラメータの設定
TOAの計算アイコンからKの計算①アイコンまで矢印を伸ばす

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

Kの計算①アイコンからLnアイコンまで矢印を伸ばす

⑦ ③と同様に，右側のラスター関数ウィンドウからバンド演算を
　 ドラッグ&ドロップで左側のラスター関数テンプレートウィンドウに
　 追加する。
　 ④と同様に，追加したバンド演算をダブルクリックしプロパティを開く。
　 一般タブの名前を「Kの計算②」
　 パラメータータブの方法を「ユーザー定義」
　 パラメータータブのバンドインデックを「1329.2405 / B1」
　 （※ここのB1はマルチバンドであっても変更しない）
　 に設定しOKを押す。
　 そして，LnアイコンからKの計算②アイコンまで矢印を伸ばす。

Kの計算②に関するパラメータの設定
LnアイコンからKの計算②アイコンまで矢印を伸ばす

⑧ ③と同様に，右側のラスター関数ウィンドウからバンド演算を
　 ドラッグ&ドロップで左側のラスター関数テンプレートウィンドウに
　 追加する。
　 ④と同様に，追加したバンド演算をダブルクリックしプロパティを開く。
　 一般タブの名前を「Tの計算」
　 パラメータータブの方法を「ユーザー定義」
　 パラメータータブのバンドインデックを「B1 - 273.15」
　 （※ここのB1はマルチバンドであっても変更しない）
　 に設定しOKを押す。
　 そして，Kの計算②アイコンからTの計算アイコンまで矢印を伸ばす。

Tの計算に関するパラメータの設定
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/）」を重ね，
レイアウトに出力すると以下のような
関東近郊の地表面温度を色分けした画像が完成する。

2024年5月18日の関東近郊の地表面温度

まとめ

以上で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日確認）



いいなと思ったら応援しよう！