GRASS

• GRASSのインストールのしかた
• 今日のGRASSコマンド
• RでGRASSを使う
• GRASSのアドオン追加

ポイントデータの読み込みとラスター化

土地被覆図精度検証などで, ポイントデータからなるラスターマップを作りたいことがよくある。例えば

105.46065139771166|15.746810408302036|1
105.5347232818659|15.749040852354911|1
105.56055831909887|15.739871092468674|1
105.58253097535079|15.731444461937985|1
... 

のように, 経度, 緯度, カテゴリIDという3つの項目がパイプ"|"で区切られたテキストファイル"vali.txt"がある。まずこれをベクターマップに読み込む:

v.in.ascii input=vali.txt output=vali x=1 y=2 cat=3 columns="x double precision, y double precision, cat integer" --o

そしてラスターに変換:

v.to.rast input=vali output=vali use=cat --overwrite

うまくできてるか確認:

r.stats vali -c
100%
1 10
2 10
3 10
...

このあと, d.rast valiとかで表示しても何も見えないかもしれないが, それはピクセルが細かすぎてつぶれて見えないだけかもしれない。不安なら, g.region res=で解像度を落としてやってみるとよろしい!

r.sun.mask

直達光が届くかどうか（自身の斜面方位だけでなく、他の斜面の影にならないか）を計算。

すげー遅い。かわりにr.sunを使うほうがよいらしい。http://osdir.com/ml/grass-user-gis/2012-10/msg00164.html

KML (Google Earthのフォーマット)で出力する

r.out.kmlというシェルスクリプトが公開されているのでそれをダウンロード。

他のマップセットのデータを使う

g.mapsets -p  いま自分が使っているマップセットのリストを表示。
motohka PERMANENT
g.mapsets -l  利用可能なマップセットのリストを表示。
kanto motohka tsukuba PERMANENT
g.mapsets addmapset=tsukuba   tsukubaというマップセットを追加。
g.mapsets -p 
motohka tsukuba PERMANENT

例えば、NDVI_2007001, NDVI_2007002,...,NDVI_2007365というNDVIデータ（365日分）から年最大値を抽出する処理

GRASSのデータディレクトリ（ロケーション名/マップセット名/cell）に移動して、
r.mapcalc out="-9999"
for i in NDVI_2007*; do
 r.mapcalc dummy="if(isnull($i), -9999, $i)"
 r.mapcalc out="if(out<dummy, dummy, out)"
done
r.mapcalc NDVI_max_2007="if(out==-9999, null(), out)"
g.remove dummy
（by T. Motohka）

現在のデフォルトのmapsetは何か?

mapsetがたくさんあると, 現在どのmapsetで仕事をしているかを知りたくなるよね
(新たに作ったマップはそこに入るので)!
g.gisenv MAPSET  （現在のデフォルトのmapset情報が表示される）
g.gisenv   （他の情報も表示される）
(by K. Nasahara)

GRASS起動時にでる邪魔なwindowを消す

GRASSをCUIモードで起動すればよい

方法1

$grass -text

方法2

$ vi /home/usrname/.grassrc6

とし,

GRASS_GUI : text   <--- と書き換える。
(by Kentaro Tanaka)

GRASS起動時のalias設定

$ cd /home/$username
$ vi .grass.bashrc     <----新規に作成し、alias等の設定を書き込む。
$ grass
(Kentaro Tanaka)

GRASS起動時のlsコマンドをカラー表示にする。

上記で作成した.grass.bashrcに以下を加える。
alias ls='ls --color'　
(Kentaro Tanaka)

ベクターデータにカラーバーを導入。

v.colors map=test3 column=dz rgb_column=rgb_dz rule=color_dz.txt

• column：色をつける際の指標となるカラム名
• rgb_column：何色かを保存するカラム名

↑以外はr.colorsとほぼ同じ。

表示方法は d.vect で -a オプションとrgb_columnを指定。

d.vect test3 rgb_column=rgb_dz -a

ベクターデータにアトリビュートを追加する。

データベースを取り込む。

db.in.ogr dsn=table.csv out=table1 key=ID  * key="name" …勝手に通し番号をふってくれる。

取り込む前に、同じな前で、拡張子 .csvt のファイル（table.csvt）を作ると取り込むときにデータの型を指定できる。

指定しないと、適当な型(私の場合は CHARACTER だった)で読み込んでくれる。

<table.csvt>
"Real","Real"      …Integer,Real,String,Date,Time,DateTimeが指定可能

<table.csv>
dz,h
0.666270,0.458210
0.427570,0.460540
0.000440,0.416240
0.000540,0.405200
0.865160,0.488950
1.554860,0.512910
2.106870,0.538120

ベクターデータにアトリビュートカラムを追加。

v.db.addtable test3 columns="dz double, h double"

データベースとベクターデータを結合。

v.db.connect test3 table=table1 key=ID -o

好きな形のマスクを作る

以下のコマンドを打ち込み、マウスでポチポチとエリアを指定する。

$ r.digit out=(出力ラスターマップ名) bgcmd="d.rast (バックに表示するマップ名)"

GMT用のgrdファイルを出力する

$ r.mapcalc test_nonull="if(isnull(test), 9999, test)"
$ r.out.bin -h input=test_nonull output=test.grd

教師なし分類する.....GRASS 6.2.3 (2010/04/15 松吉)

i.group,　i.cluster,　i.maxlik、3つのコマンドを使います。

まず、i.groupで分類に用いる指標 (NDVIとかGRVIとか) のラスターマップ (2つ以上) をグループ化。

こんな感じ。

$ i.group group=(グループ名) subgroup=(サブグループ名) input=(ラスターマップ名1,ラスターマップ名2,…)

※ラスターマップはグループの下に作られるサブグループに登録されます。

次に、i.clusterでクラスター分類法で教師なし分類し、各分類クラスの中心位置を求めます。クラスター分類法についてはgoogle先生に聞いてみましょう。

こんな感じ。

i.cluster group=(グループ名) subgroup=(サブグループ名) sigfile=(ファイル名) classes=(クラス数) 

• group、subgroupは、分類したいラスターマップを登録したグループ名とサブグループ名。
• sigfuileは、各分類クラスの中心位置を求めた結果を保存するファイル名を指定。/(LOCATION名)/(MAPSET名)/group/(グループ名)/subgroup/(サブグループ名)/sig/ に保存されます。
• classesは、分類するクラス数を指定。

• 各クラスの中心位置を求める過程でその初期位置を指定したいとき、指定したい位置が書かれたファイルを seed=(ファイル名) を追記して指定。ファイルはsigfuileと同じ書式で書かれたもの。ちなみに同じサブグループのsigフォルダに保存されたもののみ指定可能。
• 各クラスの中心位置を求める過程の履歴は、reportfile=(ファイル名)を追記するとでカレントディレクトリに保存できる。

最後に、i.maxlikで各クラスの中心位置をもとに、最尤法で各ピクセルを各グループに分類。

こんな感じ。

i.maxlik group=(グループ名) subgroup=(サブグループ名) sigfile=(ファイル名) class=(分類結果を出力するラスターマップ名)

• group、subgroupは、分類したいラスターマップを登録したグループ名とサブグループ名。
• sigfuileは、各分類クラスの中心位置を保存したファイル名。

画像中にテキストを表示する。

test.txtのテキスト情報を表示する。

d.text -b size=(フォントサイズ) color=(文字の色) at=(表示位置) < test.txt

-b : 太字表示

r.contour

DEM（ラスター）から等高線（ベクター）を作成。

GRASS7では, null値を色で塗るには, 以下のように-nとbgcolor=をつける必要あり：

d.rast -n DSM bgcolor=blue

d.mon start=PNG (2010/11/29 Kentaro Tanaka)

GRASSのデータをPNGとして吐き出す方法。d.out.png と違うのは、処理の際データを可視化する必要がない (CUI上の操作なのでshellで有効)。

使いかた

#!/bin/sh
#SCRNAME=test.sh

export GRASS_PNGFILE=/home/grass.png
export GRASS_TRUECOLOR=TRUE
export GRASS_WIDTH=900
export GRASS_HEIGHT=1200
export GRASS_PNG_COMPRESSION=1         
d.mon start=PNG
d.rast TEST_MAP
d.legend TEST_MAP
d.mon stop=PNG

を記述したシェルスクリプトの実行。

$ test.sh
$ ls  /home/
grass.png

Region情報の保存＆読み込み (2011/06/09 Kentaro Tanaka)

$ g.region save=test_area   <===現在のREGION情報がtest_areaに保存される。
$ g.region regi=test_area   <===test_areaからREGION情報を読み込む。

i.spectral

i.spectral は、GRASSコマンドの一つである。1つもしくは複数の座標について、複数の画像から画像値を取得し、その値を使ってグラフを書きたい時に使用する。出力はpng形式。

※事前にgnuplotのインストールが必要。

例) ALOS/AVNIR-2 の青、緑、赤、近赤外バンドの画像を使って、ある地点のスペクトル画像を書きたい場合。

GRASS > i.spectral -i  raster=band_1,band_2,band_3,band_4 out=test_image

実行結果

例) 森林

※うまく行かない場合: i.spectral を実行した際、下記のようなエラーが出る場合がある。(Keywordのみ抜粋)

set data style lines ~~~~~~~~~~~
~~~~~ Unrecognized option  ~~~~~

i.spectral のプログラムに書かれている内容がおかしい様ので、書き換える。

$ sudo vi /usr/lib/grass64/scripts/i.spectral

251行目を次の通り書き換える。

変更前
echo "set data style lines" >> spectrum.gnuplot

変更後
echo "set style data lines" >> spectrum.gnuplot

エラー

error setting region: child process exited abnormally
GRASS起動時にこんなエラーが出たら、regionファイルが変になったっぽいので、
g.region -d
これで解決！

GRASSコマンドを含むシェルスクリプトを, GRASSの外から実行する:

• まず, スクリプトを書き, chmod +xする。
• 次に, そのスクリプトの絶対パスを以下のように環境変数GRASS_BATCH_JOBに登録:

export GRASS_BATCH_JOB=/home/nishida/hogehoge.sh

• そしてコマンドラインからgrassを呼ぶ。その際, 以下のようにlocation, mapsetを指定する:

grass /home/nishida/latlon/PERMANENT

• そうすると, 指定されたlocation/mapsetでスクリプトが実行され, 終了するともとのシェルに戻る。