GRASSを用いた地理情報システム入門(第7回)


目    次

4.解析の基礎と応用
  4-1.ラスターの演算
   4-1-1.r.mapcalcを用いたラスター演算
   4-1-2.地形への簡単な応用
まとめ(7)


第7回目

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

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


4.解析の基礎と応用

4-1.ラスターの演算

4-1-1.r.mapcalcを用いたラスター演算

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

[基 本 式]

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

   結果 = 式

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

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

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

     def = abc * 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)

逆正接(アークタンジェント)(範囲は-90°〜90°)

 atan(x,y)

逆正接(アークタンジェント)(範囲は0°〜360°)

 

 

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

[近傍表現]

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

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

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

 例1:aというラスター地図を1セル分だけ東方向に移動した(ずらした)bという新しいラスターを作成するには,次のような式になります.

    b = a[0,1]

 例2:cというラスター地図の各セルとその周囲の4近傍のセルの値を平均(全て(5つ)を足し合わせて,セルの個数(5)で割る)して,そのセルの値にし,dという新しいラスターを作成するには,次のような式になります.

    d = (c[-1,0]+c[0,0]+c[0,-1]+c[0,1]+c[1,0])/5

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

 

4-1-2.地形への簡単な応用

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

[平滑化]

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

 

 例えば,地形のラスターファイル名をdemとします.各セルとその周辺の8近傍のセルの値(標高値)を足して,そのセル数(9)で割ります.これは近傍表現の説明時に例として用いたものを応用したものです.計算結果のラスター地図名を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

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

[接峰面図の作成]

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

 

 具体的にこの計算を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] )

[エッジの検出]

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

 

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

 

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

 

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

元の地形を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を用いることによりラスター地図の簡単な演算が行なえます.なお,これらの例は実習の方に示しました.


まとめ(7)

ラスター演算(r.mapcalc)

基 本 式: 結果 = 式

演算子・関数: 多様な演算子と関数

分解能: 現在設定されている分解能で演算

近傍表現: セルで近傍を表現したものが演算に利用可能

(以上,第7回講義)


 

第7回 実習のページへ

大阪市立大学広報ホームページ