GRASSを用いた地理情報システム入門

 第7回実習です!


 第7回目は,講義で説明した地形の解析を前々回の地形データを用いて実際に行なってみましょう.また,求めた結果を2次元や3次元で表示してみましょう.


地形を用いたラスター演算

[平滑化]

 講義で説明した原理に基づいて,実際の地形を平均して滑らかにしてみます.元の地形のラスターファイル名はdemです.まず,GRASSを立ち上げ,テストデータとなる地形のラスター地図であるdemを2次元と3次元で表示してみましょう.2次元表示のコマンドはd.rastで3次元表示のコマンドはd.3dです.

  地形図demの2次元表示例(クリックすると拡大します)

  地形図demの3次元表示例(鉛直方向に5倍)(クリックすると拡大します)

また,凡例を表示する代わりに,ヒストグラムを表示してみましょう.ヒストグラムを表示するコマンドは,d.histogramです.

  地形図demのヒストグラム(クリックすると拡大します)

 これらの図は,地形を地形らしく見えるようにカラーテーブルを別途作成しています.このカラーテーブルが必要な場合はつぎをクリックしてください.

  カラーテーブル

 平滑化を行なうために,ラスター演算コマンドr.mapcalcを起動します.

GRASS 5.0beta6 > r.mapcalc

プロンプトがつぎのように変わります.

mapcalc>

このプロンプトに続いて,式を入力します(r.mapcalcコマンドの終了する場合は<Enter>を入力する).各セルとその周辺の8近傍のセルの値(標高値)の平均を求めるには,講義で説明したように次式を入力します.計算結果を出力するラスター地図名をm1demとしました.

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

(この画面では,改行されて2行に表示されていますが,実際には続けて入力してください.)

 計算結果をd.rastコマンドで表示して比較してみましょう.

  地形の平均の結果(クリックすると拡大します)

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

mapcalc> 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とします.

mapcalc> meddem=median(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とするとr.mapcalcの計算式はつぎのようになります.max()は()内の値の最高値を出力する関数です.

mapcalc> 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])

 接峰面図では最高点を用いましたが,最低点を採用すると接谷面図というものが作成できます.これは,上記の式の関数maxの代わりに,最小値を求める関数min()を用いることで得られます.計算結果の接谷面をdemminとすると,つぎのようになります.

mapcalc> 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])

 接峰面,接谷面のそれぞれの計算結果を表示して,もとの地形図と比べてみましょう.

  接峰面(demmax)の表示例(クリックすると拡大します)

  接谷面(demmin)の表示例(クリックすると拡大します)

[起伏量図の作成]

 地形に凹凸の度合いを示すものの1つに起伏量図というものがあります.これは,ある領域内の標高の最高点と最低点を求めその差を求めることにより得られます.ここでは,各セルととその周辺の8近傍の標高値の中で最も高い値と最も低い値の差を起伏とします.すでに,接峰面図と接谷面図作成の際に最高値と最低値は求めてあるので前述のラスター地図をそのまま用いて,次のように簡単に求められます.計算結果の起伏量図をdemsaとします.

mapcalc> demsa = demmax - demmin

  起伏量図(demsa)の表示例(クリックすると拡大します)

[エッジの検出]

 講義で説明したラプラシアンをr.mapcalcで求めて,エッジの検出を行なってみましょう.ラプラシアンの計算結果のラスター地図をdemlap4とします.

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

  ラプラシアンの計算結果(4近傍;demlap4)の表示例(クリックすると拡大します)

さらに,斜め方向も考慮した8近傍を計算してみます.結果のラスター地図をdemlap8とします.

mapcalc> 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]

  ラプラシアンの計算結果(8近傍;demlap8)の表示例(クリックすると拡大します)

 エッジを求めるために4近傍の結果の絶対値を求めてみます.その結果をdmelap4aとします.

mapcalc> demlap4a = abs ( demlap4 )

  エッジの検出結果(4近傍;demlap4a)の表示例(クリックすると拡大します)

[地形面の鮮鋭化]

 上記で求めたラプラシアンを用いて地形面の先鋭化を行なってみましょう.結果のラスター地図を4近傍,8近傍でそれぞれdemshp4,demshp8とします.

mapcalc> demshp4 = dem - demlap4

mapcalc> demshp8 = dem - demlap8

計算結果を表示して,もとの地形とどのように変わったかを比較してみましょう.

  鮮鋭化の結果(4近傍;demshp4)の表示例(クリックすると拡大します)

  鮮鋭化の結果(8近傍;demshp8)の表示例(クリックすると拡大します)

[if の応用]

 if()を応用すると,色々な地形の操作ができます.例えば,対象地域の標高150m以上を削って,平らにすることや,逆に150m以下を全て埋めてしまうことなどが簡単にできます.if()はつぎの通りです.

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

例1:標高150m以上を削る.標高が150m以上であるかどうかの判断をするために,(標高-150m)を計算します.この差が,正あるいは0ならば150m以上と判断できるので,その結果は150mとします.負ならば,もとの地形の標高を代入します.その計算式をつぎに示します.計算結果のラスター地図をdemifa150とします.

mapcalc> demifa150 = if ( dem - 150, 150, 150, dem )

  ifの応用例1の結果の2次元表示(150m以上を削る)(クリックすると拡大します)

   同 3次元表示(クリックすると拡大します)

例2:標高150m以下を埋める.例1と同様な手順で判断します.その計算式をつぎに示します.計算結果のラスター地図をdemifb150とします.

mapcalc> demifb150 = if ( dem - 150, dem, 150, 150 )

  ifの応用例2の結果の2次元表示(150m以下を埋める)(クリックすると拡大します)

   同 3次元表示(クリックすると拡大します)

 ここの例では,この演算を地域全域を対象としましたが,MASKを利用することにより,特定の地域のみ(例えば,1つの山や1つの谷など)を指定することができます.

 


 第7回講義に戻る

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