時系列ラスターデータの処理について
この記事はFOSS4G Advent Calendar 2013に参加しています。
普段は自分用メモブログなので今回も雑になる可能性があります。。
GISはいうまでもなく空間データ処理にフォーカスがあたってますが、時系列ものにはめっぽう弱い印象です。唯一救いがあるのがGRASSで、時系列データを番号順に並べれて処理すれば時系列処理ができます。とはいえ処理ができるのはr.mapcalcのラスター演算とr.seriesというコマンドくらいで、時系列データの回帰分析やnull値カウントくらいしかできません。GRASS7ではt.ではじまるTemporal GIS modulesが追加されるので期待大です。
ほかになんかないですかね？というかなんでもできちゃうRでは処理できないですかねぇ？とspgrass6 libraryでGRASSデータをRにインポートするもメモリ圧迫するので数百枚の時系列ラスターデータは扱いきれません。
というわけで、なんかないかいなと調べていたことろ、ndvitsというRライブラリがtiff形式の時系列ラスターデータをRにdataframeとしてぶちこめるとのことで試してみました。これがなかなかよく、通常メモリがパンクするくらいのデータもGDALを介してRに読み込めます。
今回はGIMMS1981〜2006年のNDVI（正規化植生指数：植生量の多いところでは１に近い値を取り、少ないところでは0に近い値を取る）データをRに取り込んでみます。
データのインポートには、時系列データが記載されたフルパスのリストと対象地域のshpファイル（kmlでもよいみたい）が必要です。
今回はlist.txtの中に、
600
/media/3TB/AVHRR/EA_tif/GIMMS_1982_1_a.tif
/media/3TB/AVHRR/EA_tif/GIMMS_1982_1_b.tif
/media/3TB/AVHRR/EA_tif/GIMMS_1982_2_a.tif
/media/3TB/AVHRR/EA_tif/GIMMS_1982_2_b.tif
/media/3TB/AVHRR/EA_tif/GIMMS_1982_3_a.tif
・
・
・
・
/media/3TB/AVHRR/EA_tif/GIMMS_2006_11_a.tif
/media/3TB/AVHRR/EA_tif/GIMMS_2006_11_b.tif
/media/3TB/AVHRR/EA_tif/GIMMS_2006_12_a.tif
/media/3TB/AVHRR/EA_tif/GIMMS_2006_12_b.tif
のように、1行目にファイル数、2行目からは各フルパスを記載していきます。この順番でインポートされていくので、時系列順に並べる必要があります。
また、GIMMSデータはアジア全域のものですが、おいおいゴビ砂漠周辺を分析したいので中国・モンゴルを対象とします。境界データをGADMより入手しておきます。QGISなので一つのshpファイルに結合し、China_Mon.shpとしておきます。
ndvitsパッケージではExtractFileという関数を使用します（GIMMSにはExtractGIMMSという特別な関数もありますが汎用性をもたすためにこっちを使います）。
が、
肝心なこの関数、Ubuntuで実行する際エラーが出ます。。（Windows・Macは未確認）
そのため、１行だけ修正を加えました。
修正したものをこちらにアップしたのでお使いくださいませ
https://app.box.com/s/wkfw9zgb7a7aull6t7cs
#Rの起動
R
#ndvitsのインストール
install.packages("ndvits")
library(ndvits)
#実行前にもろもろの設定をしておきます。
#shpファイルの設定
shapefile="China_Mon"
shapedir="/media/3TB/AVHRR/CHINA_MON/SHP/" #最後の/を忘れずに！
#リストファイルの設定
listfile="/media/3TB/AVHRR/EA/list.txt"
#外部出力ファイルの設定
outfile=outfile.txt
#観測頻度の設定。ここでは月２回の観測なので年２４回
period=24
#ExtractFileならぬRxtractFile_Editを読み込み、実行します。
Source(”ExtractFileEdit.R")
GIMMS=ExtractFile_Edit(shapefile,shapedir,listfile,outfile,period,shapeext="shp")
すぐにエラーが出ていなければうまくいっています。ちなみにラスターデータの領域はshpファイルの領域に対して余裕を持たせるのがコツです。
しばらくまつと・・・
class(GIMMS)
[1] "data.frame"
となり、インポート完了です。
ちなみに、中身をみてみると
names(GIMMS)
[1] "xa" "ya" "ID_0" "ISO"
[5] "NAME_0" "ID_1" "NAME_1" "ID_2"
[9] "NAME_2" "NL_NAME_2" "VARNAME_2" "TYPE_2"
[13] "ENGTYPE_2" "ID_3" "ISO_2" "NAME_3"
[17] "ID_3_1" "NAME_3_1" "NL_NAME_1" "VARNAME_1"
[21] "TYPE_1" "ENGTYPE_1" "V23" "V24"
[25] "V25" "V26" "V27" "V28"
[29] "V29" "V30" "V31" "V32"
[33] "V33" "V34" "V35" "V36"
・
・
・
[609] "V609" "V610" "V611" "V612"
[613] "V613" "V614" "V615" "V616"
[617] "V617" "V618" "V619" "V620"
[621] "V621" "V622"
と、あたまの22列はshpの属性テーブルから引き出していることがわかります。
実際のGIMMSデータはのこりのV23〜V622を使えばよいということになります。
あとはRお得意の処理をかければよいだけ！
以上、時系列ラスターデータの処理についてでした。
ご指摘・改善点などありましたらご連絡くださいませ。
