GRASSを用いた地理情報システム入門（第７回）

第７回目

　第７回目の講義です．今回はGRASSのラスター演算の基礎について説明します．

　実習では，講義で説明したラスター演算を地形データを用いて実際に解析してみます．

４．解析の基礎と応用

４-１．ラスターの演算

４-１-１．r.mapcalcを用いたラスター演算

　前述したように，GRASSはラスター型のGISが基本であるため，多くの解析はラスター地図を用いて行なうのが一般的です．そのため，ラスター演算に関しては多様なコマンドが用意されています．これらの中で，自分で計算式を設定して，ラスターおよびラスター間の計算を行うコマンドr.mapcalcがあります．これを用いると色々な計算が簡単に行え，新しいラスター地図を容易に作成することができます．ここでは，このr.mapcalcコマンドについて説明します．

［基 本 式］

　r.mapcalcの計算式の基本形は，つぎにに示すような代入形式です．

　例えば，abcというラスター地図があったとします．このラスター地図の各セルの値を２倍して def というラスター地図を作成する場合は，次の通りです．

　なお，ラスター地図のファイル名に数字や，演算子・関数名を含んだ名前をつけていた場合は，そのままではファイルとして識別できないので，"　"で囲みます．

　例えば，500と言うファイル名とa-bというファイル名を掛け合わせ，cという結果を得る場合は，

となります．

［演 算 子］

　r.mapcalcで利用できる演算子を以下に示します．

　　例えば：

では，bの値がcより大きい場合は a は 1，それ以外は 0 となります．

　　例えば：

ではbの値が10以上で，かつ100以下の場合は，a は 1，そうでなければ 0 となります．

注３）これらは，基本的にC言語の演算子と同じです．

［関　数］

r.mapcalcで利用できる関数をつぎに示します．

［分 解 能］

　r.mapcalcは，現在設定されている分解能で計算を行い，結果のラスターを作成します．そのため，異なる分解能のラスターデータ間でも，演算を行うことができます．なお，基本的にGRASSではr.mapcalcコマンドに限らず，ほとんどのラスター地図の計算は，現在設定されている分解能（g.regionコマンド等で設定・確認可能）で結果の出力を行ないます．ただし，計算結果は元のデータの分解能に依存するため，計算時の分解能を高くしても，実際の分解能が高くなるわけではありませんので，注意してください．

［近傍表現］

　r.mapcalcでは，ラスターのマトリックスデータの各セル（点；評価点）の周辺（近傍）を配列[j，i]を用いて相対的にオフセットで指定して，利用することができます．j は縦方向下向きに＋で行数，i は横方向右向きに＋で桁数です．

　例えば，dem というファイルの各セルの１つ右下のデータは dem[１，１]です．また，周辺の各セルは下図のようになります（負の方向は−を使いますが，正の方向に＋は利用できません.-1は使えますが，+1ではなくて1としてください）．

　これを用いることにより，画像処理等でよく行う，マトリックスオペレータ的な処理が行えます（コマンドr.mfilterでも可能）．

　例１：aというラスター地図を１セル分だけ東方向に移動した（ずらした）bという新しいラスターを作成するには，次のような式になります．

　例２：cというラスター地図の各セルとその周囲の４近傍のセルの値を平均（全て（5つ）を足し合わせて，セルの個数（5）で割る）して，そのセルの値にし，dという新しいラスターを作成するには，次のような式になります．

　なお，これらの計算式を実行する場合，領域の端の部分にデータがない場合が出てきます．その場合は，値がnullとして入ります．計算結果を解釈する場合や，さらに計算を行う場合，また表示する際には，このことをよく注意しておいて下さい．

４-１-２．地形への簡単な応用

　ラスター演算の応用例として，地形の簡単な解析を説明します．なお，ここでは計算の概念と計算方法のみを示します．実際の計算結果は実習を見て下さい．

［平滑化］

　地形の小さな凹凸を減らして滑らかにします．いろいろな方法がありますが，ここでは各セルとその周囲（近傍）セルの高さデータを平均して，そのセルの値にします．その概念を断面図で下に示します．

　例えば，地形のラスターファイル名をdemとします．各セルとその周辺の８近傍のセルの値（標高値）を足して，そのセル数（９）で割ります．これは近傍表現の説明時に例として用いたものを応用したものです．計算結果のラスター地図名をm1demとし，r.mapcalcの式をつぎに示します．

　m1dem = ( dem[-1,-1]+dem[0,-1]+dem[ 1,-1]+dem[-1,0]+dem[0,0]+dem[ 1,0]+dem[-1, 1]+dem[0, 1]+dem[ 1, 1] ) / 9

　この他に，原地形に重みを付けて平均を取る方法（加重平均）や中央値（メディアン）をとる方法などがあります．

［接峰面図の作成］

　接峰面図とは，峰（尾根）に接する地形面を表したもので，地形の上に薄い布をかぶせて谷の部分を取り除いたような仮想の地形を表現します．実際には，ある一定範囲内で地形の最高点をとりその情報から再度地形図を作成します．例えば，セルの大きさで３×３の領域を考えます．各セルとその周辺の８近傍のセルの標高値の中で最も高い値をそのセル（評価点）の値として求めます．これにより，接峰面が作成できます．

　具体的にこの計算をr.mapcalcで行うには，max関数を用います．元の地形にラスターファイル名をdemとすると，次式で求められます．接峰面の計算結果のラスタ地図demmaxとします．

　demmax = max ( dem[-1,-1],dem[0,-1],dem[ 1,-1],dem[-1,0],dem[0,0],dem[ 1,0],dem[-1, 1],dem[0, 1],dem[ 1, 1] )

［エッジの検出］

　画像処理の手法ですが，濃度（ここでは高度）が急激に変化して１つの領域が終わり他の領域が始まっているのを示す場所をエッジ（edge）と言います．これはある種の不連続性（discontinutiy）を検出することになります．色々な考え方があるのですが，ここでは２次の偏微分であるラプラシアン（Laplacian）を示します．ラプラシアンは次の通りです．

これをデジタルデータにあてはめるための離散近似は，

となります．これをデジタルラプラシアンといいます．

具体的には，評価点と左右の点の差の差と評価点と上下の点の差の差の和となり，書きかえると上下左右の点の値の和と評価点を４倍した値の差となります．

元の地形をdemとすると，r.mapcalcでは次式のように表せます．結果のラスター地図をdemlap4とします．

　demlap4 = ( dem[ 1,0] + dem[-1,0] + dem[0, 1] + dem[0,-1] ) - 4 * dem[0,0]

このような方法でラプラシアンが求められます．値の絶対値が大きいところがエッジです．絶対値をとるには，abs関数を用いてつぎの式で行ないます．計算結果はdmelap4aとします．

　demlap4a = abs ( demlap4 )

［地形面の鮮鋭化］

　画像を先鋭化（sharpening）するために，原画像からラプラシアンを引く方法がよく用いられています（説明は省略します）．地形データにこの考えを導入して，地形の先鋭化を行うには，地形データから上述したラプラシアンを引き算します．これは次式のように求めます．結果のラスターデータをdemshp4とします．

　demshp4 = dem - demlap4

以上のように，r.mapcalcを用いることによりラスター地図の簡単な演算が行なえます．なお，これらの例は実習の方に示しました．

まとめ（７）

（以上，第７回講義）