2007-11-06:GRASSからRを使ってみましょう - BioEoS::環境解析学研究室::

01.イントロ

まずは三年生へ宿題です．
先日行っていただいた仕事の続きと勉強を兼ねてハザードマップについて調べてください．
具体的には，「災害学習情報などにはどういったものがあるのか」「今あるハザードマップに足りない情報は何か」などを調べてもらいたいと思います．

本日は，前回の「R」の続きとして，GRASSの情報をRから見てみたいと思います．
※今回は，最後のほうで「R」を利用するため，日本語環境を表示できる「mlterm」で作業を行うようにしてください．

02.過去に作成したlocation「Isebay」の削除とディレクトリの移動

フォルダを削除するコマンドを覚えていますか？
「rmdir」でしたね．
しかしこのコマンドの場合，フォルダの中にファイルがあると削除することが出来ません．
そこで，「rm -r」を利用します．
※このコマンドは，強制的にフォルダ内のすべてのファイルを削除します．このコマンドを入力するときは十分注意してください．

まず，ホームディレクトリにいることを確認してください．
そこで，

rm -r ~/GRASS/Isebay

と入力します．（「~」はホームディレクトリという意味です．）

ls GRASS

で「Isebay」が削除されていることを確認してください．

過去に利用した「dem.bin」があるディレクトリを検索します．

locate dem.bin

すると，「dem.bin」がある場所の一覧が表示されます．
その中の一つに移動してください．
今回はうすいさんのファイルをお借りしましょう☆

cd /home/mahito/GRASS/

「dem.bin」のファイルの大きさを確認しましょう．

ls -al

ファイルサイズが「42240000」であればそのファイルで正解です．

03.dem.binの取り込みと海岸線の作成

まずGRASSを立ち上げます．

GRASS62

「Projection values」をクリックします．
今回からは，座標系についても意識していきましょう．

①LOCATION:Isebay
②MAPSET:shibata
③Would you like create location?:「y」
④Do you have all this information?:「y」
⑤please specify the coordinate system for location:「b」→「y」
⑥Please enter a one line description for lacation:適当に(今回は「Isebay」)→「y」
⑦Do you with to specify a geodetic datum for this location?:「y」
⑧「list」→「スペース」で移動．[tokyo]があることを確認．
⑨「tokyo」と入力．
⑩「list」→「1. Used in Japan」があることを確認．
⑪「1」と入力．
⑫North : 35:35N
  South : 33:35N
  West : 135:15E
  East : 138E
  East-West : 0:0:02.25
  North-South : 0:0:01.5
 ※total rows=4800，total cols=4400となればOKです．

さっそく「dem.bin」を取り込みます．
まず，モニターを立ち上げておきます．

d.mon x0

「dem.bin」を取り込みます．

r.in.bin -b in=dem.bin out=dem n=35:35N s=33:35N e=138E w=135:15E r=4800 c=4400 by=2

「dem」が作成されているか確認します．

g.list rast

今作成した「dem」を表示させます．

d.rast dem

海岸線（coast）を作成するために「reclass」を行い，海域を「0」，陸域を「1」という画像を作ります．

r.reclass in=dem out=land_sea

• 999 thru -3 = 0
• 2 thru 2187 = 1

end

表示させて見ましょう．

d.rast land_sea

海岸線を作成します．

r.contour in=land_sea out=coast lev=1

表示させましょう．

d.vect coast

04.画像の取り込み

では，画像の取り込みを行います．
GISで扱う情報には必ず「位置情報」が含まれます．しかし，画像には一般的に位置情報は含まれません．これでは，GISに取り込むことはできません．
つまり，GISで画像を扱うためには単なる画像ファイルに位置情報を持たせる必要があります．
画像を扱う場合，一番圧縮率の低い「TIFF」画像を利用します．さらに，「位置情報」を持たした「TIFF」画像のことを「GeoTiff」といいます．
また，「jpg」から「tiff」に変換するなど，画像の圧縮方法を変更することを「convert」と言います．

単なる画像に位置座標を持たせるために変換をかけなくてはいけません．
この変換のことを「affin変換」と言います．さらに，「affin変換」の結果の係数などを並べたファイルのことを「ワールドファイル（.tfw）」と言います．
この「ワールドファイル」と「Tiff」画像を合わせて初めて「GeoTiff」画像となります．そのコマンドは例えば，

geotifcp -e 5035.tfw 5035.tif g5035.tif
※今回は行いません！！

流れを簡単に示すと以下のようになります．

jpg画像→「convert」→tiff画像→「affin変換」→ワールドファイル→「ワールドファイル」＋「tiff画像」→GeoTiff画像

今回は，先生がすでにGioTiff画像を作成してくださっているのでそれを利用させていただきます．
まず，ディレクトリを移動します．（GRASSから出る必要はありません．）

cd /home/kaoru/GIS/Geology/TIF

「list」をとり，「g○○○○.tif」があることを確認します．

ls

「g○○○○.tif」をそれぞれGRASSの中に取り込みます．

r.in.gdal -o in=g5035.tif out=g5035
※「r.in.gdal complete.」と表示されればOKです．
※順に「g5135」「g5235」「g5036」「g5136」「g5236」も行ってください．
※gdal:Geospatial Data Abstraction Library

試しに表示させて見ましょう．

d.rast g5035
※左下に小さく表示されればOKです．

これらの6つの画像を一つの画像に合成します．

r.patch in=g5035,g5135,g5235,g5036,g5136,g5236 out=geology

海岸線と一緒に表示させて見ましょう．

d.rast geology
d.vect coast

05.RとGRASSの連携

自分のホームディレクトリに移動します．（GRASSから出る必要はありません．）

cd

先週作ったディレクトリ「R_test」に移動します．

cd R_test

ディレクトリ「Isebay」を作りそこに移動します．

mkdir　Isebay
cd Isebay

先ほど「geology」を描画した際，画像を含んでいない地域も表示されていました．これでは，Rなどで解析する際，無駄な領域まで解析することになり，無駄な労力がかかります．
そこで，解析範囲を狭めてみましょう．
まず，現在の領域や分解能など（リージョン）を表示させて見ましょう．

g.region -p

変更します．

g.region n=35N s=34N e=137E w=136E ns=0:0:12 ew=0:0:18

変更できているか確認してみましょう．

g.region -p

このリージョンを保存します．

g.region save=shima

モニターの画像を消します．

d.erase

上向き矢印で「d.rast geology」と「d.vect coast」を探し決定します．

d.rast geology
d.vect coast

表示されている範囲が先ほどより小さくなっているのがわかりますか？
このように，リージョンを変更した場合，この領域でしか画像などの取り込みや解析は行われないため注意が必要です．
リージョンを最初の設定を行ったときの値に戻す方法を練習しておきましょう．

g.region -d

確認してみてください．

g.region -p

今までに作った設定ファイルを確認してみましょう．

g.list region

先ほど作ったリージョン「shima」に戻します．

g.region regi=shima

では，「R」を立ち上げてみましょう．（GRASSからは出ません．）

R

まず，GRASSと連携させるためのおまじないをいくつか入力します．

library(spgrass6)
G <- gmeta6()

grassで作った「geology」と「dem」をRで利用します．

ise <- rast.get6(c("geology","dem"), cat = c(TRUE,FALSE), ignore.stderr = TRUE)

これは，「ise」という変数を作って，そこにラスターデータである「geology」と「dem」を入れるという意味です．
後半の，cat=c(TRUE,FALSE)はカタログ番号の有無を指定しています．（有：TRUE，無：FALSEです．）カタログ番号とは，例えば以下のようなものです．

●今回は地質のデータなので
   カテゴリー名　　　カタログ番号
　　　花崗岩　　　　　　　１
　　　玄武岩　　　　　　　２
　　関東ローム　　　　　　３
といったイメージで，それぞれのカテゴリーに対して固定して割り当てられている番号のことです．
※現段階では，まだカタログ番号は割り振られていません．

「dem」の場合単なる標高値であるためカテゴリー番号はありません．

本格的な解析は来週にするとして…
今回は表示までを行いましょう．

image(ise,attr=2,col=terrain.colors(20))
title("Ise_Bay")

また「R」ではGRASSの情報も見ることが出来ます．

str(G)

では，最後に全てを終了させましょう．

q()→「y」でワークスペースを保存してください．

保存しておくことで，来週から同じ作業を行わなくて済みます．
GRASSからでます．

exit

以上で今日のゼミは終わりです．
お疲れ様でした☆