Tugas Regresi Spasial Praktikum - Interpolasi Data Meuse

Persiapan

library(sp)
library(sf)
library(terra)
library(gstat)

data(meuse)
data(meuse.grid)

d <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)

d$vble <- log(d$zinc)

coop <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)

map <- st_union(st_buffer(coop, 20, endCapStyle = "SQUARE")) |> st_sf()

summary(d$vble)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.727   5.288   5.787   5.886   6.514   7.517

plot(st_geometry(map), col = "grey90", main = "Sebaran log(zinc)")
plot(d["vble"], add = TRUE, pch = 16)

Titik dengan log(zinc) tinggi mengumpul di dekat sungai, jadi data punya pola spasial dan layak diinterpolasi.

Metode Interpolasi

Closest Observation

Nilai prediksi diambil dari satu titik pengamatan terdekat (poligon Voronoi).

v <- terra::voronoi(vect(d), bnd = vect(map))
plot(v, main = "Diagram Voronoi")
points(vect(d), cex = 0.5)

p1 <- d$vble[st_nearest_feature(coop, d)]

Nearest Neighbors

Rata-rata dari k = 5 titik terdekat, bobot semua tetangga sama.

nn <- gstat(formula = vble ~ 1, locations = d, nmax = 5, set = list(idp = 0))
p2 <- predict(nn, coop)$var1.pred

## [inverse distance weighted interpolation]

IDW

Bobot berbanding terbalik dengan jarak (pangkat 2).

idw <- gstat(formula = vble ~ 1, locations = d, nmax = nrow(d), set = list(idp = 2))
p3 <- predict(idw, coop)$var1.pred

## [inverse distance weighted interpolation]

Kriging

Ordinary kriging dengan model variogram spherical.

vemp <- variogram(vble ~ 1, d)
fv <- fit.variogram(vemp, vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1))
fv

##   model      psill    range
## 1   Nug 0.05066017   0.0000
## 2   Sph 0.59060556 897.0044

plot(vemp, fv, cex = 1.5)

kr <- gstat(formula = vble ~ 1, data = d, model = fv)
p4 <- predict(kr, coop)$var1.pred

## [using ordinary kriging]

Interpretasi: nugget = 0.051, partial sill = 0.591, dan range = 897 meter. Artinya titik yang jaraknya lebih dari range sudah tidak berkorelasi spasial, dan semakin dekat dua titik semakin mirip nilai zinc-nya.

Ensemble

Rata-rata dari prediksi NN, IDW, dan Kriging (bobot sama).

p5 <- (p2 + p3 + p4) / 3

Peta Hasil Interpolasi

xy <- st_coordinates(coop)
pred <- data.frame(x = xy[, 1], y = xy[, 2],
                   Closest = p1, NN = p2, IDW = p3,
                   Kriging = p4, Ensemble = p5)

r <- rast(pred, type = "xyz", crs = "EPSG:28992")
plot(r, nc = 3, col = hcl.colors(30, "viridis"))

Interpretasi: Closest Observation hasilnya berpetak-petak karena tiap wilayah hanya mengikuti satu titik terdekat. NN, IDW, dan Kriging lebih halus, sedangkan Ensemble merupakan gabungan ketiganya. Semua metode menunjukkan zinc tinggi di sekitar sungai.

Cross Validation (5-fold)

Data dibagi jadi 5 bagian, tiap bagian bergantian jadi data uji dan sisanya data latih.

set.seed(123)
kf <- dismo::kfold(nrow(d), k = 5)

RMSE <- function(obs, pred) sqrt(mean((obs - pred)^2))

cv <- matrix(NA, 5, 5,
             dimnames = list(NULL, c("Closest", "NN", "IDW", "Kriging", "Ensemble")))

for (k in 1:5) {
  test  <- d[kf == k, ]
  train <- d[kf != k, ]

  # closest observation
  q1 <- train$vble[st_nearest_feature(test, train)]

  # nearest neighbors
  g2 <- gstat(formula = vble ~ 1, locations = train, nmax = 5, set = list(idp = 0))
  q2 <- predict(g2, test)$var1.pred

  # idw
  g3 <- gstat(formula = vble ~ 1, locations = train, nmax = nrow(train), set = list(idp = 2))
  q3 <- predict(g3, test)$var1.pred

  # kriging (variogram di-fit ulang pakai data latih)
  fvk <- fit.variogram(variogram(vble ~ 1, train),
                       vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1))
  g4 <- gstat(formula = vble ~ 1, data = train, model = fvk)
  q4 <- predict(g4, test)$var1.pred

  # ensemble
  q5 <- (q2 + q3 + q4) / 3

  cv[k, ] <- c(RMSE(test$vble, q1), RMSE(test$vble, q2), RMSE(test$vble, q3),
               RMSE(test$vble, q4), RMSE(test$vble, q5))
}

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]

round(cv, 4)

##      Closest     NN    IDW Kriging Ensemble
## [1,]  0.6327 0.4958 0.4734  0.4281   0.4429
## [2,]  0.5987 0.4943 0.4717  0.4183   0.4368
## [3,]  0.6328 0.6516 0.6576  0.5055   0.5943
## [4,]  0.4874 0.5316 0.5158  0.4198   0.4716
## [5,]  0.5593 0.4791 0.4913  0.3840   0.4342

hasil <- data.frame(Metode = colnames(cv), RMSE = round(colMeans(cv), 4))
hasil <- hasil[order(hasil$RMSE), ]
hasil

##            Metode   RMSE
## Kriging   Kriging 0.4312
## Ensemble Ensemble 0.4760
## IDW           IDW 0.5219
## NN             NN 0.5305
## Closest   Closest 0.5822

barplot(hasil$RMSE, names.arg = hasil$Metode, las = 2, col = "steelblue",
        ylab = "RMSE", main = "Perbandingan RMSE (5-fold CV)")

Interpretasi: semakin kecil RMSE, semakin akurat metodenya. Metode terbaik adalah Kriging (RMSE = 0.4312), sedangkan yang terburuk adalah Closest (RMSE = 0.5822).

Kesimpulan

Berdasarkan RMSE hasil cross validation, Kriging paling baik untuk menginterpolasi zinc pada data meuse.