数値地図25000（地図画像）の入力

　国土地理院の刊行している数値地図25000は2万5千分の1地形図を画像化したものです．１ピクセルが0.1×0.1mm（2.5m×2.5m）のTiff形式で保存されており，1ピクセル８ビットで各ビットにそれぞれ意味がありビットごとに分けて使うと便利です．なお，下に示したように地形図そのものの画像であるため，余白や外側の線，文字等がはいっています．また，ＵＴＭ投影のため，画像は長方形ですが，地形図は下図のように台形になっています．これらの情報は，CD-ROMのDATAディレクトリの下のKANRI.CSVファイルの中に記述されています．とくに赤い字で示した４隅の情報（地図画像の４隅と実際の緯度・経度）は位置合わせ時に重要です．

（１）ｘｙ座標系での入力

　画像入力用の中間ファイル等を作成するためのｘｙ座標系をとる別のLOCATIONとMAPSETをGRASSで作成します．領域は画像より少し大きめの横5000，縦4000程度で分解能は縦横とも１にします．

　r.in.tiffコマンドによりTiff画像をラスター形式でGRASSに入力します．

　GRASS 5.0beta6 > r.in.tiff
　OPTION: Name on TIFF file to input.
　　　 key: input
　required: YES
　enter option > /cdrom/data/5538/553867.tif

　　（地図画像のファイル名をフルパスで入力，/cdromは設定により異なります）

　　----- 省 略（以降，入力の確認は省略します） -----

　Enter a new raster file name
　Enter 'list' for a list of existing raster files
　Hit RETURN to cancel request
　> 553867.img

　　（この画像用の新しいファイル名を入力する）

　FLAG: Set the following flag?
　Verbose mode on.?(y/n) [n]

　　（フラグの設定，デフォルト(リターン)でよい）

　地図画像の不要な部分を取り除きたい場合は，先ほど示した図の内側（水色の部分）を示す台形と同じ形の領域を示すマスク用のラスターを作成し不要な部分を取り除きます（領域は地図により異なります）．

（２）地図画像を８つのレイヤーに分割

　地図画像は各ビットごとに意味を持っています．必要であれば，これらを８つのレイヤーに分ける作業をr.mapcalcを用いて行ないます．ｎbit目を取り出す基本式は，（画像　MOD　２n）／２n-1　です（MODは商のあまりを求める関数）．r.mapcalcの中ではMODは％に対応します．

　GRASS 5.0beta6 > r.mapcalc

　　（コマンドr.mapcalcの起動．プロンプトがmapcalc>に変わります．）

　mapcalc> 553867b8 = ("553867.img" % 256) / 128　↓

　　（８bit目を切り出す演算，553867.imgは画像ファイル名で，マスク処理をした場合は，そのファイル名を指定します．結果は553867b8に出力されます．）

　EXECUTING 553867b8 = ... 100%
　CREATING SUPPORT FILES FOR 553867b8
　minimum value 0, maximum value 1

　（８bit目が１の場合は結果が１，それ以外は０になります）

以下同様に

　mapcalc> 553867b7=("553867.img" % 128)/64　　（７bit目の計算）
　mapcalc> 553867b6=("553867.img" % 64)/32 　　（６bit目の計算）
　mapcalc> 553867b5=("553867.img" % 32)/16　 　（５bit目の計算）
　mapcalc> 553867b4=("553867.img" % 16)/ 8　　 （４bit目の計算）
　mapcalc> 553867b3=("553867.img" % 8)/ 4　 　 （３bit目の計算）
　mapcalc> 553867b2=("553867.img" % 4)/ 2　 　 （２bit目の計算）
　mapcalc> 553867b1=("553867.img" % 2)/ 1　 　 （１bit目の計算）
　mapcalc> 　　　　　　　　　　　(r.mapcalcを終了するにはリターン．）

以上で，地図画像の各レイヤーが完成します．

（３）画像のグルーピング

　それぞれの画像を１度に座標変換したりする場合などのために，コマンドi.groupを用いてグルーピングを行ないます．

　GRASS 5.0beta6 > i.group
　LOCATION: image25000 i.group MAPSET: map
　This program edits imagery groups. You may add data layers to, or remove
　data layers from an imagery group. You may also create new groups
　Please enter the group to be created/modified
　GROUP: 553867______________ (list will show available groups)

　　（新しく作成するグループ名（ここでは，553867とした）を入力します）

　AFTER COMPLETING ALL ANSWERS, HIT <ESC> TO CONTINUE
　(OR <Ctrl-C> TO EXIT)

　ESCキーを押します．

　553867 - does not exist, do you wish to create a new group? (y/n) [n] y

　　（はじめて作成するので．ｙを入力します）

　LOCATION: image25000 GROUP: 553867 MAPSET: map
　Please mark an 'x' by the files to be added in group [553867]

　　（画像の一覧が表示されるので，グループとして登録する画像の前にxを入れ，リターンを押します．）

　MAPSET: map
　　x_ 553867b1
　　x_ 553867b2
　　　　　：省略
　　x_ 553867b8
　　（すべて指定し終えたら，ESCキーを押してください．）

　Group [553867] references the following raster files
　-----------------------　　（指定した画像の一覧が表示されます．）
　553867b1 in map
　553867b2 in map
　　　　　：省略
　553867b8 in map
　-----------------------
　Look ok? (y/n) y　　　　　　　　　（問題がなければｙを入力します．）
　Group 553867 created　　　　　　　（これでグループが完成します．）

　　（つぎに，下のような画面が表示されるので，リターンを押し終了します．）

　LOCATION: image25000 GROUP: 553867 MAPSET: map
　 1. Select a different group
　 2. Edit group title
　 3. Include new raster files in the group
　 　 or remove raster files from the group
　 4. Assign colors to the group
　 5. Create a new subgroup within the group
　 RETURN exit

　>

（４）UTM座標系への投影

　ｘｙ座標系で作成した地図画像25000をUTMへ座標変換し移動させるためには，UTM座標系のデータベース・LOCATION等が必要です．これらはあらかじめ作成しておいてください．

　グルーピングされた画像を別の座標系に投影するには，i.targetおよびi.rectifyの２つのコマンドを用います．まず，UTM座標系をもつLOCATIONとMAPSET名（座標変換のターゲットと言う）をi.targetコマンドを用いて指定します．

　GRASS 5.0beta6 > i.target
　This program targets an imagery group to a GRASS database
　Enter group that needs a target
　Enter 'list' for a list of existing imagery groups
　Enter 'list -f' for a verbose listing
　Hit RETURN to cancel request
　> 553867　↓　　（i.groupで作成したグループ名を入力します．）
　Please select the target LOCATION and MAPSET for group <553867>
　CURRENT LOCATION: image25000 　 （現在のLOCATIONとMAPSETが表示されます．）
　CURRENT MAPSET: map

　
　TARGET LOCATION: topomap____ 　（目的とするターゲットのLOCATIONを入力します）
　TARGET MAPSET: map________　 （ターゲットのMAPSETを入力し，ESCキーを押します．）
　group [553867] targeted for location [topomap], mapset [map]

　つぎに，座標変換のためのパラメータを指定します．前述したように地図画像25000は，４隅の座標値が明確であるために，数値入力で座標変換時に使用する座標変換POINTSファイルを作成します．緯度経度→UTM座標値変換コマンドm.ll2uを用いて，この地図画像の４つの隅のUTM座標値を求めておきます．

　GRASS 5.0beta6 > m.ll2u bessel

　　（コマンドm.ll2uを起動し，楕円体のオプションbesselを入れておきます）

　Enter lon lat, one coordinate pair per line, in the format
　 ddd:mm:ss{E|W} dd:mm:ss{N|S}
　Enter the word <end> when done
　> 138:52:30E 37:10:00N　 　　　 （４隅の経度緯度座標値を入力します）
　311350.12100612 4115065.29092433 54　（UTM座標値の結果が表示されます）
　　　　　：　省略
　> end　　 　　　　　　　　　　　　 （コマンドの終了はendと入力する）

座標変換のためのPOINTSファイルは，グループディレクトリの中にある必要があります．そこで，現在のディレクトリを移動します．

　GRASS 5.0beta6 > cd $LOCATION/group
　GRASS 5.0beta6 > ls

　　（先ほど作成したグループ名（例では，553867）のディレクトリがあります．）

　GRASS 5.0beta6 > cd 553867

　　（グループ名の下に移動します．ここにはすでに，REFとTARGETが存在します）．

　viまたはxeditを用いて，上で求めた情報から下記のようなPOINTSファイルを作成します．

　GRASS 5.0beta6 > vi POINTS
　# image target status
　# east north east north (1=ok)
　#
　　73.000000 　75.000000 311350.121006 4115065.290924 1
　　75.000000 3774.000000 311557.573253 4124310.789833 1
　4511.000000 3774.000000 322643.673839 4124069.129211 1
　4513.000000 　75.000000 322448.439130 4114823.826961 1

（左側はGRASS画像のｘｙ座標値で，右側はそれぞれに対応したUTM座標値です．注意しなければならないのは，もとのtiff画像の座標（上で示した図中の赤い数字）は左上が原点ですが，GRASSのｘｙ座標では左下が原点となることです．このため，ｙ座標は下から上へと座標値を読み変えなければなりません．）

　最後に， i.rectifyあるいは，i.rectify2コマンドを用いて座標変換を実行します．これらの違いは，マニュアルを参照して下さい．なお，ここでは，i.rectifyコマンドの例を示します．

　GRASS 5.0beta6 > i.rectify
　Enter the group containing files to be rectified
　Enter 'list' for a list of existing imagery groups
　Enter 'list -f' for a verbose listing
　Hit RETURN to cancel request
　> 553867　　　　（変換する画像のグループ名を入力します）
　----------------------------------------------------------
　Please select the file(s) you wish to rectify by naming an output file

　　　　　（変換先の新しいファイル名(前と同じでも可)を入力します．）

　553867b1@map . . . . . . . . . . . . 553867b1______
　553867b2@map . . . . . . . . . . . . 553867b2______
　　　　　　　　　　　　：　省略
　553867b8@map . . . . . . . . . . . . 553867b8______
　(enter list by any name to get a list of existing raster files)

　　（入力後，ESCキーを押します．）

　Please select one of the following options
　1. Use the current region in the target location
　2. Determine the smallest region which covers the image
　> 1

　　（オプションの選択，ここでは１を選びます．）

　You will receive mail when i.rectify is complete

　　（場合により，かなり時間がかかるため，バックグランドで計算されます．sendmailデーモンが動いているマシンではメールで終了の返事が来ます．）

（５）画像ファイルの圧縮と再合成

　変換が正常に行なわれたことを確認するために，GRASSを終了し，変換先（ターゲット）のLOCATION，MAPSETで再度GRASSを起動します．変換が終了したかどうかはg.listコマンド等でリストを表示して確認して下さい．

　これらの画像は変換先の分解能によりかなり大きなファイルになることがあります．そこでr.compressコマンドを用いて圧縮します．ビットに分解した画像は圧縮が非常に有効です．

　GRASS 5.0beta6 > r.compress
　Enter the name of an existing raster file
　Enter 'list' for a list of existing raster files
　Enter 'list -f' for a list with titles
　Hit RETURN to cancel request
　> 553867b1　↓　　　　　　　（圧縮するファイル名を入力します）
　<553867b1>
　553867b1 is uncompressed. Would you like to compress it? (y/n) [y] y　↓

　　（圧縮の確認．ｙを入力します）

　COMPRESS [553867b1]
　DONE: compressed file is 68006199 bytes smaller　　（圧縮された結果を表示します．）

　各ビットに分解したレイヤーをもとに戻したい場合はr.mapcalcコマンドを用いて以下の計算を行ないます．

　GRASS 5.0beta6 > r.mapcalc
　mapcalc> topomap = 128*5538b8+64*5538b7+32*5538b6+16*5538b5
　　　　　　　　　　　+8*5538b4+4*5538b3+2*5538b2+1*5538b1

　　（計算結果のラスター名をtopomapとした例です．）

　EXECUTING topomap = ... 100%
　CREATING SUPPORT FILES FOR topomap
　minimum value 0, maximum value 255

最後に変換したラスター画像を表示し，確認します．

戻る