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

![rectangle_large_type_2_3ce092ed02fa2d938296553f5ca15a32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/d9c68f0e-de2e-d72c-1018-127fc3f77450.png)

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

https://blog.esrij.com/2014/08/26/arcgis-landsat-890b/

QGISによる地表面温度の算出について参考になるサイトには，
「地表面温度プロダクトを表示する方法（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

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

https://www5.hp-ez.com/hp/tottyiwata/page22/11

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

https://sorabatake.jp/737/

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

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

# 1. Landsat8の衛星画像の取得
## 1.1 Landsatとは？
そもそもLandsatとは何なのか？
Landsatは，アメリカ合衆国の地質調査所（**USGS**：United States Geological Survey）とアメリカ航空宇宙局（**NASA**：National Aeronautics and Space Administration）が管理する**地球観測衛星**である。毎日地球を周回しながら，中分解能マルチスペクトル画像及び，熱赤外画像の取得・保存を行っている。Landsatシリーズは，1972年に1号機が打ち上げられ，2024年現在までで9号機まで打ち上げられている。
![1717408552831-7McoTEUF5v.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/f231015e-fc38-89d6-ae81-d9a8a5c34eef.png)
出典：Landsatの衛星概念図（Credit: NASA）

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

![1717417992272-5Ge8SQu6Ux.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/9b1e8ed2-d21e-b455-1f2b-24b3377d5f13.png)

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

![1717424475200-puPSO5wzoG.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/db582555-ef75-e5d9-9fc0-207da41a1ddc.png)
出典：光の波長と波長による観測しやすい観測対象を表すグラフ（Credit：JAXA）
（https://sorabatake.jp/364/ ）

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

![1717515247043-HdjXg2VdGH.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/fdb7f82a-8f2d-4d59-625c-c92c23fa8dbe.png)

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

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

https://earthexplorer.usgs.gov/ 

![1717418439019-M1Kgufxm2i.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/c5ef9440-afed-3951-16ba-c09e54a6250f.png)

② 検索範囲を設定する。地図上で数点を左クリックしてエリアを選ぶか，座標を直接入力する。今回は熊谷市を含むエリアにする。
![1717418955155-v4W6VzJaFp.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/b456c366-4c85-ee00-2d32-1b26e7b4c843.png)

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

![1717419918026-nBV5JWWxFu.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/bc02f9b8-68bc-58fd-57fe-5de9fdade32e.png)

| タブ | 項目 | 入力値 | 選択 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**」を押して条件に沿った画像を検索する。衛星プロダクトの種類についての説明は割愛する。

![1717421134600-tQIAhZi0kL.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/51eaa6d8-7c3e-221f-9514-3d5e480e48f0.png)

今回の条件だと7枚の画像が合致した。検索結果から，画像の取得日を確認することができる。一番上の画像は**2024年5月18日**であった。7つのアイコンから様々な差操作を行える。
* 左から1つ目のアイコン：画像の範囲の可視化
* 左から2つ目のアイコン：**画像の可視化**
* 左から4つ目のアイコン：**メタデータの表示**
* 左から5つ目のアイコン：**ダウンロード方式の選択**

![1717422376573-0dlI1EfsPU.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/00739e60-b2f7-3b4c-dca3-ae1f6f6a1d7f.png)

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

![1717423031247-FEluHj1q5S.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/fa5a3acd-471d-a3b7-60a9-f20a44b41e92.png)

:::note warn
画像の取得開始時間を示す「**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**」
前者は，先ほど説明した<b>解像度30 m </b>の<b>Band 10（TIRS）</b>である。
後者は，**メタデータが記載されたテキストファイル**であり，後の解析で必要となってくる。

![1717424253057-FHWRQ36eis.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/3ef338dc-75e4-f439-ec37-48c27a7bf478.png)

<br>

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

<img width="450" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/dcb9c54b-dcf5-648e-b0bf-579145a1733a.png">


# 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**

![1717507069514-DeJ6Y1BXZ0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/b3efba4e-8c8b-3930-ad08-bf1d64b7629f.png)

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

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

![1717508847499-JkZFcxx5yn.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/3949eb18-478d-7252-1bd7-3c2ffc8ace3b.png)

## 2.3 カスタム関数の作成

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

![1717505414212-5AoleiKIo6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/7e0a7188-914f-979b-6956-b24b7ea3b364.png)

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

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

![1717509507771-dzX4cPoVyg.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/ff290b74-6485-3e15-9e88-2a55f381e9ef.png)
![1717517323026-y1pAif02c1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/58482e04-aac3-26f9-8539-6ee9f4924514.png)

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

![1717518567576-uikqVFiVnr.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/14b22486-b597-e1aa-8bcc-14eec4b04a7a.png)

④ **追加したバンド演算をダブルクリック**し，プロパティを開く。
| タブ | 項目 | 入力値 | 選択 or 入力 |
|:-----------|:-----------|:------------|:------------|
| 一般 | 名前  | **TOAの計算** | 入力 |
| パラメーター | 方法　| **ユーザー定義** | 選択 |
| パラメーター | バンドインデック　| **0.00038 * B1 + 0.10000** | 入力 |
| 変数 | Rasterのパブリック | **☑チェック**を付ける| 選択 |

![1717519061492-pkaMcqtJRL.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/4d496497-266c-798a-769f-180c70fbc7dc.png)

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

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

<img width="350" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/7c3be3f0-cd82-f949-b7c0-fb64bfdb5fba.png">

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

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

![1717520246429-IIluJBc2LT.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/b1c30092-6f2c-526c-83d0-15575acacf35.png)

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

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

<img width="350" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/c4f0f94d-5a3c-2109-9081-fb0b5af737bf.png">

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

<img width="500" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/8ce8d10f-93e2-0c2c-0f7f-d338dc8551de.png">

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

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

![1717521273748-BgRqeDgq9G.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/0388af2a-e243-0032-f7bc-d0acfd843bea.png)

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

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

<img width="500" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/db1350d6-7ee4-58b2-05da-12d35e391c5b.png">


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

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

![1717521458892-zktCrlDWuj.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/4c29788c-eaa6-2441-8ae5-561c185dbac5.png)

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

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

<img width=600" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/dd22426f-4c10-1a61-a4bd-183f6e31645e.png">

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

<img width=600" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/26a52f85-df68-3d29-eaed-73d00444ebc8.png">

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

<img width=300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/d5a67370-a34a-7a5b-9554-e0b2e1e406e1.png">



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

![1717524158283-3aXbGnDvke.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/5ac55d35-d1b0-8e9a-f507-843a57abf24e.png)

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

![1717524882386-mlSrXAFdr9.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1049445/b762196d-565a-282c-161d-e57cc1a7bcf5.png)

# まとめ
以上で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日確認）