2014年11月21日金曜日

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))

2014年11月14日金曜日

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")
}

2014年11月7日金曜日

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))

#エクスポート(略)