[English page]

GRASSによるリモートセンシング画像処理:演習3

演習内容:

  1. UTM座標系の衛星画像(Landsat ETM+)の表示、植生指数(NDVI)画像の作成GUIを用いた操作の説明
  2. 等経緯度座標系の標高データ(SRTM)の表示、モザイク合成、3次元表示
  3. ASTER画像からのNDVI、輝度温度の計算、標高データとの3次元表示
  4. 植生・水指数の時系列変化を用いたトンレサップ湖周辺およびメコン川氾濫域の抽出(牧雅康先生)

演習2を踏まえて、今度はUTM座標系のASTER画像を使用して、NDVIと輝度温度を計算し、同じUTM座標系のASTER DEMという標高データに重ね合わせる。

処理の流れ

  1. ASTER画像のインポート
    • バンド2と3 (data1.l3a.vnir2, data1.l3a.vnir3n):NDVI用
    • バンド10〜14 (data1.l3a.tir10 ...):輝度温度用
    • 標高データ (data1.l3a.demzv)

  2. NDVIの計算:DN値から放射輝度へ

  3. 輝度温度の計算:DNから放射輝度を求め、プランクの式を用いて輝度温度へ

  4. 標高データとの重ね合わせ(nvizの起動:NDVIと標高、輝度温度と標高を重ね合わせる)

  5. (余力があれば)バンド10〜14のうち異なるバンドを用いて輝度温度を推定して、標高データとの重ね合わせ。熱バンドの違いによる輝度温度の違いを確認する。

これまでに演習してきた内容の応用なので、細かい手順はここでは示さない。特に必要な注意点を下記に示す。

ASTER画像用のフォルダの作成

演習1は大阪を中心とするETM+画像であり、今回使用するASTER画像は琵琶湖周辺の画像で、当然範囲が異なる。そのため、自分自身で範囲を設定する必要があるだけでなく、ETM+画像用とは異なるフォルダが必要になる。

まずフォルダの設定に関しては、LOCATION, MAPSET, DATABASEの指定において、下記の例のように設定してみよう。

「utm2」というLOCATIONを新規に作成

LOCATION:  utm2
MAPSET:    PERMANENT
DATABASE:  /home/student

次に範囲の設定について述べる。様々な考え方があると思うが、まずは範囲や画像解像度の情報が記載されている"data1.l3a.gh"というファイルを開いて確認しよう。この情報から大体の範囲を予想できるが、厳密な範囲はよく分からない。よって、

適当に範囲を指定(あるいは指定しない) → 後で画像が持つ範囲の情報をGRASS上で読み込む(g.regionを使って)

としても、画像を閲覧することはできる。この方法を採って、今回は下記のように設定してみる。

測地系の指定:

(中略)

Enter Zone: 53(と入力してEnter)

Is this South Hemisphere (y/n) [n] n(と入力してEnter)


                  DEFINE THE DEFAULT REGION


             ====== DEFAULT REGION ======
             | NORTH EDGE: 100     |
             |                          |
  WEST EDGE  |                          |EAST EDGE
  0          |                          |100
             | SOUTH EDGE: 0     |
             ============================

  PROJECTION: 1 (UTM)                   ZONE: 53

                    GRID RESOLUTION
                        East-West:     15
                      North-South:     15

現在の範囲を表示させる方法

> g.region -p

インポートした画像から範囲を取得する方法(例えば、vnir1という名前でインポート済みの状況を想定)

> g.region rast=vnir1

範囲が変更されたか確認

> g.region -p

ASTER画像のDN値から放射輝度への変換

可視から熱バンドの変換係数が、佐賀大学の新井康平先生のPDFファイルに記載されている。こちらの情報を抜き出しておく。(本来は画像のヘッダファイルを参照すべきであるが、ここではこの表の値を用いて演習する)

DN値を放射輝度Rへ変換する式

R = (DN - 1) UCC

  UCC: Unit Conversion Coefficients(変換係数)

UCCの値 (W/m2/str/μm)

注意: バンドごとにHigh gainモードで観測されたか、Normal gainモードかは異なる。その情報はdata1.l3a.ghに記載されているので、ファイルを開いて確認しよう。

Band High gain Normal gain
1 0.676 1.688
2 0.708 1.415
3N/3B 0.423 0.862
4 0.1087 0.2174
5 0.0348 0.0696
6 0.0313 0.0625
7 0.0299 0.0597
8 0.0209 0.0417
9 0.0159 0.0318
10 N/A 0.006822
11 N/A 0.00678
12 N/A 0.00659
13 N/A 0.005693
14 N/A 0.005225

放射輝度から輝度温度への変換

輝度温度:観測物体と等しい放射エネルギーを放射する黒体の温度

黒体:観測物体と等しい放射エネルギーを放射する黒体の温度

(日本リモートセンシング研究会編「図解リモートセンシング」より)

プランクの放射法則:放射輝度と輝度温度の関係

式(1)がプランクの放射法則と呼ばれ、それを変形した式(2)を用いることで、黒体の絶対温度つまり輝度温度が得られる。

Bλ 黒体の分光放射輝度 (W/m2/str/μm)
T 黒体の絶対温度 (K)
λ 波長 (μm)
c 光速度: 2.988 × 108(m/s)
h プランク定数: 6.626 × 10-34(J s)
k ボルツマン定数: 1.380 × 10-23(J/K)

ASTERの熱バンドの波長帯

(注意:下表のK1, K2の値は、単に演習用の数値として、波長帯の中間の値を代入して得られた値にすぎない。実際の各バンドのK1, K2の値とは多少異なる)

バンド 観測波長帯 (μm) 波長帯の中間 (μm) K1: 2hc25 (W/m2/str/μm) K2: hc/k λ (K)
10 8.125 - 8.475 8.300 3.00×103 1.73×103
11 8.475 - 8.825 8.650 2.44×103 1.66×103
12 8.925 - 9.275 9.100 1.90×103 1.58×103
13 10.25 - 10.95 10.60 8.84×102 1.35×103
14 10.95 - 11.65 11.65 6.42×102 1.27×103

カラーテーブルの作成

輝度温度画像にカラーテーブルを作成してみよう。下記に一例が示してあるが、色の割当は自由に変更して試して欲しい。 画像と同じ名前で下記のファイルを作り、下記の場所に保存する。

保存先:

/home/student/utm2/PERMANENT/colr
(例えば「T_tir10」という画像に対してカラーテーブルをつけるのであれば、/home/student/utm2/PERMANENT/colr/T_tir10 という名前をつけて保存)

ファイルの中身:

% 0 300
0:0
270:0:0:0 280:0:0:255
280:0:0:255 290:200:200:0
290:200:200:0 300:255:0:0

[演習1:UTM座標系の衛星画像(Landsat ETM+)の表示、NDVI画像の作成]
[演習2:等経緯度座標系の標高データ(SRTM)の表示、モザイク合成、3次元表示]
[演習3:ASTER画像からのNDVI、輝度温度の計算、標高データとの3次元表示]
[演習4:植生・水指数の時系列変化を用いたトンレサップ湖周辺およびメコン川氾濫域の抽出(牧雅康先生)]
[須崎純一のページへ]
須崎純一 京都大学大学院 工学研究科社会基盤工学専攻 空間情報学講座