[English page]

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

ここでは、GRASSというGISソフトウェアを使った簡単な演習を説明しています。GRASSではリモートセンシング画像も表示、処理できます。GRASSはGUIでの操作もできますが、細かい処理に関しては単に私が知らないだけかもしれませんが、CUIでのコマンドによる操作を使わざるを得ない時もあります。このWebページではコマンドを用いた操作を説明していきます。

演習内容：

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

GRASSのインストール

Cygwinのsetup.exeファイルのダウンロード：http://cygwin.com/より
setup.exeの起動：インストールに必要なファイルのダウンロード
setup.exeをダブルクリック　→　「Download Without Installing」　→　保存先の指定　→　ミラーサイトの指定
- Cygwin本体：例えばftp://ring.nict.go.jp（どこでもいい）
- GRASS：http://grass.ibiblio.org
（ダウンロード終了後）setup.exeの起動：インストール
setup.exeをダブルクリック　→　「Install from Local Directory」　→　インストールディレクトリの指定　→　ダウンロードしたファイルのディレクトリの指定　→　インストールするパッケージの指定（分からなければ全てを選択）

GRASSの起動

Cygwinの起動：
スタート　→　プログラム　→　Cygwin
X-windowの起動：
$ startx
GRASSの起動：
$ grass62

UTM座標系用のデータフォルダの作成

GRASS起動後にWindow右下の「Projection values」をクリック

LOCATION, MAPSET, DATABASEの指定：

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

と入力したら、Escを押しながらEnter

座標系の指定：UTMを指定

Please specify the coordinate system for location <utm>

A   x,y
B   Latitude-Longitude
C   UTM
D   Other Projection
RETURN to cancel

>C（と入力してEnter）

測地系の指定：

Do you wish to specify a geodetic datum for this location? (y/n) [y] y

Please specify datum name
Enter 'list' for the list of available datums
or 'custom' if you wish to enter custom parameters
Hit RETURN to cancel request

>list（と入力してEnter）

（種々の測地系が表示される）
Short Name  Long Name / Description

（中略）
tokyo   Tokyo       (bessl ellipsoid)
wgs72   WGS_1972    (wgs72 ellipsoid)
wgs84   WGS_1984    (wgs84 ellipsoid)

Please specify datum name
Enter 'list' for the list of available datums
or 'custom' if you wish to enter custom parameters
Hit RETURN to cancel request

>wgs84（と入力してEnter）

Now select Datum Transformation Parameters
Please think carefully about the area covered by your data
and the accuracy you require before making your selection.

Enter 'list' to see the list of available Parameter sets
Enter the corresponding number, or <RETURN> to cancel request

>list（と入力してEnter）

Number Details
---
1    Used in whole wgs84 region
     (PROJ.4 Params towgs84=0.000,0.000,0.000)
     Default 3-Parameter Transformation (May not be optimum for older datums;
 use this only if no more appropriate options are available.)
---

（中略）

>1（と入力してEnter）

Enter Zone: 53（と入力してEnter）

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


                  DEFINE THE DEFAULT REGION


             ====== DEFAULT REGION ======
             | NORTH EDGE: 4000000     |
             |                          |
  WEST EDGE  |                          |EAST EDGE
  380000    |                          |630000
             | SOUTH EDGE: 3700000     |
             ============================

  PROJECTION: 1 (UTM)                   ZONE: 53

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


と入力した後で、Escを押しながらEnter

→　GRASSのWindowが起動する

画像ファイル（Landsat ETM+）の移動

/home/student/data/etmの下にある全てのファイル（p110r036_7t20011015_z53_nn30.tifなど）を、/home/student/utm/PERMANENTの中に移動

p110r036_7t20011015_z53_nn30.tif
p110r036_7t20011015_z53_nn40.tif
p110r036_7t20011015_z53_nn50.tif
p110r036_7x20011015.met.txt

※Windowsのエクスプローラーを使っても、あるいはCygwinのWindow内でコマンド(mv)で指定しても可能

※ETM+ (Enhanced Thematic Mapper Plus) 画像のダウンロード：環境リモートセンシングの演習では既にダウンロード済みでローカルディスクにコピーしてあるが、今後は下記のWebサイトからダウンロードするといい。

Global Land Cover Facility (GLCF) Earth Science Data Interface, University of Maryland

画像ファイル（Landsat ETM+）のインポートと表示

UTM座標系のLandsat ETM+画像を使用して演習を行う。

作業ディレクトリの変更：画像ファイルのあるディレクトリへ移動（GUIを使用する場合には不要）
```
> cd /home/student/utm/PERMANENT
```

画像ファイル（TIFFフォーマット）のインポート（GUIの方が簡単に指定可能）

> r.in.gdal -o input=p110r036_7t20011015_z53_nn30.tif output=b3
> r.in.gdal -o input=p110r036_7t20011015_z53_nn40.tif output=b4
> r.in.gdal -o input=p110r036_7t20011015_z53_nn50.tif output=b5

インポートした画像の表示（GUIの方が簡単に指定可能）

> d.mon x0
> g.region rast=b3
> d.rast b3　（バンド3のみの表示）

> d.rgb red=b3 green=b4 blue=b5　（RGB合成表示）

画像の拡大・縮小表示（GUIの方が簡単に指定可能）
> d.zoom
と入力すれば、後は左クリックで範囲の左上を指定し、範囲の右下の場所で真ん中クリックをすれば拡大される。コマンドを用いたより詳しい操作については、
> d.zoom -help
で確認するように。

画像の演算（NDVIの計算）

植生の粗密と相関のあるNDVI (Normalized Difference Vegetation Index) と呼ばれる植生指数を表示させる。Landsat ETM+画像の場合、バンド3と4から計算できる。

ただし衛星画像に記録されているのはDN (Digital Number) と呼ばれる単なる数値であり、DNから放射輝度、反射率を経由してNDVIを計算するのが一般的である。また放射輝度あるいは反射率は、そのままでは大気上端での放射輝度や反射率を意味するので、大気補正を実施してからNDVIを計算することもある。ここでは大気補正をせず、DNから放射輝度を求め、NDVIを計算する手順を説明する。

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

/home/student/utm/PERMANENTの下にある「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
（中略）
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
（中略）
END_GROUP = MIN_MAX_PIXEL_VALUE

DNから放射輝度への変換（GUIの方が簡単に指定可能）：
DN値が0の時は放射輝度（L3やL4, L5）を0に設定する条件処理を行う
```
> r.mapcalc "L3=if(b3==0, 0, (b3-1)*(152.9+5.0)/(255.0-1.0)-5.0)"
> r.mapcalc "L4=if(b4==0, 0, (b4-1)*(157.4+5.1)/(255.0-1.0)-5.1)"
> r.mapcalc "L5=if(b5==0, 0, (b5-1)*(31.06+1.0)/(255.0-1.0)-1.0)"
```
（GUIでの操作）放射輝度が0になっているかを確認するには、画像を表示させた上で、をクリックしてから画像上でクリックすれば、

569932.356846|3922792.861||0||0||0|

のように、座標とRGBの値が表示される。
放射輝度からNDVIの計算（GUIの方が簡単に指定可能）：
```
> r.mapcalc "ndvi=if(L3==0 && L4==0, -1, (L4-L3)/(L4+L3))"
> d.rast ndvi
```
一面緑色の意味のないような画像が表示されるだろう。しかし、画像上でクリックすれば-1から1までの値が入っていることが分かる。この範囲に合わせたカラーテーブルを作成しよう。
色付けの変更（GUIの方が簡単に指定可能）：
```
> r.colors map=ndvi rules=ndvi
> d.rast map=ndvi
```
上記の1行目のうち、2番目の「ndvi」は予め用意されたカラーテーブルを指定することを意味している。自分で任意のカラーテーブルをファイルとして用意することも可能だが、ここでは説明しない。

最終的に得られるNDVI画像を下記に示す。

参考にしたWebサイト

筑波大学　奈佐原（旧姓西田）顕郎先生のGIS概論

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

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