GSMaPのHDFファイルからGeotiffの作成

GSMaPのHDFファイル(h5)からRにてGeotiffを作成しました。
HDFファイルはGRASSや他のリモセンソフトでGISに取り込むことができるようですが、GSMaPで試してもうまく行かなかったので、Rのrhdfパッケージでデータを読み込み、rasterパッケージでGeotiffを作成しました。
スクリプトは全球合成降水マップレベル3（月平均)の任意の地点の降水量(mm/month)をGeotiffに変換できるように作成したので、他の種類の画像の場合、改編する必要があります。
####必要なライブラリの読み出し####
#source("http://bioconductor.org/biocLite.R")
#biocLite("rhdf5")
library(rhdf5)
library(raster)
####設定値####
#データの場所
in.data.dir         <- "datafiledir/"
in.data.file.common <- "GPMMRG_MAP_"
in.data.file.para   <- c("1403_M_L3S_MCM_03A", "1404_M_L3S_MCM_03A", "1405_M_L3S_MCM_03A", "1406_M_L3S_MCM_03A", "1407_M_L3S_MCM_03A", "1408_M_L3S_MCM_03B", "1409_M_L3S_MCM_03B")
in.data.file.ext    <- ".h5"
in.data.name.para   <- c("1403_M_L3S_MCM_03A", "1404_M_L3S_MCM_03A", "1405_M_L3S_MCM_03A", "1406_M_L3S_MCM_03A", "1407_M_L3S_MCM_03A", "1408_M_L3S_MCM_03B", "1409_M_L3S_MCM_03B")  # in.data.file.paraと同じ長さのベクトルであること。
out.data.dir <- "datafiledir/"
out.data.file.common <- "GPMMRG_MAP_"
out.data.file.para   <- c("1403_M_L3S_MCM_03A", "1404_M_L3S_MCM_03A", "1405_M_L3S_MCM_03A", "1406_M_L3S_MCM_03A", "1407_M_L3S_MCM_03A", "1408_M_L3S_MCM_03B", "1409_M_L3S_MCM_03B") # in.data.file.paraと同じ長さのベクトルであること。
out.data.file.ext    <- ".tif"
#月別降雨量にするため、各月の日数を記載
data.unit.conversion.days <- c(31, 30, 31, 30, 31, 31, 30)
data.unit.conversion.hour <- 24
#HDFファイルの内、利用するデータ
in.data.dataset <- "Grid/monthlyPrecipRateGC"
#データの開始緯度経度
origin.lon <- -180
origin.lat <- -90
in.data.size.lon <- 3600
in.data.size.lat <- 1800
#解像度
resolution.lon <- 0.1
resolution.lat <- 0.1
#ほしい地域の緯度経度
target.lonlat.sw <- c(-18, 12)
target.lonlat.ne <- c(-11, 17)
####################################
for(i in 1:length(in.data.file.para)){
#データの読み出し
    filename <- paste(in.data.dir, in.data.file.common, in.data.file.para[i], in.data.file.ext, sep='')
    in.data.h5 <- h5read(paste(in.data.dir, in.data.file.common, in.data.file.para[i], in.data.file.ext, sep=""), in.data.dataset)
#データの内、必要な範囲を抽出
#hdfファイルをhdfviewで見ると、１行に南から北へデータが格納されているが、Rで読み取るとマトリクスの列に南から北へ格納される。
#rastarへの書き出し時には、マトリクスの行を北から南に記載する必要が有るため、行番号は反転させる。
    #データが入っている行
        target.data.row <- c(((target.lonlat.ne[2] - origin.lat) / resolution.lat):((target.lonlat.sw[2] - origin.lat) / resolution.lat + 1))
    #データが入っている列
        target.data.col <- c(((target.lonlat.sw[1] - origin.lon) / resolution.lon + 1):((target.lonlat.ne[1] - origin.lon) / resolution.lon))
     #必要な範囲のマトリクス作成
        target.data <- in.data.h5[target.data.row, target.data.col]
*データの特定の範囲を抽出するには、h5readにindex=list(2:3,2:3)オプションを付けることで部分的に抽出できるようです。
#ターゲットエリアのテーブルをgeotiffに書き出し
    #データ格納用空のラスタ作成  
        out.nrow <- (target.lonlat.ne[2] - target.lonlat.sw[2]) / resolution.lat
        out.ncol <- (target.lonlat.ne[1] - target.lonlat.sw[1]) / resolution.lon
        (x <- raster(ncol=out.ncol, nrow=out.nrow, xmn=target.lonlat.sw[1], xmx=target.lonlat.ne[1], ymn=target.lonlat.sw[2], ymx=target.lonlat.ne[2]))
    #データをラスタ内に格納
        values(x) <- target.data * data.unit.conversion.days[i] * data.unit.conversion.hour
    #Geotiffとして書き出し
        writeRaster(x, filename=paste(out.data.dir, out.data.file.common, out.data.file.para[i], out.data.file.ext, sep=""),format="GTiff", overwrite=TRUE)
}

RでGeotifの面積を出す

#ラスタデータの面積計算
#GISでできると思ったのにうまくいかないため、
#Rでセルの数を数えて面積に変換する。
#傾斜(度)のラスタから、任意の区分で区切った面積を%で出す。
#GeotiffのプロジェクションやセルサイズのデータをRから呼び出す方法が不明。

####parameter####
in.file.path         <- "C:/path/hoge.tif"

pixel.size.x        <- 30.31054
pixel.size.y        <- 30.30613

#特定の値のピクセル面積を出す場合
pixel.value.target  <- c(1:14)

#特定の範囲の値のピクセル面積を出す場合
pixel.value.target.min <- c(0, 5.710593137, 14.03624347, 16.69924423, 30.96375653)  #min, maxのi番目が、目的の階級に対応するように。
pixel.value.target.max <- c(5.710593137, 14.03624347, 16.69924423, 30.96375653, 90)
pixel.value.target.lab <- c("0-10","10-25%", "25-30%", "30-60%", "OVER60" )
#################

library(raster)
#ファイルを読み込み
i <- 1
target.tif <- raster(in.file.path)

#パラメータから変数生成
pixel.area <- pixel.size.x * pixel.size.y

#ベクトルデータに変換 <-これ必要？<-計算が早く終わる気がする
target.v <- as.vector(target.tif)
#NAのないベクトルデータに変換 <-これ必要？<-計算が早く終わる気がする
target.v.nna <- target.v[!is.na(target.v)]

#任意の値のセルの面積
    #各値をもつ面積をデータフレーム形式で出力する
    i <- 1
    tmp.area <- 0
    tmp.name <- 0
    for(i in 1:length(pixel.value.target)){
        tmp.area[i] <- sum(target.v.nna == pixel.value.target[i]) * pixel.area
        tmp.name[i] <- paste("pv", pixel.value.target[i], sep="")
    }
    (area.spec.range <- data.frame(area = tmp.area, row.names=pixel.value.target.lab))

#条件に合うピクセルの面積
    #条件に当てはまる値をもつピクセルの合計面積をデータフレーム形式で出力する
    i <- 3
    tmp.area <- 0
    tmp.name <- 0
    for(i in 1:length(pixel.value.target.min)){
        tmp.area[i] <- sum(target.v.nna > pixel.value.target.min[i] & target.v.nna <= pixel.value.target.max[i]) * pixel.area
        tmp.name[i] <- pixel.value.target.lab[i]
        #tmp.name[i] <- paste(pixel.value.target.min[i], "-", pixel.value.target.max[i], sep="")
    }
    (area.spec.value <- data.frame(area = tmp.area, row.names=tmp.name))

Editing geitif file with R script

Rでgeotifファイルの値を書き換え
インストールしているGRASSの調子が悪いため、Rを使ってgeotiffフィアルの値を変更した。
今回は、NA値担っているセルをすべて0に変更した。
####parameter####
in.file.dir          <- "C:/Users/XXXX/indir"
in.file.name.common  <- "example_"
in.file.name.para    <- c("001", "002", "003")
in.file.name.ext     <- ".tif"
out.file.dir         <- "C:/Users/XXXX/outdir"
out.file.name.common <- "example_NNA"
out.file.name.ext    <- ".tif"
#################

library(raster)
i <- in.file.name.para[1]

for(i in in.file.name.para){
    target.tif <- raster(paste(in.file.dir, "/", in.file.name.common, i, in.file.name.ext,sep=""))
    target.tif[is.na(target.tif)] <- 0
    writeRaster(target.tif, filename=paste(out.file.dir, "/", out.file.name.common, i, out.file.name.ext, sep=""),format="GTiff")
}

Frequency Tabulations from geotif (with R)

GISで、ポテンシャル地域を評価した後、ポテンシャルインデックスの累積累積度数分布表を作成するために、Rにてtiffファイルを操作した。
#ライブラリの読み込み
library(raster)

#GISで出力したtiffファイルの読み込み
pot.tif <- raster('XXXXXX.tif')

#画面で確認
plot(pot.tif)

#方はSで読み込まれているので、数値に変換する。(必要か？)
#また、その際データなしの場所を飛ばす。（必要？）
mode(pot.tif)
pot.num <- as.vector(pot.tif[!is.na(pot.tif)])
#モードの確認。
mode(pot.num)
#ヒストグラム作成
hist(pot.num)

#数値データをカテゴリに分ける
pot.cat <- cut(pot.num, breaks=c(14:30)/2, labels=(c(15:30)/2))
(pot.table <- table(pot.cat))

#エクスポート(略)

QGISでの表示言語の変更(QGIS 2.4)

海外の人にQGISの使い方を示すとき、メニューが英語表記であると説明しやすいため、UIの言語設定を変更しました。
参考　http://deerfoot.exblog.jp/12372794/
設定(s)→オプション→ロケールから、「システムロケールを上書きする」にチェックをいれて、変更したい言語を選択。
OKを押して、QGISを再起動すると、表示言語が変更されます。

Visibility on SAGA GIS

SAGA GISの機能で、特定の地点からどこの場所が見ることができるか計算するモジュールがあります。このモジュールを使えば、たとえば特定の山からどの範囲が展望できるかわかります。利用するモジュールは、Terrain Analysis – Lighting, Visibility / Visibility (Single Point) [Interactive]です。

いつも通り、DEMのグリッドを指定します。

実行すると、下のメッセージに、”Interactive module execution has been started” と表示されます。Interactiveのモジュールは、計算起点等を指定してやる必要があります。地図を表示し、地図の上の矢印(Action)ツールで、地図の適当な地点を選択すると、計算が始まります。

 

一回のクリックで一か所からのVisibilityが計算され、再び待機状態になるので、Moduleメニューから、Visivility モジュールを停止させます。

Visivilityグリッドが作成されるので、地図に重ねてみます。

白くなっているところが、選択した箇所から見えるところです。3D表示すると以下の通り。

できました。

ImageMagickでGIFアニメーション

SAGA GISでの3D表示機能には、任意の視点をつなげた鳥瞰図をアニメーション表示する機能があります。SAGA GISではそれらを連続したスナップショットとしてpng出力ができるので、これをImageMagickでGIFアニメとして結合します。

まずは、出力されたpng画像をそのままつなげるとサイズが大きくなるのでリサイズします。

convert out*.png -resize 50% resize_%06d.png

新しくサイズが半分になったpngが作成されます。これを、以下のコマンドでGIFアニメーションに変換します。

convert -delay 5 -loop 0 resize_*.png movie.gif

できました。

Watershed Basin in SAGA GIS

Sink Removal SAGA GISで、集水域図を作成しました。Saga Users Guide vol.2を参考に。

使うツールは以下の通り

• Terrain Analysis - Preprocessing/Sink Removal
• Terrain Analysis – Hydrology/Catchment Area(Parallel)
• Terrain analysis – Channels/Channel Network
• Terrain analysis – Channels/Watershed Basins

作業の流れとしては、１）DEMから窪地を削除、２）各セルにおける集水面積の計算、３）河川流路の作成、４）集水域グリッドの作成といった感じです。

データは、ASTER GDEMを利用しましたが、面積の計算等があるため、プロジェクションはUTMに変換して置きます。

Sink Removal

ModuleのTerrain Analysis - Preprocessing内にあるSink Removalツールにより、標高グリッドから窪地を削除します。ツールを選択すると、Decriptionに概要が表示される。このツールはグリッドを指定すると、窪地を埋めたグリッドを出力します。また、Thresholdを指定すると、どれくらいの標高差まで乗り越えを許容するか指定できます。

標高グリッドを指定して実行します。実行すると、河川が流下できるように、窪地や山岳地域で沢がうまく流下できない箇所が、埋められます。下の図は、左が処理前、右が処理後です。

湖周辺部や、火口も以下のように埋められます。

 

 

Catchment Area(Parallel)

上で作ったグリッドに対し、各セルへの集水面積が入ったグリッドが作成されます。利用するのは、Terrain Analysis – Hydrology/Catchment Area(Parallel) です。標高データとしては、先ほど作成したno sinkのグリッドを指定します。

OptionのMethodでは、河川の流れがDEMの中でどう流れる仮定をするか設定できます。 河川の分岐を許可するかどうかや、流下方向を8方位に限定するか、などが指定できます。どの仮定を用いるかは、やってみて不自然でないものを選べばよいようです。

実行結果すると、それぞれのセルにそのセル自身の面積とそのセルの常流側に位置するセルの合計の面積が入力されます。色づけを調整(Logarithmic(up) 10000)して、Hillshadingレイヤを重ねて表示してみると、以下のようになります。

さらにズームして、セルに入っている値を表示すると、各セルの値が確認できます。

 

Channel Network

Sink RemovalとCatchment Area(Parallel)を用いて作成したグリッドにから、河川ネットワークのグリッドとシェープを作製します。モジュールは、Terrain Analysis – Channels / Channel Networkのモジュールを利用します。

Initiation Gridには、河川が始まる起点を指定します。今回は、集水域が1km以上になったら、河川が発生すると仮定します。Initiation GridにCatchment Area(Parallel)で作成したグリッドを指定し、Greater than 1000,000とします。また、最短河川長をセル数として与えます。今回は10としました。

実行すると、河川ネットワークを表したグリッドとシェープ、各セルの流下方向を示したグリッド、が作成されます。

グリッドの各セルには、河川の識別番号が与えられ、各河川の合流点には-1が与えられます。下の図は河川ネットワークのグリッドとシェープを重ねて表示しました。河川9と河川1が各端点-1のグリッドで合流して河川10になっています。

Watershed Basins

Terrain Analysis – Channels / Water Shed Basinsモジュールを使って、これまで作成したグリッドから集水域を作成します。標高グリッドと、先ほど作成したChannel Networkを指定します。

OptionのMin. Sizeは、最少集水面積をセルの数で指定します。今回は1000セルとしました。

できました。

SAGA GISでのグリッドの計算

複数のグリッドを用いて、計算を行います。
今回は一番簡単な二つのグリッドの差分を取る計算をします。Water Basinを作成する際に作成したfill sinkを行ったDEMと、オリジナルのDEMを比較しました。
モジュールは、Grid – Calculation / Grid Defference　を利用します。
窪地を埋められたDEMをA、オリジナルのDEMをBとして指定します。

差分のグリッドが出力されました。色をGraduated Colorでwhite-Redに指定し、陰影図と重ねてみました。

オリジナルと比べると、どこがどれくらい埋められているかよくわかります。

SAGA GIS 3D View 2

前回、3D Viewが失敗してしまいましたが、座標系が緯度経度であったため、できなかったみたいです。UTM座標へ変換したら、表示されました。

3D表示をやめるには、地図を右クリックし、Show 3D Viewのチェックを外します。

デフォルトの設定では荒いので、3Dの設定画面のResolusionを1000くらいにしてみます。


きれいになりました。
また地形を強調するため、縦に引き伸ばします。

この右から3つ目のアイコンで、強調されます。

そのほか、赤青メガネで立体視できる画像を作ったり、

ムービーもつくれます。

SAGA GISでCRSの変更

TASTER GDEMはWGS84/EGM96で作成されている緯度経度座標のデータのため、これをUTMに変換します。
DEMをSAGA GISに読み込みます。
変換に用いるのは、Projection – Proj.4 のProj.4 (Dialog, Grid)モジュールです。
モジュールを起動すると、変換元(Source)と変換後(Target)のCRSを指定するダイヤログが出てきます。

Source Projection Parameters の28 Parameters …をクリックすると、変換元のCRSの詳細を入力するダイヤログが現れます。
変換元のProjectionをlat/long(Geodetic)、DatamをWGS84と指定します。

OKを押すと、ひとつ前のダイヤログにもとるので、次にTarget Projection Parameterを指定します。今回はUTMに変換するので、UTMを指定します。

OKを押すと、再度元のダイヤログに戻ります。

プロジェクションの設定ができたため、OKを押すと、変換後のUTMゾーンを聞かれるため、該当するZONE 54を指定します。

OKを押します。
グリッドを描写する空間の設定が聞かれます。

OKを押すと、変換が開始されます。

ワークスペースに、新しいグリッドが作成されました。

SAGA GIS 3D View

前回作った陰影図をSAGA GISで3D表示してみます。
3D表示をするのは、3Dと書かれたアイコンをクリックするだけ。

DEMを指定するダイヤログが出てくるので、指定します。ほかのパラメータはデフォルトのままで大丈夫です。

すると、、、

あれ、失敗です。

Ubuntuでキーの入れ替え

Ubuntu 14.04をインストールしたので、まずはキーボードのカスタマイズを行いました。
ネットで検索すると、GUIによる設定もできると書かれていましたが、14.04にはそのメニューがなかったので、設定ファイルによる変更を行いました。

参考
http://d.hatena.ne.jp/adragoona/20130518/1368850862
http://bearmini.hatenablog.com/entry/2013/07/12/161637


まず、物理的なCapslockキーのキーコードを以下のコマンドで確認します。
$ xev

イベントテスターなるウィンドウにフォーカスした状態で、任意のキーをタイプすると、キーコードを含む出力がされます。
例えば、以下の通り。

KeyPress event, serial 38, synthetic NO, window 0x1600001,
    root 0xa0, subw 0x0, time 4787925, (144,65), root:(1494,296),
    state 0x0, keycode 66 (keysym 0xff30, Eisu_toggle), same_screen YES,
    XLookupString gives 0 bytes:
    XmbLookupString gives 0 bytes:
    XFilterEvent returns: False

Capslockキーはkeycode 66であることがわかります。


以下のコマンドで現在のキーバインドを確認します。d

$ xmodmap -pke

すると

keycode  24 = q Q q Q

と行ったものが出力されます。
４つ表示されるのは、修飾キー無し、Shiftキーによる修飾、Mode_switchによる修飾、Mode_switch+Shiftキーによる修飾を意味しています。Mode_switchは、以下で確認できます。

$ xmodmap -pke | grep -i mode_switch
keycode 203 = Mode_switch NoSymbol Mode_switch

キーボードによってはキーコード203入力用のキー（AltGr）があるものもあるようです。

モディファイヤキーの確認
Ctrlキーなどの修飾キーは、モディファイヤキーとして特別に登録されている。モディファイヤキーは、以下で確認できます。

$ xmodmap
xmodmap:  up to 4 keys per modifier, (keycodes in parentheses):

shift       Shift_L (0x32),  Shift_R (0x3e)
lock        Eisu_toggle (0x42)
control     Control_L (0x25),  Control_R (0x69)
mod1        Alt_L (0x40),  Alt_R (0x6c),  Meta_L (0xcd)
mod2        Num_Lock (0x4d)
mod3     
mod4        Super_L (0x85),  Super_R (0x86),  Super_L (0xce),  Hyper_L (0xcf)
mod5        ISO_Level3_Shift (0x5c),  Mode_switch (0xcb)

Ctrlキーがモディファイヤキーに指定されているため、このキーを入れ替えるためには一度モディファイヤキーの指定を解除する必要があります。それには、以下を.Xmodmapの先頭に追加します。

以上をまとめると、以下を~/.Xmodmapに記述します。

! Ctrlキー、CapsLockをモディファイヤキーから解除
remove control = Control_L
remove lock = Eisu_toggle

! 左Ctrlであるキー(37)をCapslockに、Capslockキー(66)を左Ctrlに割り当て
keycode 37 = Eisu_toggle Caps_Lock Eisu_toggle Caps_Lock
keycode 66 = Control_L NoSymbol Control_L

! 再度Ctrlキー、CapsLockをモディファイヤキーへ登録
add control = Control_L
add lock = Eisu_toggle

$ xmodmap .Xmodmap
で反映します。

Contour Map on SAGA GIS

SAGA GIS等高線図を作成しました。いつものASTER GDEM¹を使います。

SAGA GISでのコンターマップの作成

等高線を引くためのモジュールは、Shapes – Grid / Contour Lines from Gridです。

設定項目は、

• 標高データのグリッド
• 最少標高
• 最高標高

を指定します。最高標高と最低標高の指定は、特定の高さの範囲の等高線を引きたいときや、海底地形を含むグリッドデータで、海底のみのコンターを引きたいときとかに便利です。今回はASTER GDEMの画像一枚すべてを作成してみました。
できました。

ズームしてみるとこんな感じです。

国土地理院のウォッちずで、正解の地図を確認しようとしたところ、

平素より、地図閲覧サービス（ウォッちず）をご利用頂きありがとうございます。
ウォッちずは平成26年4月1日以降地図を更新しておらず、平成26年10月1日をめどにサイトを停止する予定です。
なお、最新の地図は「地理院地図」でご覧いただけます。

とのことです。
地理院地図は、ウォッちずで提供されていた1:25,000の地図のほか、基準点や土壌図、断層、火山情報等のさまざまな情報をオーバーレイできるように整備したサイトのようです。しかも、webサイトへの埋め込み機能があるようです。
上の地図は北岳周辺なので、地理院地図と比較してみました。
1:25,000の地図はさすがに地表の細かい崖まで書かれていますが、20万分の1程度の地図であれば、そこそこの再現率ではないでしょうか。

---

1 ASTER GDEMはMETI/NASAに帰属します。^↩

SAGA GISで陰影図の作成

DEMから陰影図を作成します。DEMはASTER GDEM^*1を利用しました。日本で一番ダウンロードされているだろうと思われる、富士山を含む地域のDEMです。

DEMの読み込み

まずはDEMをSAGA GISに読み込みます。

workspaceのDataタブにDEMデータをドラッグアンドドロップします。ASTEER GDEMをダウンロードした際に同梱されるnumファイルは、標高データを作成するために何枚の衛星画像を用いたかや、欠損値や異常値の補完に用いた画像に関するデータが書かれているため、今回は利用しません。
データツリーに追加されたDEMファイルをダブルクリックすると、ワークスペースのMapsタブにDEMが追加され、地図が表示されます。

陰影図の作成

読み込んだDEMにTerrarian Analysis – Lightning, Visibility 内のAnalytical Hillshadingモジュールを利用します。

入力としてDEMを指定す他はデフォルトのパラメータで、実行します。(上記画像参照)。
ワークスペースのデータタブにAnalytical Hillshadingのグリッドデータが追加されるため、ダブルクリックしてDEMと同じ地図に追加します。

色づけ

色づけしたいレイヤをMapsタブで選択します。まずはDEMから作業をするので、Analytical Hillshadingを、ダブルクリックして非表示にします。非表示レイヤはmapsタブ内では[ ]でくくられます。
DEMのレイヤを選択すると、右側のウィンドウ(Object Propery Window)のsettingタブでDEMの色を設定できるようになります。colorsをクリックして出てくる色設定画面から、Presetを選択、defoltを指定します。

Objet Property Windowに戻り、Applyを押すと、DEMに色が付きます。

Analytical Hillshadingを再度ダブルクリックして表示させ、colorsのTypeでShadeを選択し、Applyします。

色合いはあんまりよくないけれど、完成です。
*1:ASTER GDEMはMETI/NASAに帰属する

R コマンドメモ

Rで使った、今後もちょくちょく使いそうなコマンドをメモ。自分で使うときに見返すためのメモのため、使い方の中のコマンドには、事前にxやiに適当な変数が代入されていたり、for loop内での使用や直前にグラフの描写がされていることをを想定していたりします。

カテゴリ	説明	コマンド	使い方
作業準備	ワークスペースのオブジェクトをリストに出力	ls()	ls()
	ワークスペース内のオブジェクトを削除	rm(list=ls())	rm(list=ls())
	グラフを消す	graphics.off()	graphics.off()
	作業ディレクトリの表示	getwd()	getwd()
	作業ディレクトリの変更	setwd()	getwd()
データの読み込み	テキストファイルのデータを読み込む	read.table()	read.table("data\\rain.tsv", row.names=1, header=TRUE)
テキスト操作	変数に入っている文字列、変数等をつなぎ合わせる	paste()	paste("Type", x, sep= "-")
	各列に同じ操作を繰りかえし行う.	apply()	apply(x, 2, sd) または apply(x, 2, function(i){         runs.test(as.factor(i < mean(i)))     } )
オブジェクトの確認	代入と同時に値を表示	()	(x <- 3 +5)
	データの最初の部分を確認	head()	head(dataframe)
分析	主成分分析	prcomp()	pc.df <- prcomp(df, scale=TRUE)
	主成分分析結果の表示	summary()	summary(pc.df)
	クラスター分析	hclust()	x <- hclust(x, method="ward") plot(x)
	クラスター分析のデンドログラムを資格で囲む	rect.hclust()	rect.hclust(x, k=i, border=rainbow(i))
グラフ描写	描写画面の分割	par()	par(mfcol=c(2,2))
	レーダーチャートの描写	radarchart()	radarchart(rbind(max=a, min=b), x, maxmin=TRUE)
出力	画像としてグラフを保存	dev.copy()	dev.copy(png,file=paste("out/", i, ".png", sep="")) dev.off()
	テキストファイルに出力	print() sink()	sink(paste("out/", dir, "Summary_", case,".txt", sep="")) print(summary(x)) sink()
	テキストファイルにデータフレームの出力	write.table()	write.table(x, paste("out/", i,".txt", sep=""), append=FALSE, quote=FALSE,row.names=FALSE, sep="\t")
その他	ヘルプを表示	help()	help()

複数のデータを一つのグラフに重ねて描写する場合、軸は繰り返し描写されないよう、最後のデータ描写時に描写する。

for (i in 1:n-1){
    plot(x[i], type="l", col=rainbow(n)[i], ylim=c(a, b), xlab="", ylab="", axes=FALSE)
    par(new="TRUE")
}
plot(x[n], type="l", col=rainbow(n)[n], ylim=c(a, b), xlab="day", ylab="price")

解析を行う場合、その解析が適応できるデータ分布となっているか事前に判定すること。(正規分布を仮定している解析ではないか、独立標本を仮定していないか、など)

SAGA GIS

QGISをインストールしたら、SAGA GISも同時にインストールされていました。以前はなかった気がするので、スタンドアロンのインストーラではなく、ネットワークインストーラを利用したからでしょうか。

大学ではGRASSを教えてもらったのですが、GRASSはマップセットの設定とかがうまく理解できていないため、GUIで何とかなりそうなSAGA GISでまずはいろいろやってみようかと思います。

QGISはUIが多言語化されているし、書籍やネット上の情報もほかのフリーのGISソフトと比べて入手しやすいですが、地形解析のプラグインや3D viewの機能があまりよくありません。プラグインは2.0になってから対応するものが減ったのでしょうか。

QGISのメニューの中に、GRASSのマップセットを開くメニューがあるので、QGISの中でGRASSと連携して解析できるのだろうけれど、ひとまずこちらで。

とりあえず、以下を目標に調べてみます。

• 鳥瞰図の作成
• 水系図の作成
• 特定の地点でのCatchment Areaの計算
• 降水量を与えた時の流量の推定

河川流量は、どういった流出モデルやパラメータを使うかによって全然違うだろうけれど、どうせ正確な雨量データなんて入手できないので、桁が違わなければいいかなと考えています。