二次元データの(線形)補間を行う
一次元のケースは二年も前に書いた
でよいとして、二次元のケースをなんとかしたい。したいというか、しなければならない状況になった。
ググりまくったところ
というパッケージに入っている関数を使えばよさげだということがわかった。
この２つのパッケージは、手元の補間したいデータの構造に応じて使い分ければよさそうだった。具体的には、以下のように適当にデータを作っておく。
library(dplyr) library(tidyr) #適当なサンプルデータ X0 <- 1:10 Y0 <- 1:10 df <- expand.grid(x=X0, y=Y0) %>% mutate(z=x^2-y^2+y) ls <- list(x=X0, y=Y0, z=matrix(df$z, nrow=length(X0), ncol=length(Y0)))
これは、dfがデータフレームで、lsがリストのデータ構造になっている、データ形式だけが違う同一のデータだ。中身の意味は見ればだいたいわかるように、x, yという二次元座標値に値zが割り振られている、そういうもんだ。
> df %>% head x y z 1 1 1 1 2 2 1 4 3 3 1 9 4 4 1 16 5 5 1 25 6 6 1 36 > ls %>% head $x [1] 1 2 3 4 5 6 7 8 9 10 $y [1] 1 2 3 4 5 6 7 8 9 10 $z [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 1 -1 -5 -11 -19 -29 -41 -55 -71 -89 [2,] 4 2 -2 -8 -16 -26 -38 -52 -68 -86 （... 途中省略 ...） [9,] 81 79 75 69 61 51 39 25 9 -9 [10,] 100 98 94 88 80 70 58 44 28 10
このような２つの異なるデータ形式の同一のデータがあり、それらの二次元補間をしようと思った場合には、
- データフレーム(data.frame)の場合はakimaパッケージ
- リスト(list)のの場合はfieldsパッケージ
という使い分けで補間してやるのがよさそうだった。実際は、データ形式を揃えればどちらか片方のパッケージだけでいいんだけども。rlistパッケージあたりを使えば簡単にデータの変換できるのか？
実際にたしかめるために、まず、補間する点(x, y)の組を以下のように与えておく。そして、各点に対する線形な補間値を計算させてみる。また、あわせて答えの一致も確認する。
#補間する点 x <- c(2, 4, 8) y <- c(3, 7, 2)
akimaパッケージの場合
akimaパッケージは不規則なグリッド*1上にあるデータのための補間パッケージだ。このパッケージで二次元線形補間をするためには
- interp関数
- inerpp関数
を用いる。interpp関数の方が"点ごと"の補間に対応して、interp関数は、下記コードのコメントに書いたようにx×yの直積集合の全ての元に対して補間を実行してくれる。
library("akima") #点(2,3)(4,7)(8,2)での線形補間 inter_pp <- interpp(df$x, df$y, df$z, x, y) #直積集合(2,4,8)×(3,7,2)での線形補間 inter_p <- interp( df$x, df$y, df$z, x, y)
結果はリストで帰ってきて、zに補間値が入っている。
> inter_pp$z [1] -2 -26 62 > inter_p$z [,1] [,2] [,3] [1,] -2 -38 2 [2,] 10 -26 14 [3,] 58 22 62
ちなみにx, yには補間した位置情報が入ってる。また、interp/interpp関数ともにlinear引数があって、これをFALSE(デフォルトはTRUE)にするとスプライン補間にかわるはずなんだけど、このデータだとなぜか全部NAになった、なぜ？？？関数を掘っていったら、"sdsf3p"なるFortranの関数がよばれていたので、これを読まないとだめだな・・・
fieldsパッケージの場合
fieldsパッケージは空間統計のためのツールが格納されたパッケージとのことだ。このパッケージでは
という２つの関数が2次元線形補間を実行してくれる関数だ。これもakimaパッケージのinterp & interppの関係に近くて、点の補間かグリッドまるごと(直積集合)の補間かで関数が異なっているだけ。
> library(fields) > #点(2,3)(4,7)(8,2)での線形補間 > interp.surface(ls, cbind(x, y)) [1] -2 -26 62 > #直積集合(2,4,8)×(3,7,2)での線形補間 > interp.surface.grid(ls, list(x=x, y=y)) $x [1] 2 4 8 $y [1] 3 7 2 $z [,1] [,2] [,3] [1,] -2 -38 2 [2,] 10 -26 14 [3,] 58 22 62
結果はakimaパッケージのものとちゃんと一致している。どちらのパッケージを使うのかは、データ形式に応じて使い分けるか、変換の手間を惜しまず、どちらか一方だけ使っててもいい感じかな。
*1:均一なグリッドじゃないという意味かと
