結局、振り出しに戻った解析作業。
　とりあえず、データに戻って、ぶち込む変数を再度計算させようとしているわけですが、それにはラスタデータの集計が必要です。

　といっても、あんな腐れ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さんに泣きついたもんですが、今回は頼れるのは自分だけ。ま、メモリたらんくなったら、分散処理にするとか、メッシュサイズ大きくするとか、対策はある。

　ので、細かいことはそのとき考えるか。