2014年12月4日木曜日

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で見ると、1行に南から北へデータが格納されているが、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)
}