植生指数 (NDVI) の計算、表示

ここでは、GRASSというGISソフトウェアを使った簡単な演習を説明しています。GRASSではリモートセンシング画像も表示、処理できます。このページは、GRASSのGUIでの操作を説明していきます。

[GRASSのインストール、標高データを用いた地滑り危険度マップの作成]
[植生指数 (NDVI) の計算、表示]
[標高データ (SRTM)の表示、植生指数 (NDVI) の3次元表示]
[反射率、輝度温度、標高データを用いた土地被覆分類]
[QGIS, Rを用いた公示地価データの空間統計分析：空間的自己回帰モデル]
[QGIS, Rを用いた公示地価データの空間統計分析：静的な時空間モデリング]

本演習内容：

この演習で最終的に得られる画像

GRASSのインストール

QGISとセットになったGRASSのインストーラーをダウンロードする。

インストーラーのダウンロード：download.osgeo.org/osgeo4w/osgeo4w-setup.exeより
osgeo4w-setup.exeの起動：インストールに必要なファイルのダウンロード
osgeo4w-setup.exeをダブルクリック
→　エクスプレス　デスクトップ　インストール
→　ミラーサイトの指定（例えば、「http://download.osgeo.org」）
→　パッケージの選択：「QGIS」「GDAL」「GRASS GIS」

GRASSの起動と設定

プログラムの中の「OSGeo4W」の中に「GRASS GIS」のアイコンが生成されている。それを選択して、GRASSを起動する。

GRASS起動後

GRASS GIS database directoryの指定：
データを保存するディレトクリ（フォルダ）を指定。任意の場所で構わないが、スペースなしの英語名のディレトクリを指定する必要がある。一例を下記に示す。
```
C:\GRASS_data
```
データベースとロケーションの指定：
```
GISデータディレクトリ:   C:\GRASS_data
プロジェクトロケーション:ETM
Location Title:          PERMANENT
```
「デフォルト領域と解像度をセット」のボックスにはチェックを入れない。
新しいロケーションの作成方法を選択
EPSG (European Petroleum Survey Group) コードを指定
→　EPSGコード：32653（WGS84楕円体、UTM投影法、ゾーン番号53）を入力

「Start GRASS session」というボタンをクリックすれば、GRASSの操作画面が表示される。

使用する衛星画像

本講義用の衛星画像のファイル（p110r036_7t20011015_z53_nn30.tifなど）を、自分のGRASS用のディレクトリ（例えば、C:\GRASS_data\ETM\PERMANENTの中にコピーする。

2001/10/15に観測された画像（衛星名：Landsat7, センサ名：ETM+）
- p110r036_7t20011015_z53_nn30.tif (バンド3の画像)
- p110r036_7t20011015_z53_nn40.tif (バンド4の画像)
- p110r036_7t20011015_z53_nn50.tif (バンド5の画像)
- p110r036_7x20011015.met.txt (メタデータ)
1989/05/31に観測された画像（衛星名：Landsat5, センサ名：TM）
- p110r36_5t19890531_nn3.tif (バンド3の画像)
- p110r36_5t19890531_nn4.tif (バンド4の画像)
- p110r36_5t19890531_nn5.tif (バンド5の画像)
- p110r36_5t19890531.met (メタデータ)

※ETM+ (Enhanced Thematic Mapper Plus) 画像のダウンロード：本演習では既にダウンロード済みで上記URLにコピーしてあるが、今後は下記のWebサイトからダウンロードするといい。

Global Land Cover Facility (GLCF) Earth Science Data Interface, University of Maryland
（上記がアクセスできない時は、FTPサイトを利用する。但し、対象地域のパス、ロウ番号を事前に把握しておく必要がある） ftp://ftp.glcf.umd.edu/glcf/Landsat/WRS2/, University of Maryland

衛星画像のインポート

Landsat ETM+画像の、バンド3、4、5を、それぞれ「b3」「b4」「b5」という名前でインポートする。下記にコマンドを示すが、GUI上で処理する。ちなみに、ETM+画像のバンド3、4、5は、それぞれ緑、赤、近赤外の波長に対応している。

「ファイル」→ 「ラスターデータのインポート」 → 「Import of common raster formats [r.in.gdal]」

現れた画面上で下記のように、バンド3の画像を「b3」としてインポートする。

[入力項目]

Name of raster file to be imporated: p110r036_7t20011015_z53_nn30.tif
出力ラスターマップ名: b3
（既存のファイルを上書きする場合）「オプション」→「既存のファイルに上書きする」をチェック

Name of raster file to be imporated: p110r036_7t20011015_z53_nn40.tif
出力ラスターマップ名: b4
Name of raster file to be imporated: p110r036_7t20011015_z53_nn50.tif
出力ラスターマップ名: b5

範囲の設定

画像ファイルの範囲をデフォルト領域に設定する。

「設定」→「Computataional region」→「領域設定 [g.region]」

[入力項目]

「既存領域」→ Set region to match raster map(s): （例えば）b3@PERMANENT
「プリント」→ 「現在の領域を表示」をチェック

注意：この範囲設定を正しく行わないと、領域の範囲が行方向（column）に1画素、列方向（row）に1画素しか存在しない設定になっているため、以下の処理が実行できません。

放射輝度の計算

衛星画像に記録されているのはDigital Number (DN) と呼ばれる単なる数値である。DN値自体には意味がなく、DN値から変換される放射輝度、反射率は地表面の反射特性を表す物理量であり、地表面の種類や状態を推定するのに活用される。また、植生の密生度合いを反映する植生指数の一つであるnormalized difference vegetation index (NDVI) は放射輝度や反射率から計算するのが一般的である。また放射輝度あるいは反射率は、そのままでは大気上端での放射輝度や反射率を意味するので、大気補正を実施してからNDVIを計算することもある。ここではまずDN値から放射輝度を求める手順を説明する。

メタファイル（ヘッダファイル）から変換係数の抽出：

「p110r036_7x20011015.met.txt」というファイルをテキストエディタで（Wordpadでも何でも）開いてみる。すると、下記の記述を見つけるだろう。

GROUP = MIN_MAX_RADIANCE 
（中略）
	LMAX_BAND3 = 152.900
	LMIN_BAND3 = -5.000
	LMAX_BAND4 = 157.400
	LMIN_BAND4 = -5.100
	LMAX_BAND5 = 31.060
	LMIN_BAND5 = -1.000
	LMAX_BAND61 = 17.040
	LMIN_BAND61 = 0.000
	LMAX_BAND62 = 12.650
	LMIN_BAND62 = 3.200
（中略）
END_GROUP = MIN_MAX_RADIANCE 
GROUP = MIN_MAX_PIXEL_VALUE
（中略）
	QCALMAX_BAND3 = 255.0
	QCALMIN_BAND3 = 1.0 
	QCALMAX_BAND4 = 255.0
	QCALMIN_BAND4 = 1.0
	QCALMAX_BAND5 = 255.0
	QCALMIN_BAND5 = 1.0
	QCALMAX_BAND61 = 255.0
	QCALMIN_BAND61 = 1.0
	QCALMAX_BAND62 = 255.0
	QCALMIN_BAND62 = 1.0
（中略）
END_GROUP = MIN_MAX_PIXEL_VALUE

DN値から放射輝度への変換：
Landsat ETM+の場合、以下の式に従ってDN値から（大気上端の）放射輝度へ変換する（一般的に下記の式のようにgainとoffsetを用いた線形式で記述される）。

上記の式をGRASS上で記述し、変換する。DN値が0の時は放射輝度（L3やL4, L5）を0に設定する条件処理を行う。

「ラスター」→「ラスターマップカリキュレータ」→「ラスターマップカリキュレータ [r.mapcalc]」
[入力項目]
- 出力：L3
- 式：if(b3==0, 0, (b3-1.0)/(255.0-1.0)*(152.9+5.0)-5.0)
- 出力：L4
- 式： if(b4==0, 0, (b4-1.0)/(255.0-1.0)*(157.4+5.1)-5.1)
- 出力：L5
- 式： if(b5==0, 0, (b5-1.0)/(255.0-1.0)*(31.06+1.0)-1.0)

反射率の計算

Landsat ETM+に限らず、全ての衛星画像において、以下の式に従って大気上端の放射輝度から大気上端の反射率へ変換できる。大気上端 (top of atmosphere) ということで、TOA radiance, TOA reflectanceと表現されることもある。

ここでは、バンド3, 4, 5の反射率r3, r4, r5を求めてみる。

放射輝度から反射率への変換式

ρ: 反射率
d: earth-sun distance in astronomical units（地球―太陽間天文距離）：2001/10/15の時は、d=0.997197
DOY: 1月1日を1とする通算の日：2001/10/15なので、DOY=288
E_λ: mean solar exoatmospheric irradiances（大気圏外太陽照度）
θ: solar zenith angle in degrees（太陽天頂角）：計算が少し大変なので、ここではp110r036_7x20011015.met.txtに記載されているデータを使って計算した、θ = 46.68 (deg)という値を使用する。

表1: Landsat7 ETM+の大気圏外太陽照度E_λ
バンド	E_λ (W/m²/μm)
1	1970
2	1842
3	1547
4	1044
5	225.7
7	82.06

上記のパラメータを利用して、放射輝度と同様に「r.mapcalc」という関数を用いて処理する。

「ラスター」→「ラスターマップカリキュレータ」→「ラスターマップカリキュレータ [r.mapcalc]」

[入力項目]

出力：r3
式：3.14159*L3*0.997197^2/1547/cos(46.68*3.14159/180)

出力：r4
式：3.14159*L4*0.997197^2/1044/cos(46.68*3.14159/180)
出力：r5
式：3.14159*L5*0.997197^2/225.7/cos(46.68*3.14159/180)

NDVIの計算

Normalized Difference Vegetation Index (NDVI) という植生の密生度合いを反映した指標を計算する。NDVIは-1から1の値を取り、植生が多いほど1に近い値を取る。

「ラスター」→「ラスターマップカリキュレータ」→「ラスターマップカリキュレータ [r.mapcalc]」

[入力項目]

出力：ndvi
式：if(r3==0 && r4==0, -1, (r4-r3)/(r4+r3))

この処理後に、GRASSに既に用意された、NDVIに対応したカラーテーブルを適用する。

色付けの変更：

「ラスター」→「カラー調整」→「カラーテーブル [r.colors]」

[入力項目]

「Map」→ Name of raster map(s)：ndvi
「Define」→ 「Name of color table」：「ndvi」を選択

最終的に得られるNDVI画像を下記に示す（クリックすると拡大画像が表示される）。

課題

これまで2001年の植生指数NDVIを計算してきた。今度は、1989年の画像を使って、1989年のNDVIを計算せよ。
但し、1989年の画像（Landsat5/TM）のメタデータには、DN値から放射輝度に変換する変換係数が記されていない。そこで、関連する論文（Chander and Markham, IEEE Trans Geoscience and Remote Sensing, 2003）から得た変換係数を下記に記す。これを用いて、1989年のNDVIを計算せよ。

表2: Landsat5 TMのDN値から放射輝度への変換係数
バンド	L_min,λ (W/m²/μm)	L_max,λ (W/m²/μm)	Q_calmax	Q_calmin
1	-1.52	193.0	255	0
2	-2.84	365.0	255	0
3	-1.17	264.0	255	0
4	-1.51	221.0	255	0
5	-0.37	30.2	255	0
7	-0.15	16.5	255	0

表3: Landsat5 TMの大気圏外太陽照度E_λ
バンド	E_λ (W/m²/μm)
1	1957
2	1826
3	1554
4	1036
5	215.0
7	80.67

1989年のNDVIに比べて、2001年のNDVIが大きく減少した地域を知りたい。「ラスターマップカリキュレータ [r.mapcalc]」と「カラーテーブル [r.colors]」の機能を利用して、0.2以上減少した地域を赤色、0.1以上0.2未満減少した地域を青色として出力した画像を生成せよ。
上記のNDVIが大きく減少した要因について考察せよ。

[講義「リモートセンシングと地理情報システム(GIS)」のトップページ]

[須崎純一のページへ]

須﨑純一　京都大学大学院工学研究科社会基盤工学専攻空間情報学講座