結局、振り出しに戻った解析作業。
とりあえず、データに戻って、ぶち込む変数を再度計算させようとしているわけですが、それにはラスタデータの集計が必要です。
といっても、あんな腐れArcGIS使ってたら精神衛生上悪くってしょうがない。
ということで、ラスタをテキストで吐き出して、そのデータでラスタ演算させるRスクリプト書いてみました。
やっていることは、ターゲットセルから円形のバッファを発生させて、バッファ内の特定セル値をカウントさせるという、ラスタ演算入門編にでも出てきそうなネタ。
んが、こんな単純処理でも腐れGISソフトってば、理不尽なエラー吐きまくりやがるのですな。多分。
それに比べれば、Rで書けば、ミスしても「なんて僕はお馬鹿さんなんだろう」と自責の念に駆られるだけですむので気が楽。
しかも、とても親切なエラーメッセージもついてくるので、見通しも立てやすいです。
ってわけで、書いたスクリプトはこんな感じ。
in.raster <- sample(1:10,1000000,replace=TRUE)
dim(in.raster) <- c(1000,1000)
buffer.length <- 5
x.square.buffer <- rep((buffer.length *-1):buffer.length,times=(buffer.length*2+1))
y.square.buffer <- rep(c((buffer.length *-1):buffer.length),each=11)
distance.buffer <- sqrt(x.square.buffer^2 + y.square.buffer^2)
fun <- function(x,y){
if(x <= buffer.length){
y
}
}
x.buffer <- unlist(mapply(fun,distance.buffer,x.square.buffer))
y.buffer <- unlist(mapply(fun,distance.buffer,y.square.buffer))
out.raster <- in.raster * 0
for(out.y in (buffer.length + 1):(nrow(in.raster) – buffer.length)){
for(out.x in (buffer.length + 1):(nrow(in.raster) – buffer.length)){
check.vector <- in.raster[cbind(x.buffer + out.x,y.buffer + out.y)]
out.raster[out.x,out.y] <-length(check.vector[check.vector == 5])
}
}
細かくチェックはしとらんが、大枠はいけてるはず。
今回は、シンプルを心がけてみました。というのも、今回の処理は結構時間がかかりそうなので、無駄なループとかは極力排除してベクトル演算中心のコーディングです。
最初は、apply系をメイン処理に使おうと考えていたのですが、どうやらapply系は速度が遅いんだそうな。
perlのmapとかもそうみたいだけど、forループのラッパーな関数って、総じて遅いってことなんだろうか。
ちなみに、apply系関数の源流ってlispにあるみたいですな。つってもlispなんてエディタいぢるときくらいにしか見ないけど。
たしかに、Cとかlispとかいろんな言語の影響を受けているであろうRの言語仕様って、非常に複雑なうえに構文糖的な仕様もあるんだけど、perlなんかと比べても見やすい気が立てやすい・・・気がするだけかも。あ、でも結構過去の遺産をうまく使えてたりするんですよね。
さあ、後はほんちゃんのデータにこのスクリプトを突っ込むだけだ。ターゲットは6億セルだが、どんなもんだか。
昔、perlで同じようなスクリプト書いたときは、謎のメモリリーク起こして、ODさんに泣きついたもんですが、今回は頼れるのは自分だけ。ま、メモリたらんくなったら、分散処理にするとか、メッシュサイズ大きくするとか、対策はある。
ので、細かいことはそのとき考えるか。