はじめに
点データから、内挿補間を行い、面データやコンター（等値線）を作るという作業、特にクリギングという手法が昔から好きでした。どういうわけか、最近また急に、クリギングをやってみたくなりました（別に必要に迫られたわけではないのですが）。
過去、R 上でいろいろと試行錯誤していた結果、automap ライブラリに入っている autoKrige() を使うと、簡単にクリギングをすることができたので、使い方を記事に残します。今回は、座標値 (x, y) と、目的の数値 (value) というシンプルな CSV データをクリギングで補間する一連のやり方となります。
なお、クリギングの概要は、昔から Esri さんの以下のサイトがわかりやすいです。
データの準備等
パッケージの読み込み
最初に R で automap パッケージを読み込んでおきます。
library(automap)
サンプルデータ
サンプルデータとして、東京周辺の標高を使ってクリギングしてみます。標高値と経緯度は、地理院地図から取得します。
また、経緯度で表示すると、同じ経度差であっても、緯度によって実距離が異なりますので、クリギングをする場合は、平面直交座標系等に置き換えるとよいです。今回は、QGIS で経緯度を平面直交座標系9系へ変換しました。直感で操作していたら、なかなか思った通りに変換できずに大変でした。
平面直交座標系の場合、X 軸は真北が正、Y 軸は真東が正なので注意が必要です。（ただし、QGIS で表示や処理をしたら、平面直交座標系とは、X 座標と Y が逆転して表示されていることに気が付きました。）
サンプルデータは、平面直交座標系の X 軸と Y 軸を入れ替えて、x が平面直交座標系の Y 座標、y が 平面直交座標系の X 座標、value が標高値となるような CSV にして、次のステップの入力とします。名前は、test.csv としましょう。ヘッダの名前も後のコードで参照するため、変えないでください。
x,y,value
-4490.344528,-28548.142,1.72
268.1160572,-30613.14644,0
-5503.139564,-31568.54818,18.45
-8764.076397,-32155.14559,5.59
（以下省略）
データの読み込み
上記で作成したデータを読み込み、SpatialPointsDataFrame へと変換しておきます。
# データ読み込み
df <- read.csv("test.csv", header=TRUE)
dat.k <- as.data.frame(df)
coordinates(dat.k) = ~x+y #SpatialPointsDataFrame 化
# 確認
head(dat.k)
plot(dat.k)
グリッドの作成
次に、補間するポイント群となるグリッドデータを作成します。このグリッド上の各格子点に kriging で補間されたデータが格納されるようです。読み込んだデータの x, y のそれぞれの最大値と最小値をもとに、seq() でその間を100分割して、100×100のグリッドとします。最後に、gridded() で、x と y の組み合わせでグリッドを作成します。
coord <- coordinates(dat.k)
x.grid <- seq(min(coord[,1]), max(coord[,1]), length=100)
y.grid <- seq(min(coord[,2]), max(coord[,2]), length=100)
grid <- expand.grid(x.grid, y.grid)
colnames(grid) <- c('x','y')
gridded(grid) = ~x+y
クリギングの実行
autoKrige について
クリギングの際に使うことになる autoKrige() の使い方を簡単に説明します。
autoKrige(formula, data, grid, model=c(候補1, 候補2, ...), verbose=true)
-
formula：モデル式。普通クリギングはvalue~1、普遍クリギングはvalue~x+y。（value, x, y は読み込んだデータのカラム名に依存。） -
data：上記で CSV から 読み込んだデータ。SpatialPointsDataFrame 形式で。 -
grid：上記で作成した補間地点を示すグリッド。 -
model：バリオグラムにフィットさせるためのモデルの候補。例えば、球面（Sph）, 指数（Exp）, ガウス（Gau）, 線形（Lin）など。 -
verbose：TRUE にすると、詳細情報が出力される。
※他にもとれる引数はありますので、以下のマニュアルでご確認ください。
普通クリギング/Ordinary Kriging
まずは、普通クリギングを以下の通り実行してみます。モデル式は value~1、モデル候補は、球面（Sph）, 指数（Exp）, ガウス（Gau）, 線形（Lin） から最も良いものを選択します。
kriging_o <- autoKrige(value~1, dat.k, grid, model=c('Sph', 'Exp', 'Gau', 'Lin'), verbose=TRUE)
parameters_o <- kriging_o$var_model
以下の通り、今回は Gau モデルの平方和 (the sums of squares, SSerror) が最も小さいので、これが最も良いと判定されました。
Selected:
model psill range
1 Nug 32.24346 0.000
2 Gau 76.93488 4123.647
Tested models, best first:
Tested.models kappa SSerror
3 Gau 0 0.01213574
1 Sph 0 0.01253477
2 Exp 0 0.01285483
4 Lin 0 0.02059230
[using ordinary kriging]
次に、autoKrige() の結果をそのまま plot() に渡すと、補間値の様子、標準誤差の様子、そしてセミバリオグラムと選択されたモデルの様子が表示されます。
plot(kriging_o)
さらに、クリギング結果上に、元のデータの位置とコンターを表示してみます。ここでは、raster() 関数を使うために、raster ライブラリも読み込んでいます。
library(raster)
krig_o <- kriging_o$krige_output
r_o <- raster(krig_o['var1.pred'])
krig_o@proj4string <- CRS('+proj=tmerc +lat_0=36 +lon_0=139.833333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs')
plot(r_o, main="Ordinary kriging")
contour(r_o, col='white', add=T)
points(dat.k[,1], dat.k[,2])
普遍クリギング/Universal Kriging
次に、普通クリギングを同じように通り実行してみます。普通クリギングのモデル式は、value~x+y です。
kriging_u <- autoKrige(value~x+y, dat.k, grid, model=c('Sph', 'Exp', 'Gau', 'Lin'), verbose=TRUE)
parameters_u <- kriging_u$var_model
以下の通り、普遍クリギングでも Gau モデルが最も良いと判定されました。
Selected:
model psill range
1 Nug 26.76034 0.000
2 Gau 36.34114 2902.752
Tested models, best first:
Tested.models kappa SSerror
3 Gau 0 0.006642060
1 Sph 0 0.006769627
2 Exp 0 0.007182246
4 Lin 0 0.009727579
[using universal kriging]
補間値の様子、標準誤差の様子、そしてセミバリオグラムを表示してみます。
plot(kriging_u)
クリギング結果上に元のデータの位置とコンターを表示します。
krig_u <- kriging_u$krige_output
krig_u@proj4string <- CRS('+proj=tmerc +lat_0=36 +lon_0=139.833333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs')
r_u <- raster(krig_u['var1.pred'])
plot(r_u, main="Universal kriging")
contour(r_u, col='white', add=T)
points(dat.k[,1], dat.k[,2])
交差検証
automap パッケージには、交差検証を行うための関数 autoKrige.cv() が準備されています。これを用いて、普通クリギングと普遍クリギングを比較してみます。すでに、どちらも Gau モデルのフィッティングが良さそうとわかっていますので、モデルには Gau のみ指定します。
cv.o <- autoKrige.cv(value~1, dat.k, model = c('Gau'), verbose=c(TRUE,TRUE))
cv.u <- autoKrige.cv(value~x+y, dat.k, model = c('Gau'), verbose=c(TRUE,TRUE))
複数の交差検証の結果を比較するのに便利な関数 compare.cv() が用意されているので、これに上記の結果を渡します。設定で、bubbleplots=TRUE とすることで、各点の残渣を表示させることができます。
compare.cv(cv.o, cv.u, bubbleplots=TRUE)
結果は以下の通りです。
cv.o cv.u
mean_error -0.02803 -0.1119
me_mean -0.003054 -0.0122
MAE 5.301 4.914
MSE 50.87 39.71
MSNE 0.9934 0.8619
cor_obspred 0.7899 0.8396
cor_predres 0.08478 -0.04251
RMSE 7.132 6.302
RMSE_sd 0.6048 0.5344
URMSE 7.132 6.301
iqr 6.123 7.613
荒木ほか (2011) によれば、
- 平均誤差（mean error, 0 に近いほど良い，以下ME）
- 平均二乗平方根誤差（Root mean squared error, 小さいほど良い，以下RMSE)
- 標準化誤差（mean squared normalized error, 1に近いほど良い，以下MSNE)
をモデル選択の精度評価指標として利用しています。この３指標（mean_error、RMSE、MSNE）の結果を見ると、RMSE は逆転していますが、そのほかは、普通クリギングの方が良さそうです。
クリギング結果を出力して GIS で使う
せっかく作成したクリギング結果を他の GIS でも使うために、ファイル出力する方法も紹介します。
GeoTiff で出力
上記のコンター表示の際に作成したラスタデータを出力します。以下の関数は raster パッケージに入っているものです。冗長になるので、普通クリギングの例だけ示しますが、普遍クリギングも全く同じです。raster パッケージの詳細は以下の通りです。
まず、目的に応じて projectRaster() で投影変換します。
# 普通クリギング
r_o_ll <- projectRaster(r_o, crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
その後、writeRaster() でラスタデータを出力します。ここでは、GeoTiff で出力しています。
# 普通クリギング
writeRaster(r_o_ll, "ordinaly_kriging_result.tif", "GTiff")
QGIS で表示させた結果は以下の通りです。背景地図は地理院タイル（淡色地図）を使っています。
CSV で出力
以下のように、補間値をグリッドの座標値とともに出力します。ここでも普通クリギングの例だけ示しますが、変数が変わるだけで普遍クリギングも全く同じです。
# 普通クリギング
p_o <- krig_o@data$var1.pred
gc <- krig_o@coords
res_o <- cbind(gc, p_o)
colnames(res_o) <- c("x", "y", "pred")
write.csv(res_o, "ordinary_krige.csv", row.names=FALSE)
以下は、出力した CSV を QGIS で読み込んだ例です。背景地図は同じく地理院タイル（淡色地図）です。
地理院地図で実際の標高と比べてみます。今回の地形は、崖を挟んで低地と台地が広がり、標高 5-20 m にかけての領域が極めて少ないです。高低差の傾向はきちんと示せているかと思いますが、台地から低地へだいぶなだらかに補間されています。崖の部分のサンプリングが少ないため、崖の予想はなかなか難かったと思います。
最後に
これまでいろいろと試行錯誤に時間をかけましたが、思い起こすと、結局クリギングを実務で使ったことはないです。重要なのは元となる点データであって、そこを集めないとクリギングまでたどり着きません。やっぱり、足を使ってデータを集めるのが一番大変です。