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

0 件のコメント:

コメントを投稿