2.数値地図25000(地図画像)の入力

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

(1)xy座標系での入力

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

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

 grass4.2.1> 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]

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

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

(2)地図画像を8つのレイヤーに分割

 地図画像は各ビットごとに意味を持っています.必要であれば,これらを8つのレイヤーに分ける作業をr.mapcalcを用いて行ないます.nbit目を取り出す基本式は,(画像 MOD 2n)/2n-1 です(MODは商のあまりを求める関数).r.mapcalcの中ではMODは%に対応します.

 grass4.2.1> r.mapcalc 

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

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

  (8bit目を切り出す演算,553867.imgは画像ファイル名で,

   マスク処理をした場合は,そのファイル名を指定します.

   結果は553867b8に出力されます.)

 EXECUTING 553867.cut = ... 100%

 CREATING SUPPORT FILES FOR 553867b8

 minimum value 0, maximum value 1

 (8bit目が1の場合は結果が1,それ以外は0になります)

以下同様に

 mapcalc> 553867b7=("553867.img" % 128)/64  (7bit目の計算)

 mapcalc> 553867b6=("553867.img" % 64)/32  (6bit目の計算)

 mapcalc> 553867b5=("553867.img" % 32)/16  (5bit目の計算)

 mapcalc> 553867b4=("553867.img" % 16)/ 8  (4bit目の計算)

 mapcalc> 553867b3=("553867.img" % 8)/ 4  (3bit目の計算)

 mapcalc> 553867b2=("553867.img" % 4)/ 2  (2bit目の計算)

 mapcalc> 553867b1=("553867.img" % 2)/ 1  (1bit目の計算)

 mapcalc>            (r.mapcalcを終了するにはリターン.)

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

(3)画像のグルーピング

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

 grass4.2.1> 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

  (初めて作成するので.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        (問題がなければ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

 > 

(4)UTM座標系への投影

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

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

 grass4.2.1> 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を用いて,この地図画像の4つの隅のUTM座標値を求めておきます.

 grass4.2.1> 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 ↓    (4隅の経度緯度座標値を入力します)

 311350.12100612 4115065.29092433 54 (UTM座標値の結果が表示されます)

     : 省略

 > end                (コマンドの終了はendと入力する)

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

 grass4.2.1> cd $LOCATION/group 

 grass4.2.1> ls 

  (先ほど作成したグループ名(例では,553867)のディレクトリがあります.)

 grass4.2.1> cd 553867 

  (グループ名の下に移動します.ここにはすでに,REFとTARGETが存在します).

 viまたはxeditを用いて,下記のようなPOINTSファイルを作成します.

 grass4.2.1> 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画像のxy座標値で,右側はそれぞれに対応したUTM座標値です.注意しなければならないのは,もとのtiff画像の座標(上で示した図中の赤い数字)は左上が原点ですが,GRASSのxy座標では左下が原点となることです.このため,y座標は下から上へと座標値を読み変えなければなりません.)

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

 grass4.2.1> 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を選びます.)

 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デーモンが動いているマシンではメールで終了の返事が来ます.)

(5)画像ファイルの圧縮と再合成

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

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

 grass4.2.1> 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 

  (圧縮の確認.yを入力します)

 COMPRESS [553867b1]

 DONE: compressed file is 68006199 bytes smaller

  (圧縮された結果を表示します.)

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

 grass4.2.1> 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

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


戻る