R Space interpolation 空間補間

前回、3次元データ(3列のデータ)の2次元カーネル推定を行ったが、別の手段を考える必要が出てきた。2次元のカーネル推定と言っても、

• 点(x,y)の分布から点の密度を滑らかな関数で推定する場合
• 点(x,y,z)の分布が与えられている場合に滑らかな分布の形状を推定する場合

では、Rの操作は異なる。今回は後者をやりたいのだが、ここに載っているkde2d関数を使う方法は実質zが整数値かつ（整形後データの規模が）Rで扱える程の小さい値であることが前提なので、今回はkde2d関数は使えないと考えた。

そもそも、こちらを見るとわかるように後者の場合は「カーネル推定」というよりは「空間補間」の問題であるのではないか。…というわけでRの関数を探した。krigeという関数が「クリギング」という手法で補間をやってくれる。

空間補間をするにはpackage 'sp' 'gstat' を使う。

まずは下準備。

library(gstat)
library(sp)

#適当に点(x,y,z)の組を発生させる
x <- seq(1,100,by=1)
y <- rnorm(100,50,20)
z <- rnorm(100)
dat <- matrix(cbind(x,y,z),ncol=3)
colnames(dat) <- c("x","y","z")
dat <- data.frame(dat)

#データ形式をSpatialPointsDataFrameに変換する
dat.coords <- cbind(dat$x,dat$y) #x,y座標をペアにする
dat <- SpatialPointsDataFrame(dat.coords,dat)

次に、バリオグラムを調べる。これは、全点について2点間距離を集計し、ある程度の区間毎にセミバリアンス(距離の2乗の平均のようなもの?)を計算したもの。全体の形状を近似するためのinputとするらしい。

#バリオグラムを図示(乱数なので、場合によって変わります）
dat.vario <- variogram(z~x+y,dat,cressie=T)
plot(dat.vario,pch=1,cex=1.2)

次に、上記を近似するモデルを作る。fit.variogram関数を使う。

#バリオグラムの近似曲線を決定
my.psill = 1.2
my.range = 50
my.nugget = 0.60
dat.fit = fit.variogram(dat.vario,vgm(model="Sph",psill=my.psill,range=my.range,nugget=my.nugget),fit.method=1)

dat.fit
f.psill = dat.fit$psill[2]
f.range = dat.fit$range[2]
f.nugget = dat.fit$psill[1]

plot(dat.vario,dat.fit,pch=20,cex=1.5,col=1)
m <- vgm(f.psill, "Sph", f.range, f.nugget)

以下のように近似される。正直微妙だが、とりあえず先に進む。

次に、空間補間の状況を図示するためにグリッド(描画時の升目)を用意する。

x.grid <- seq(1,100,by=1)
y.grid <- seq(1,100,by=1)
dat.grid <- matrix(c(x.grid,y.grid),ncol=2)
dat.grid <- data.frame(dat.grid)
colnames(dat.grid) <- c("x.grid","y.grid")
coordinates(dat.grid) = ~x.grid+y.grid
#データ型をSpatialPoints型に変換
dat.grid <- SpatialPoints(dat.grid, proj4string=CRS(as.character(NA)))
#データ型をSpatialPixel型に変換
dat.grid <- spsample(dat.grid,type="regular",cellsize=1,offset=c(1,1))
gridded(dat.grid) = TRUE

元のデータ(dat)、バリオグラムの近似モデル(m)、グリッド(dat.grid)を合わせて空間補間を行う。krige関数を用いる。

dat.krige <- krige(z~1, dat, dat.grid, model = m)
spplot(dat["z"]) #inputデータの分布図
spplot(dat.krige["var1.pred"]) #空間補間後

inputデータの分布はこちら。

【気づき】

inputデータにおける赤い点(z座標は小さい)の合間にある青い点(z座標は大きい)がうまく再現出来ていない。実際の地形のように、高低差の傾向が空間全体に一貫している場合に有効な方法であると認識した。とりあえず、GCMSのデータで使ってみるとどうなるか…、試してみる。