Hyetgraph TRMM 3B42 with R

RにてTRMMの3B42ファイルから、特定の地点のハイエトグラフを作成しました。
http://disc.sci.gsfc.nasa.gov/giovanniからダウンロードできる3B42はHDF4なので、h4toh5convertでHDF5に変換します。
> for %i in (*.hdf) do h4toh5convert %i
以降Rで作業します。
HDFファイルを扱うために、rhdf5ライブラリを利用します。rhdf5パッケージについては、以下ドキュメントがあります。
http://www.bioconductor.org/packages/release/bioc/vignettes/rhdf5/inst/doc/rhdf5.pdf

> library(rhdf5)

hdfファイルの構造を確認します。
> h5ls("3B42_daily.2013.08.27.7.G3.h5")

TRMM 3B42のHDFファイルには、西経180°南緯50°から、西経180°北緯50°へ、続いて西経179.75°南緯50°から北へ、といった順で、経緯0.25°ずつデータが格納されています。そのため、取得したい地点の緯度経度から、データの位置を求めます。

> origin.lon <- -180
> origin.lat <- -50
> resolution.lon <- 0.25
> resolution.lat <- 0.25

> in.data.dir             <- "in/3B42_daily/"
> in.data.name.common     <- "3B42_daily."
> in.data.name.common2    <- ".7.G3.h5"
> out.data.dir            <- "out/3B42_hyet/"
> in.data.name.hdf        <- "precipitation"

ハイエトグラフを作成したい地点の緯度経度のリストをテキストで準備しておきます。
> in.coodination.filename <- "in/3B42_daily/hyet_coodination.csv"

対象となる期間を指定します。
> target.period.year.begin <- 2012
> target.period.year.end   <- 2014
> target.period.begin.date <- "06-01"
> target.period.end.date   <- "09-30"

サイクルの日数の計算
> dulation <- as.numeric(as.Date(paste(target.period.year.begin,"-",target.period.end.date, sep="")) - as.Date(paste(target.period.year.begin,"-",target.period.begin.date, sep="")) + 1 )
> target.period.years <- target.period.year.end - target.period.year.begin + 1

期間の日付を作成
> target.date <- matrix(data = , nrow = dulation, ncol = target.period.years)
> for(y in 1:target.period.years){
>   day.begin <-  as.Date(paste(target.period.year.begin + y -1,"-",target.period.begin.date, sep=""))
>   day.end   <-  as.Date(paste(target.period.year.begin + y -1,"-",target.period.end.date, sep=""))
>   target.date[,y] <- day.begin:day.end
>   }

期間のファイル名を作成
>target.filenames <- matrix(paste(in.data.dir, in.data.name.common, format(as.Date(target.date, origin="1970/1/1"), "%Y.%m.%d"),
>                                 in.data.name.common2, sep=""),  nrow = dulation, ncol = target.period.years)

緯度経度のリストを読み込み
> location.table <- read.table(in.coodination.filename, header = TRUE, sep = ",")

データの位置の設定
> tada.matrix.row <- (location.table$Lon - origin.lon) %/% resolution.lon + 1
> tada.matrix.col <- (location.table$Lat - origin.lat) %/% resolution.lat + 1

ファイルの読み込み
> precipitation.table <- data.frame(location.table)
> for (y in 1:target.period.years){
>   for (d in 1:dulation){
>     precipitation.day <- h5read(target.filenames[d,y], name = in.data.name.hdf)[matrix(c(tada.matrix.row,tada.matrix.col),ncol = 2)]
>     precipitation.table  <- cbind(precipitation.table, precipitation.day)
>   }
> }
> for (y in 1:target.period.years){
>   precipitation.table <- cbind(precipitation.table, apply(precipitation.table[ncol(location.table)+dulation * (y - 1) + (1:dulation)],MARGIN=1,sum))
> }
> colnames(precipitation.table) <- c(colnames(location.table),format(as.Date(target.date, origin="1970/1/1"), "%Y-%m-%d"), target.period.year.begin:target.period.year.end)

グラフの描写及び書き出し
> for (p in 1:nrow(location.table)){
>   win.graph(width = 700, height = 500)
>   for (y in 1:target.period.years){
>     plot(as.Date(target.date[,y], origin="1970/1/1"), precipitation.table[ncol(location.table)+dulation * (y - 1) + (1:dulation)][p,],
>          type="l", lwd=2, ylim=c(0,100), xlab="Date", ylab="Precipitation (mm)", main=paste("Precipitation", location.table$SN[p], location.table$Name[p]),
>          col=rainbow(target.period.years)[y])
>     par(new="TRUE")
>   }
>   legend("topleft", legend=paste((target.period.year.begin:target.period.year.end),"Total",
>                                  sprintf("%5.1f",precipitation.table[p,as.character(target.period.year.begin:target.period.year.end)]), "mm"),
>          col=rainbow(target.period.years), lty=1, lwd=2, par(mar = c(1, 1, 1, 1)))
>   dev.copy(png, file=paste(out.data.dir,"plot_", location.table$SN[p], "_", location.table$Name[p], ".png", sep=""),width=700, height=500)
>   dev.off()
>   graphics.off()
> }

データの書き出し
> write.table(x = precipitation.table, paste(out.data.dir,"Precipitation.csv", sep=""), quote=FALSE,row.names=FALSE, sep="\t")

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

#エクスポート(略)

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

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