地球科学におけるGRASS GIS入門 (第6回:10月6日)

 実習です!!


地質モデルの作成

 本講座の第6回目は,前々回入力した地形データを用いて地質モデルを作成してみましょう.地質体は架空のものです.地質関係の方に怒られそうなモデルですが,理論を理解するためのものとしてください.ここでは,3つの地質体から構成されるモデルを作成してみます.


3次元地質モデルの構築

 例として,講義の図−15で示したものと同じモデルで地質モデルを作成してみましょう.必要なものは,講義で述べたように以下の3種類です.

・境界面

・関数 g2 :規則とその実現方法

・関数 g1 :関係コード表とその実現方法

これらに関して,順に説明します.

境界面:境界面S1に関しては,アスキー形式のラスターファイルを準備しましたので,下記のS1.tar.gzをクリックしてダウンロードして,解凍し,ラスターレイヤーを作成してください.データは以前の地形と同じ領域および設定になっています.(前回の実習を行っていない方にはヘッダー部分がそのまま答えになります).GRASSでのラスター名はファイル名と同じく"S1"にして下さい.アスキー形式のラスターを入力するコマンドは,r.in.asciiです.

S1.tar.gz (2KByte)

(ブラウザの設定によりファイル名が変わった場合は,変更してください)

tar -zxvf S1.tar.gz      (解凍方法.S1というファイルが出来上がります)

前々回と異なり,GRASSが設定さえできていればr.in.asciiで入力できる形式になっています.

GRASS 4.2.1 > r.in.ascii ↓  (は,リターンを示す)

 

OPTION: Ascii raster file to be imported

key: input

required: YES

enter option > S1 ↓

 

You have chosen:

input=S1

Is this correct? (y/n) [y]

 

OPTION: Name for resultant raster map

key: output

required: YES

 

Enter raster file name

Enter 'list' for a list of existing raster files

Hit RETURN to cancel request

> S1 ↓

<S1>

 

OPTION: Title for resultant raster map

key: title

format: "phrase"

required: NO

enter option >

 

OPTION: Multiplier for ascii data

key: mult

default: 1.0

required: NO

enter option >

 

CREATING SUPPORT FILES FOR S1

Mapset <map> in Location <sample>

GRASS 4.2.1 >

入力後,d.rastなどで表示して確認して下さい.領域の中心部の標高200mを通る走向EW,傾斜5°Nの平面です.

境界面S2に関しては,高さ180mの水平面として,各自で作成しましょう.作成は,r.mapcalcを用います.次のようにすると作成できます.

GRASS 4.2.1 > r.mapcalc ↓

mapcalc> S2=180 ↓

EXECUTING S2 = ... 100%

CREATING SUPPORT FILES FOR S2

minimum value 180, maximum value 180

mapcalc>

Mapset <map> in Location <sample>

GRASS 4.2.1 >

これも,d.rastで表示してみましょう.平面なので1色で表示されるはずです.

関数g2GRASSの演算機能r.mapcalcを用いてこの関数を実現します.r.mapcalcには以下の様な比較の機能があります.

     if(x,a,b,c)

ここで, x>0 ならば a

     x=0 ならば b

     x<0 ならば c

これは,xの値が正ならばaという値,xの値が0ならばbという値,xの値が負ならばcという値になります.これを用いて関数g2を実現します.ただし,講義ではある地点Pと境界面の上下を比較すると説明しましたが,ここでは,ある地点Pがラスター形式で対象地域全面にあると考え,これを地質情報を得る目的のラスター面So(以降,目的面So)とします.地点Pで考えずにラスター形式の面で考える事により,ある地点Pが境界面のどちらにあるかを判断するのに,各格子点での目的面Soの値と各境界面の値を比較することで可能になります.例えば,ある地点での目的面Soと境界面S1の上下判定は,So-S1を求め,この値が正ならばSoの方が上にあり,0ならば同じ,負ならばSoの方が下にあると判断できます.この計算を実現するのが,ここで述べたr.mapcalcのif(x,a,b,c)です.この「x」の部分を「So-S1」に置き換え,aを1,bを0,cを0とすれば,minsetの各0か,1かが求められます.なお,ここでは,So-S1=0(いいかえれば,SoとS1が等しい)もSoの方が下にあるということにしておきます.

 このままでは,minsetが2進数で得られ,扱いにくいので,普通に扱える10進数に変えることにしましょう.

 2進数を10進数に変えるのは簡単です.例えば,2進数で「100」は,

  ×22×21×20×4+×2+×1=

というように,2進数の下位の桁から順に2の(桁数-1)乗をかけて,それらをたすことにより求められます.これと先ほどのifを組み合わせることにより,10進数のminsetが得られます.例えば,n枚の境界面があり,それぞれの境界面を下位からS1,..Sn,地質情報を得る目的面をSoとします.r.mapcalcでは次のようにします.

 minset=if(So-S1,1,0,0)*2 n-1 + if(So-S2,1,0,0)*2 n-2 +

      .... + if(So-Sn,1,0,0)*2 0

例に示した境界面2枚の場合では,次のようになります.

 mapcalc> minset = if(So-S1,1,0,0)*2 + if(So-S2,1,0,0)*1

なお,ここでは,計算結果の入るラスターをminsetとしています.

関数g1minsetと地質体の関係を表す関係コード表の左側の2進数のminsetも10進数で表すことにします.また,関係コード表の右側の地質体名も番号で表すことにします(地質体b1は"1",b2は"2"...).これにより,図−17で示した関係コード表は次のようになります(赤字の部分).

 0 → 1   m 00 b 1

 1 → 3   m 01 b 3

 2 → 2   m 10 b 2

 3 → 3   m 11 b 3

関数g1は,この表をもとにr.reclassにより再分類して実現します.この表をr.reclassの再分類テーブル形式で書き改めると次のようになります.

    0 = 1

    1 = 3

    2 = 2

    3 = 3

    end

なお,最後の行にはendを付けます.これをviやxeditを用いてファイルとして作成します.ここでは,このファイル名をlogic.mdlとしました.これを用いて,先ほど求めたminsetというファイルをgeologyというファイルに再分類するには,次のようにします.

r.reclass input=minset output=geology < logic.mdl ↓

これで,geologyという地質の番号でできたファイルが完成します.

これを表示すれば目的面Soの地質図を見ることができます.

3次元地質モデルの利用

 目的の表示面を地形面とした実行例を以下に示します.

GRASS 4.2.1 > r.mapcalc ↓   (r.mapcalcコマンドを起動する)

mapcalc> So = dem ↓      (目的の面Soを地形面demからコピーする)

EXECUTING So = ... 100%

CREATING SUPPORT FILES FOR So

minimum value 0, maximum value 654

 

mapcalc> minset = if(So-S1,1,0,0)*2 + if(So-S2,1,0,0) ↓

EXECUTING minset = ... 100%

CREATING SUPPORT FILES FOR minset

minimum value 0, maximum value 3

 

mapcalc>   (演算を終了する(リターン))

 

GRASS 4.2.1 > r.reclass input=minset output=geology < logic.mdl ↓ 

   ここで,inputは入力のファイル名(ここでは,minset),

       outputは再分類結果が入るファイル名(ここでは,geology;地質図),

       < はUnixのリダイレクションで,再分類テーブルとしてlogic.mdlを入力する事を示します.

 Soとして色々な面を作成することにより,多様な面での地質図が作成できます.例えば,標高100mの水平断面図が作成したければ,r.mapcalcの最初にSo=demとしたところを,So=100とします.

 以上で,3次元地質モデルを用いて目的とした面での地質図(geology)が完成します.このラスターファイルをd.rastやd.3dを用いて表現してみましょう.

 なお,ここで示したモデルでは,地表上の空気の部分を考慮に入れていないため,空気の部分も地質体が表示されてしまいます.いわば,侵食される前の地質図??を表示することになります.これを改良するためには,空気は最も新しい地質体として考えることにより,対応できます.各自で考えてみてください.また,S2を180mとしましたが,これも変えてみるのも面白いです(レポートではありませんが,わかった方はMLに報告してください).


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


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

第6回講義に戻る

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