実習です！！ 実習です！！本講座の第５回目は，前回入力した地形データで簡単な計算を行ってみます．使うコマンドは基本的に１つだけです．色々試してみましょう．なぜか，講座の方がマニュアルのようで，実習のほうが本文の様な感じになってきました．
GRASSでラスターおよびラスター間の計算を行うコマンドとしてr.mapcalcがあります．これを用いると色々な簡単な計算が行え，新しいラスターデータを作成することができます．
計算式の基本形は，下記に示す様な代入形式です．
結果 ＝ 式
結果
：ラスターのファイル名です．計算実行後には，ラスターファイルのみではなく，サポートファイル等も作成されます．式
：ラスターのファイル名と定数，およびつぎに示す演算子・関数が利用できます．ラスターのファイル名は，直接変数名となります．例えば，dem の高さを２倍して dem2 というラスターを作成する場合は，次の通りです．
dem2 = dem * 2
なお，ラスターのファイル名に数字や，演算子・関数名を含んだ名前をつけていた場合は""で囲みます．
例えば，500と言うファイル名とa-bというファイル名を掛け合わせ，cという結果を得る場合は，
c = "500" * "a-b"
となります．
r.mapcalcで利用できる演算子を以下に示します．
|
演算子 |
意 味 |
形 式 |
説 明 |
|
% |
剰 余 |
算術演算 |
割り算の余り |
|
/ |
徐 算 |
算術演算 |
割り算（÷） |
|
* |
乗 算 |
算術演算 |
かけ算（×） |
|
+ |
加 算 |
算術演算 |
足し算（＋） |
|
- |
減 算 |
算術演算 |
引き算（−） |
|
== |
等しい |
関係演算 |
左と右の値が等しい |
|
!= |
等しくない |
関係演算 |
左と右の値が等しくない |
|
> |
大きい |
関係演算 |
左の値が右の値より大きい |
|
>= |
大きいか等しい |
関係演算 |
左の値が右の値より大きいか等しい |
|
< |
小さい |
関係演算 |
左の値が右の値より小さい |
|
<= |
小さいか等しい |
関係演算 |
左の値が右の値より小さいか等しい |
|
&& |
論理積(and) |
論理演算 |
左 かつ 右 |
|
|| |
論理和(or) |
論理演算 |
左 または 右 |
注１）関係演算の結果は，式が正しい場合は真（１），そうでなければ偽（０）となります．
例：
a = b > c bの値がcより大きい場合は a は 1，それ以外は 0 となります．注２）論理演算の用いる値は，基本的に真（１），か偽（０）として下さい．これらは関係演算子を結ぶときにつかいます
例：
a = ( 10 <= b ) && ( b <= 100 ) b の値が10以上で，かつ100以下の場合は，a は 1，そうでなければ 0 となります．注３）これらは，基本的に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 でなければ１，他は０ |
|
if(x,a) |
もし x が 0 でなければ a，他は０ |
|
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 というファイルの各点の１つ右下のデータは dem[１，１]です．また，周辺の各点は下図のようになります（負の方向は−を使いますが，正の方向に＋は利用できません.-1は使えますが，+1ではなくて1としてください）．
これを用いることにより，画像処理等でよく行う，マトリックスオペレータ的な処理が行えます（コマンドr.mfilterでも可能）．例えば，demの各点とその周囲の８近傍を平均して，その点の値にし，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で計算し，簡単な地形の解析を行ってみましょう．また，計算結果を画面に表示して，考え方が正しく反映しているか見て検討して下さい．
実行方法の例をリストに示します
（<-クリックしてください）．
地形面の平滑化：
地形面の小さな凹凸を減らして滑らかにします．各点とその周辺の８近傍の値を足して，その数（９）で割ります．これは平均を行った先ほどの例をそのまま利用します．計算結果のラスターデータを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
平滑化が強すぎて，地形の特徴がわかりにくくなる場合は，原地形に重みを付けて平均を取ります（加重平均）．例えば，評価点の影響を２倍にするには，次のようにします．計算結果のラスターデータを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] )
接峰面図の作成：
地形面の小さな凹（谷）を減らして接峰面を作成します．各点とその周辺の８近傍の値の中で最も高い値をその点（評価点）の値とします．計算結果のラスターデータを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] )
接谷面図の作成：
各点とその周辺の８近傍の値の中で最も低い値をその点（評価点）の値とします．計算結果のラスターデータを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] )
起伏量図の作成：
各点とその周辺の８近傍の値の中で最も高い値と低い値の差を起伏とします．すでに，最高値と最低値は求めてあるので，次のように簡単に求められます．計算結果のラスターデータをdemsaとします．demsa = demmax - demmin
エッジの検出：
画像処理の手法ですが，濃度（ここでは高度）が急激に変化して１つの領域が終わり他の領域が始まっているのを示す場所をエッジ（edge）と言います．これはある種の不連続性（discontinutiy）を検出することになります．色々な考え方があるのですが，ここでは２次の偏微分であるラプラシアン（Laplacian）を示します．ラプラシアンは次の通りです． 
これをデジタルデータにあてはめるための離散近似は，
となります．これをデジタルラプラシアンといいます．具体的には，評価点と左右の点の差の差と評価点と上下の点の差の差の和となり，書きかえると上下左右の点の値の和と評価点を４倍した値の差となります．
r.mapcalcでこれを計算し，結果のラスターデータをdemlap4とします．
demlap4 = ( dem[ 1,0] + dem[-1,0] + dem[0, 1] + dem[0,-1] ) - 4 * dem[0,0]
さらに，斜め方向も考えに入れて（８近傍）計算し，結果のラスターデータを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）するために，原画像からラプラシアンを引く方法がよく用いられています（説明は省略します）．地形データにこの考えを導入して，地形面の先鋭化を行って見ましょう．結果のラスターデータを４近傍，８近傍でそれぞれdemshp4，demshp8とします．demshp4 = dem - demlap4
demshp8 = dem - demlap8
以上の結果を，d.rastやd.3dを用いて２次元や３次元で表示し，どの様な結果になったかを検討してください．ここでは，３×３の９個の点についての操作でしたが，５×５やもっと多くの周辺を用いて計算することが必要な場合もあります．この他にも，色々な計算が考えられます．面白いあるいは実用的な計算方法を考え，試して見て下さい．
なお，これらの計算式をそのまま実行すると，データのいちばん外側の部分はデータがないために値ゼロが用いられます．そのため，計算で求めたラスター画像の周囲は変な結果になります．r.mapcalcで計算後，表示する際にg.regionを用いて領域を少し狭めるなどの工夫が必要です．また，色の設定等も工夫してみてください．
注意：今回作成したファイルやMAPSETは今後も利用しますので，終了時に削除しないで下さい！！
