地球科学におけるGRASS GIS入門 (第5回:8月25日)

 実習です!!


GRASSの操作

 本講座の第5回目は,前回入力した地形データで簡単な計算を行ってみます.使うコマンドは基本的に1つだけです.色々試してみましょう.なぜか,講座の方がマニュアルのようで,実習のほうが本文の様な感じになってきました.


ラスター画像の計算(r.mapcalc)

 GRASSでラスターおよびラスター間の計算を行うコマンドとしてr.mapcalcがあります.これを用いると色々な簡単な計算が行え,新しいラスターデータを作成することができます.

基本式

 計算式の基本形は,下記に示す様な代入形式です.

   結果 = 式

結果:ラスターのファイル名です.計算実行後には,ラスターファイルのみではなく,サポートファイル等も作成されます.

:ラスターのファイル名と定数,およびつぎに示す演算子・関数が利用できます.ラスターのファイル名は,直接変数名となります.

   例えば,dem の高さを2倍して dem2 というラスターを作成する場合は,次の通りです.

     dem2 = dem * 2

 なお,ラスターのファイル名に数字や,演算子・関数名を含んだ名前をつけていた場合は""で囲みます.

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

     c = "500" * "a-b"

となります.

演算子

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

演算子

意 味

形  式

説  明

 %

剰 余

算術演算

割り算の余り

 /

徐 算

算術演算

割り算(÷)

 *

乗 算

算術演算

かけ算(×)

 +

加 算

算術演算

足し算(+)

 -

減 算

算術演算

引き算(−)

       

 ==

等しい

関係演算

左と右の値が等しい

 !=

等しくない

関係演算

左と右の値が等しくない

 >

大きい

関係演算

左の値が右の値より大きい

 >=

大きいか等しい

関係演算

左の値が右の値より大きいか等しい

 <

小さい

関係演算

左の値が右の値より小さい

 <=

小さいか等しい

関係演算

左の値が右の値より小さいか等しい

       

 &&

論理積(and)

論理演算

左 かつ 右

 ||

論理和(or)

論理演算

左 または 右

注1)関係演算の結果は,式が正しい場合は真(1),そうでなければ偽(0)となります.

  例: a = b > c  bの値がcより大きい場合は a は 1,それ以外は 0 となります.

注2)論理演算の用いる値は,基本的に真(1),か偽(0)として下さい.これらは関係演算子を結ぶときにつかいます

  例: a = ( 10 <= b ) && ( b <= 100 ) b の値が10以上で,かつ100以下の場合は,a は 1,そうでなければ 0 となります.

注3)これらは,基本的にC言語の演算子と同じです.

関 数

r.mapcalcで利用できる関数を以下に示します.

 関 数

説  明

 abs(x)

絶対値 |x|

 exp(x)

指数

 exp(x,y)

x の y 乗

 log(x)

自然対数

 log(x,b)

底 b の対数

 sqrt(x)

√x(平方根)

 float(x)

整数を実数へ変換

 int(x)

整数化(切り捨て)

 round(x)

整数化(四捨五入)

 sin(x)

正弦(サイン)

 cos(x)

余弦(コサイン)

 tan(x)

正接(タンジェント)

 atan(x)

逆正接(アークタンジェント)

   

 max(x,y,z,…)

x,y,z,…の中の最大値

 min(x,y,z,…)

x,y,z,…の中の最小値

 median(x,y,z,…)

x,y,z,…の中のメディアン(中央値)

 rand(x,y)

x〜yの間の乱数

   

 if(x)

条件判断,もし x が 0 でなければ1,他は0

 if(x,a)

もし x が 0 でなければ a,他は0

 if(x,a,b)

もし x が 0 でなければ a,他は b

 if(x,a,b,c)

もし x > 0 ならば a,x = 0 ならば b,x < 0 ならば c

 

分解能

 r.mapcalcは,現在設定されている分解能で常に計算を行い,結果のラスターを作成します.そのため,異なる分解能のラスターデータ間でも,演算を行うことができます.

 

近傍表現

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

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

これを用いることにより,画像処理等でよく行う,マトリックスオペレータ的な処理が行えます(コマンドr.mfilterでも可能).例えば,demの各点とその周囲の8近傍を平均して,その点の値にし,mdemという新しいラスターを作成するには,次のような式になります.

    mdem = ( 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

この計算により,凹凸を少し減らした滑らかな地形を計算で求めることができます.


 では,実際にGRASSを起動して,前回作成した新しい地域の地形データdemを用いてr.mapcalcで計算し,簡単な地形の解析を行ってみましょう.また,計算結果を画面に表示して,考え方が正しく反映しているか見て検討して下さい.

実行方法の例をリストに示します(<-クリックしてください).

地形面の平滑化:地形面の小さな凹凸を減らして滑らかにします.各点とその周辺の8近傍の値を足して,その数(9)で割ります.これは平均を行った先ほどの例をそのまま利用します.計算結果のラスターデータをm1demとします.

 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

 平滑化の断面イメージ

平滑化が強すぎて,地形の特徴がわかりにくくなる場合は,原地形に重みを付けて平均を取ります(加重平均).例えば,評価点の影響を2倍にするには,次のようにします.計算結果のラスターデータをm2demとします.

 m2dem = ( dem[-1,-1]+dem[0,-1]+dem[ 1,-1]+dem[-1,0]+2*dem[0,0]+dem[ 1,0]+dem[-1, 1]+dem[0, 1]+dem[ 1, 1] ) / 10

このほかに,中央値を利用する,次のもなども考えられます.計算結果のラスターデータをmeddemとします.

 meddem = mdeian ( 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] )

接峰面図の作成:地形面の小さな凹(谷)を減らして接峰面を作成します.各点とその周辺の8近傍の値の中で最も高い値をその点(評価点)の値とします.計算結果のラスターデータを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] )

 接峰面の断面イメージ

接谷面図の作成:各点とその周辺の8近傍の値の中で最も低い値をその点(評価点)の値とします.計算結果のラスターデータをdemminとします.

 demmin = min ( 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] )

 接谷面の断面イメージ

起伏量図の作成:各点とその周辺の8近傍の値の中で最も高い値と低い値の差を起伏とします.すでに,最高値と最低値は求めてあるので,次のように簡単に求められます.計算結果のラスターデータをdemsaとします.

 demsa = demmax - demmin

 起伏量図の断面イメージ

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

   

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

   

となります.これをデジタルラプラシアンといいます.具体的には,評価点と左右の点の差の差と評価点と上下の点の差の差の和となり,書きかえると上下左右の点の値の和と評価点を4倍した値の差となります.

r.mapcalcでこれを計算し,結果のラスターデータをdemlap4とします.

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

さらに,斜め方向も考えに入れて(8近傍)計算し,結果のラスターデータをdemlap8とします.

 demlap8 = ( dem[-1,-1]+dem[0,-1]+dem[ 1,-1]+dem[-1,0]+dem[ 1,0]+dem[-1, 1]+dem[0, 1]+dem[ 1, 1] ) - 8 * dem[0,0]

このような方法でラプラシアンが求められます.値の絶対値が大きいとことがエッジです.このため,絶対値をそれぞれとって,dmelap4aとdmelap8aも求めてみましょう.

 demlap4a = abs ( demlap4 )

 demlap8a = abs ( demlap8 )

地形面の先鋭化:画像を先鋭化(sharpening)するために,原画像からラプラシアンを引く方法がよく用いられています(説明は省略します).地形データにこの考えを導入して,地形面の先鋭化を行って見ましょう.結果のラスターデータを4近傍,8近傍でそれぞれdemshp4,demshp8とします.

 demshp4 = dem - demlap4

 demshp8 = dem - demlap8

 

 以上の結果を,d.rastやd.3dを用いて2次元や3次元で表示し,どの様な結果になったかを検討してください.ここでは,3×3の9個の点についての操作でしたが,5×5やもっと多くの周辺を用いて計算することが必要な場合もあります.この他にも,色々な計算が考えられます.面白いあるいは実用的な計算方法を考え,試して見て下さい.

 なお,これらの計算式をそのまま実行すると,データのいちばん外側の部分はデータがないために値ゼロが用いられます.そのため,計算で求めたラスター画像の周囲は変な結果になります.r.mapcalcで計算後,表示する際にg.regionを用いて領域を少し狭めるなどの工夫が必要です.また,色の設定等も工夫してみてください.

注意:今回作成したファイルやMAPSETは今後も利用しますので,終了時に削除しないで下さい!!


 以上で,今回の実習は終了です.

第5回講義に戻る

大阪市立大学インターネット講座'99