R の automap パッケージでクリギング

## はじめに

点データから、内挿補間を行い、面データやコンター（等値線）を作るという作業、特にクリギングという手法が昔から好きでした。どういうわけか、最近また急に、クリギングをやってみたくなりました（別に必要に迫られたわけではないのですが）。

過去、[R](https://www.r-project.org/) 上でいろいろと試行錯誤していた結果、[automap](https://cran.r-project.org/web/packages/automap/index.html) ライブラリに入っている `autoKrige()` を使うと、簡単にクリギングをすることができたので、使い方を記事に残します。今回は、座標値 (x, y) と、目的の数値 (value) というシンプルな CSV データをクリギングで補間する一連のやり方となります。

なお、クリギングの概要は、昔から Esri さんの以下のサイトがわかりやすいです。

https://desktop.arcgis.com/ja/arcmap/10.5/tools/spatial-analyst-toolbox/how-kriging-works.htm

## データの準備等

### パッケージの読み込み

最初に R で automap パッケージを読み込んでおきます。

```R
library(automap)
```

### サンプルデータ

サンプルデータとして、東京周辺の標高を使ってクリギングしてみます。標高値と経緯度は、[地理院地図](https://maps.gsi.go.jp/)から取得します。

また、経緯度で表示すると、同じ経度差であっても、緯度によって実距離が異なりますので、クリギングをする場合は、平面直交座標系等に置き換えるとよいです。今回は、[QGIS](https://qgis.org/ja/site/about/index.html) で経緯度を[平面直交座標系9系](https://www.gsi.go.jp/sokuchikijun/jpc.html)へ変換しました。直感で操作していたら、なかなか思った通りに変換できずに大変でした。

:::note warn
平面直交座標系の場合、X 軸は真北が正、Y 軸は真東が正なので注意が必要です。（ただし、QGIS で表示や処理をしたら、平面直交座標系とは、X 座標と Y が逆転して表示されていることに気が付きました。）
:::

サンプルデータは、平面直交座標系の X 軸と Y 軸を入れ替えて、x が平面直交座標系の Y 座標、y が 平面直交座標系の X 座標、value が標高値となるような CSV にして、次のステップの入力とします。名前は、test.csv としましょう。**ヘッダの名前も後のコードで参照するため、変えないで**ください。

```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 へと変換しておきます。

```R
# データ読み込み
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 の組み合わせでグリッドを作成します。

```R
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()` の使い方を簡単に説明します。

```R
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 にすると、詳細情報が出力される。

※他にもとれる引数はありますので、以下のマニュアルでご確認ください。

https://cran.r-project.org/web/packages/automap/automap.pdf

### 普通クリギング/Ordinary Kriging

まずは、普通クリギングを以下の通り実行してみます。モデル式は `value~1`、モデル候補は、球面（`Sph`）, 指数（`Exp`）, ガウス（`Gau`）, 線形（`Lin`） から最も良いものを選択します。

```R
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)
```

![ordinary-kriging](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/b2bac082-9d7b-e62d-48d3-912ee7e59782.png)

さらに、クリギング結果上に、元のデータの位置とコンターを表示してみます。ここでは、`raster()` 関数を使うために、[raster](https://cran.r-project.org/web/packages/raster/index.html) ライブラリも読み込んでいます。

```R
library(raster)
```

```R
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])
```

![Ordinary Kriging Contour](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/a0b6ba2e-da65-dec3-0471-69cbb0fcf62a.png)

### 普遍クリギング/Universal Kriging

次に、普通クリギングを同じように通り実行してみます。普通クリギングのモデル式は、`value~x+y` です。

```R
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]
```

補間値の様子、標準誤差の様子、そしてセミバリオグラムを表示してみます。

```R
plot(kriging_u)
```

![universal-kriging](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/a581c267-5802-6c38-4cfd-5cd6b0563c48.png)


クリギング結果上に元のデータの位置とコンターを表示します。

```R
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])
```

![Universal Kriging Contour](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/5635bb4f-3ff6-abb2-202e-53acd8742b67.png)

### 交差検証

automap パッケージには、交差検証を行うための関数 `autoKrige.cv()` が準備されています。これを用いて、普通クリギングと普遍クリギングを比較してみます。すでに、どちらも Gau モデルのフィッティングが良さそうとわかっていますので、モデルには `Gau` のみ指定します。

```R
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` とすることで、各点の残渣を表示させることができます。

```R
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
```

![kriging cross validation](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/418ee8b1-14a8-ec22-8fec-3a8ab50187ee.png)


[荒木ほか (2011)](https://www.jstage.jst.go.jp/article/taiki/46/5/46_241/_pdf) によれば、

* 平均誤差（mean error, 0 に近いほど良い，以下ME）
* 平均二乗平方根誤差（Root mean squared error, 小さいほど良い，以下RMSE)
* 標準化誤差（mean squared normalized error, 1に近いほど良い，以下MSNE)

をモデル選択の精度評価指標として利用しています。この３指標（mean_error、RMSE、MSNE）の結果を見ると、RMSE は逆転していますが、そのほかは、普通クリギングの方が良さそうです。

:::note warn
出力結果の MSNE は、[マニュアル](https://cran.r-project.org/web/packages/automap/automap.pdf)では、"Mean Squared Normalized Error, mean of the squared z-scores. Ideally small." とあります。一方、`autoKrige.cv()` は内部的に gstat パッケージの `krige.cv()` を使っているようですが、その[マニュアル](https://cran.r-project.org/web/packages/gstat/gstat.pdf)を見ると、MSNE は "Mean square normalized error, ideally close to 1" とあります。
:::

## クリギング結果を出力して GIS で使う

せっかく作成したクリギング結果を他の GIS でも使うために、ファイル出力する方法も紹介します。

### GeoTiff で出力

上記のコンター表示の際に作成したラスタデータを出力します。以下の関数は raster パッケージに入っているものです。冗長になるので、普通クリギングの例だけ示しますが、普遍クリギングも全く同じです。raster パッケージの詳細は以下の通りです。

https://cran.r-project.org/web/packages/raster/raster.pdf

まず、目的に応じて `projectRaster()` で投影変換します。
```R
# 普通クリギング
r_o_ll <- projectRaster(r_o, crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
```

その後、`writeRaster()` でラスタデータを出力します。ここでは、GeoTiff で出力しています。

```R
# 普通クリギング
writeRaster(r_o_ll, "ordinaly_kriging_result.tif", "GTiff")
```

QGIS で表示させた結果は以下の通りです。背景地図は[地理院タイル（淡色地図）](https://maps.gsi.go.jp/development/ichiran.html#pale)を使っています。
![ordinary kriging points plot on QGIS](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/88708b0c-65bc-20dc-89c2-f89b8315ec3f.png)

### CSV で出力

以下のように、補間値をグリッドの座標値とともに出力します。ここでも普通クリギングの例だけ示しますが、変数が変わるだけで普遍クリギングも全く同じです。

```R
# 普通クリギング
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 で読み込んだ例です。背景地図は同じく地理院タイル（淡色地図）です。

![ordinary kriging points plot on QGIS](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/930108c2-0e87-379d-254c-e4f225a51954.png)


地理院地図で実際の標高と比べてみます。今回の地形は、崖を挟んで低地と台地が広がり、標高 5-20 m にかけての領域が極めて少ないです。高低差の傾向はきちんと示せているかと思いますが、台地から低地へだいぶなだらかに補間されています。崖の部分のサンプリングが少ないため、崖の予想はなかなか難かったと思います。

![実際の標高](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/4ea9fc10-28f4-ab8c-a52c-018781722d1d.png)

## 最後に

これまでいろいろと試行錯誤に時間をかけましたが、思い起こすと、結局クリギングを実務で使ったことはないです。重要なのは元となる点データであって、そこを集めないとクリギングまでたどり着きません。やっぱり、足を使ってデータを集めるのが一番大変です。