GRASSを用いた地理情報システム入門
第９回目は，講義で説明したマルチスペクトルを理解するために，前回ダウンロードしていただいたＴＭ画像のバンド３とバンド４を用いて実習します．また，前回は利用しませんでしたが，その中にいれてあったＤＥＭ（地形標高データ；ただし，全域ではありません）も用います．なお，以下の実習の結果として作成される画像は，基本的に講義で紹介したものと同じです．したがって，ここでは結果を示しません．実行結果は講義の方を参照して下さい．
下に示したLOCATIONとMAPSET名でGRASSを立ち上げます．なお，DATABASEの/home/grassの部分は，前回と同様に各自でファイルを解凍したディレクトリ名に変更してください．また，MAPSETも同じです．
grass5.0beta6
画像を表示します．
d.monコマンドでモニターを開きます．GRASS 5.0beta6 > d.mon x1
d.rastコマンドを用いて，バンド３とバンド４の画像（tmb3，tmb4）を表示してみます．
GRASS 5.0beta6 > d.rast tmb3
Step１．画像の型変換
バンド比の演算結果は実数で取り扱う必要があります．しかし，配布したTM画像は整数型のため，実数型に変換する必要があります．そこで，つぎのように，
r.mapcalcを用いてバンド３（tmb3）とバンド４（tmb4）の画像を実数型に変換します．変換結果はそれぞれ，tmb3.0とtmb4.0とします．GRASS 5.0beta6 >
r.mapcalcStep２．比演算の実行
r.mapcalcコマンドを用いて比演算を行ない，その結果の新しい画像ratioを作成します．ただし，計算結果の値が，小さいので10倍しておきます（これは後で，利用するr.scaleコマンドが，実際の値を整数化したカテゴリー値で分類するためです）．
GRASS 5.0beta6 >
r.mapcalcStep３．カラーテーブルの割り当て
作成した画像ratioに，r.colorsコマンドを用いてカラーテーブルを割り当てます．ここでは，『ryg』を割り当てた例を示します．rygは値の低いものから順に，赤色（r）→黄色（y）→緑色（g）を割り当てます．
GRASS 5.0beta6 >
r.colors>
<Enterキー>d.rastコマンドを用いて，画像ratioを表示します．
GRASS 5.0beta6 >
d.rast ratioStep４．特徴の抽出
比演算の結果から，水域，土壌，植生を区分して抽出した画像を作成します．講義で説明した様に，比率の計算結果は，おおまかに，水は１より小さい値，土壌は１を少し超える値，植物では１を大きく超える数になっていることが予想されます．また，上で行なった実際の計算でも，0.2から4.8の範囲であることがわかったので，ここでは以下のように分類してみました．なお，実際はここで示したよりもより適切な値があると考えられます．
水 ：１以下 （ratioでは，10倍しているので10未満）
土壌：１以上，１.５未満（ratioでは，10倍しているので10以上15未満）
植生：１.５以上 （ratioでは，10倍しているので15以上）
これらの値を用いて，r.rescaleコマンドにより分類します．r.rescaleコマンドで用いたオプションは，つぎの通りです．
r.rescale input=
name [from=min,max] output=name [to=min,max]input=
name ：入力するファイル名水の場合：水の範囲は
ratioでは10未満と想定したので，カテゴリ（整数値）としては，2（値の最小値）〜9に対応します．そこで，出力ファイル名をwaterとし，分類後の値を１として，以下のように実行します．GRASS 5.0beta6 >
r.rescale input=ratio from=2,9 output=water to=1,1土壌の場合：土壌の範囲は
ratioでは10以上15未満と想定したので，カテゴリ（整数値）としては，10〜14に対応します．そこで，出力ファイル名をsoilとし，分類後の値を１として，以下のように実行します．GRASS 5.0beta6 >
r.rescale input=ratio from=10,14 output=soil to=1,1植生の場合：植生の範囲は
ratioでは15以上と想定したので，カテゴリ（整数値）としては，15〜49（値の最大値）に対応します．そこで，出力ファイル名をvegetationとし，分類後の値を１として，以下のように実行します．GRASS 5.0beta6 >
r.rescale input=ratio from=15,49 output=vegetation to=1,1Step５．データの重ね合せ
GRASSでは，多様な方法での重ね合わせが可能です．その１つの方法として，
Hue（色合い）, Intensity（明るさ），およびSaturation（鮮やかさ）を用いたd.hisコマンドがあります．植生と地形との関係を可視化するためにこのd.hisコマンドを用いて，画像の重ね合わせを行なってみましょう．ただし，地形をそのまま表したのでは平面的にはわかりにくいので，地形の傾斜方位（aspect）を利用します．aspectは，地形の効果を表すこのような場合に良く利用されます．aspectは以前にも説明した様に，DEM（地形標高データ）から，r.slope.aspectコマンドを用いて作成します．ただし，ここでは残念ながら全域のDEMがないため，DEMがある領域に限定して作業を進めます．なお，DEMのラスターファイル名はdemで，人工衛星画像と同じ場所に入っています．まず，作業領域の限定ですが，これには，g.regionコマンドを用います．選択する領域の範囲としてdemのラスターファイルがある領域をオプションで指定します．GRASS 5.0beta6 >
g.region raster=demこれにより，
demのある領域に作業が限定されました．つぎに，aspectをr.slope.aspectコマンドを用いて作成します．aspectの出力ファイルをaspとしました．GRASS 5.0beta6 >
r.slope.aspect elevation=dem aspect=aspつぎに，d.hisコマンドを用いて画像を作成します．Hue（色合い）に植生分布であるvegetationを指定し（h_map=vegetation），Intensity（明るさ）にaspを指定します（i_map=asp）．出力ファイルをhveg+iaspという名前にします（out=hveg+iasp）．Saturation（鮮やかさ）は指定しません．
GRASS 5.0beta6 >
d.his h_map=vegetation i_map=asp out=hveg+iaspこれで，ファイルが作成されると同じに画面にも出力されます．同様に，Hue（色合い）を土壌分布であるsoilに変えて，土壌と地形との関係を表す画像を作成してみます．出力ファイルをhsoil+iaspという名前にします．
GRASS 5.0beta6 >
d.his h_map=soil i_map=asp out=hsoil+iaspこのように，重ね合せはリモートセンシングで複数の情報を可視化する際に有効な方法です．
Step１．正規化植生指標の計算
正規化植生指標は，講義で説明したように，つぎの式から計算します．
NDVI＝(Band4-Band3)/(Band4+Band3)
この計算もやはり，
r.mapcalcコマンドを用いてつぎのように求めます．なお，計算結果の出力ファイル名をndviとしました．GRASS 5.0beta6 >
r.mapcalc比演算と同じようにr.colorsコマンドを用いて，この計算結果をrygのカラーテーブルに割り当てます．
GRASS 5.0beta6 >
r.colors map=ndvi color=rygStep２．傾斜方位と正規化植生指標の関係
求めた正規化植生指標が地形の傾斜方位の向き（aspect）とどのような関係にあるかを簡単に計算してみましょう．傾斜方位は，東から度単位で反時計まわりに360°まであります．これを90°ごとのつぎの４つ（東西南北）の方向に分けてみます．
東（East） -- 0〜45 と 316〜360
北（North）-- 46〜135
西（West） -- 136〜225
南（South）-- 226〜315
前に求めた傾斜方位のファイル（
asp）をr.rescaleコマンドを用いて，上の４つにそれぞれ分けます．ただし，東の値は2つの範囲になるので，別々に求め，r.patchコマンドを用いて１つにまとめます．GRASS 5.0beta6 >
r.rescale input=aspect from=0,45 output=east1 to=1,1これで，それぞれ東（
east_slope），北（north_slope），西（west_slope），南（south_slope）の方位のみを示すファイルがそれぞれ作成されます．正規化植生指標（ndvi）がそれぞれの地域でどうなっているかの簡単なレポートをr.reportコマンドを用いて作成します．GRASS 5.0beta6 >
r.report -n map=east_slope,ndvi units=p nsteps=10この結果から，傾斜の方位とNDVIの間に関係があるか考えてみましょう．また，NDVIの最大の値を示す方位を探してみましょう．
また，ここでは説明しませんが，地形の３次元表示（実習５）で説明した
d.3dコマンドを用いて，ndviを地形（dem）に張り合わせて表示してみましょう．
第９回講義に戻る
大阪市立大学広報のホームページ