実習です！！ 実習です！！本講座の第６回目は，前々回入力した地形データを用いて地質モデルを作成してみましょう．地質体は架空のものです．地質関係の方に怒られそうなモデルですが，理論を理解するためのものとしてください．ここでは，３つの地質体から構成されるモデルを作成してみます．
例として，講義の図−１５で示したものと同じモデルで地質モデルを作成してみましょう．必要なものは，講義で述べたように以下の３種類です．
・境界面
・関数
g2 ：規則とその実現方法・関数
g1 ：関係コード表とその実現方法これらに関して，順に説明します．
境界面：
境界面S1に関しては，アスキー形式のラスターファイルを準備しましたので，下記のS1.tar.gzをクリックしてダウンロードして，解凍し，ラスターレイヤーを作成してください．データは以前の地形と同じ領域および設定になっています．（前回の実習を行っていない方にはヘッダー部分がそのまま答えになります）．GRASSでのラスター名はファイル名と同じく"S1"にして下さい．アスキー形式のラスターを入力するコマンドは，r.in.asciiです． S1.tar.gz （２KByte）(ブラウザの設定によりファイル名が変わった場合は，変更してください)
tar -zxvf S1.tar.gz （解凍方法．S1というファイルが出来上がります）
前々回と異なり，GRASSが設定さえできていればr.in.asciiで入力できる形式になっています．
GRASS 4.2.1 >
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などで表示して確認して下さい．領域の中心部の標高２００ｍを通る走向EW，傾斜５°Nの平面です．
境界面S2に関しては，高さ180ｍの水平面として，各自で作成しましょう．作成は，r.mapcalcを用います．次のようにすると作成できます．
GRASS 4.2.1 >
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で表示してみましょう．平面なので１色で表示されるはずです．
関数
g2：GRASSの演算機能r.mapcalcを用いてこの関数を実現します．r.mapcalcには以下の様な比較の機能があります．ｉｆ（ｘ，ａ，ｂ，ｃ）
ここで， ｘ＞０ ならば ａ
ｘ＝０ ならば ｂ
ｘ＜０ ならば ｃ
これは，ｘの値が正ならばａという値，ｘの値が0ならばｂという値，ｘの値が負ならばｃという値になります．これを用いて関数
g2を実現します．ただし，講義ではある地点Ｐと境界面の上下を比較すると説明しましたが，ここでは，ある地点Ｐがラスター形式で対象地域全面にあると考え，これを地質情報を得る目的のラスター面So（以降，目的面So）とします．地点Ｐで考えずにラスター形式の面で考える事により，ある地点Ｐが境界面のどちらにあるかを判断するのに，各格子点での目的面Soの値と各境界面の値を比較することで可能になります．例えば，ある地点での目的面Soと境界面S1の上下判定は，So-S1を求め，この値が正ならばSoの方が上にあり，0ならば同じ，負ならばSoの方が下にあると判断できます．この計算を実現するのが，ここで述べたr.mapcalcのif(x,a,b,c)です．この「x」の部分を「So-S1」に置き換え，ａを１，ｂを0，ｃを0とすれば，minsetの各0か，1かが求められます．なお，ここでは，So-S1=0（いいかえれば，SoとS1が等しい）もSoの方が下にあるということにしておきます．このままでは，minsetが２進数で得られ，扱いにくいので，普通に扱える１０進数に変えることにしましょう．
２進数を１０進数に変えるのは簡単です．例えば，２進数で「１００」は，
１×２2＋０×２1＋０×２0＝１×４＋０×２＋０×１＝４
というように，２進数の下位の桁から順に２の（桁数-１）乗をかけて，それらをたすことにより求められます．これと先ほどのifを組み合わせることにより，１０進数の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
例に示した境界面２枚の場合では，次のようになります．
mapcalc> minset = if(So-S1,1,0,0)*2 + if(So-S2,1,0,0)*1
なお，ここでは，計算結果の入るラスターをminsetとしています．
関数
g1：minsetと地質体の関係を表す関係コード表の左側の２進数のminsetも１０進数で表すことにします．また，関係コード表の右側の地質体名も番号で表すことにします（地質体b1は"１"，b2は"２"．．．）．これにより，図−１７で示した関係コード表は次のようになります（赤字の部分）．０ → １ （m 00 → b 1）
１ → ３ （m 01 → b 3）
２ → ２ （m 10 → b 2）
３ → ３ （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の地質図を見ることができます．
目的の表示面を地形面とした実行例を以下に示します．
GRASS 4.2.1 >
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として色々な面を作成することにより，多様な面での地質図が作成できます．例えば，標高100ｍの水平断面図が作成したければ，r.mapcalcの最初にSo=demとしたところを，So=100とします．
以上で，３次元地質モデルを用いて目的とした面での地質図（geology）が完成します．このラスターファイルをd.rastやd.3dを用いて表現してみましょう．
なお，ここで示したモデルでは，地表上の空気の部分を考慮に入れていないため，空気の部分も地質体が表示されてしまいます．いわば，侵食される前の地質図？？を表示することになります．これを改良するためには，空気は最も新しい地質体として考えることにより，対応できます．各自で考えてみてください．また，S2を180ｍとしましたが，これも変えてみるのも面白いです（レポートではありませんが，わかった方はMLに報告してください）．
注意：今回作成したファイルやMAPSETは今後も利用しますので，終了時に削除しないで下さい！！
