Interpolasi Spasial

Interpolasi spasial adalah proses menduga nilai suatu peubah pada lokasi yang tidak tersampel berdasarkan nilai pada lokasi-lokasi yang telah diamati. Pada tugas ini digunakan data meuse (paket sp), yaitu data kandungan logam berat tanah di dataran banjir Sungai Meuse, Belanda (155 titik pengamatan). Peubah yang diinterpolasi adalah log(zinc), karena sebaran zinc menjulur ke kanan dan transformasi log membuat sebaran lebih simetris serta lebih sesuai untuk kriging.

Metode yang digunakan dibagi menjadi dua kelompok:

  • Metode Deterministik: Closest Observation, Inverse Distance Weight (IDW), Nearest Neighbors (NN), dan Pendekatan Ensembel.
  • Metode Stokastik: Kriging.

Kinerja seluruh metode dibandingkan menggunakan RMSE hasil k-fold Cross Validation.

library(sp)       # sumber data meuse & meuse.grid
library(sf)       # data spasial
library(gstat)    # idw, variogram, kriging
library(ggplot2)  # visualisasi
library(dplyr)
library(tidyr)
library(knitr)
data(meuse)
data(meuse.grid)

# Data titik pengamatan
meuse_sf <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992, remove = FALSE)
meuse_sf$log_zinc <- log(meuse_sf$zinc)

# Grid lokasi prediksi
grid_sf <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992, remove = FALSE)
coords  <- st_coordinates(grid_sf)

cat("Jumlah titik pengamatan:", nrow(meuse_sf), "\n")
## Jumlah titik pengamatan: 155
cat("Jumlah titik grid prediksi:", nrow(grid_sf), "\n")
## Jumlah titik grid prediksi: 3103
head(st_drop_geometry(meuse_sf)[, c("x", "y", "zinc", "log_zinc")])
summary(meuse_sf$log_zinc)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.727   5.288   5.787   5.886   6.514   7.517
hist_df <- st_drop_geometry(meuse_sf) |>
  select(zinc, log_zinc) |>
  pivot_longer(everything(), names_to = "Peubah", values_to = "nilai")

ggplot(hist_df, aes(nilai)) +
  geom_histogram(bins = 20, fill = "steelblue", color = "white") +
  facet_wrap(~ Peubah, scales = "free") +
  labs(title = "Sebaran zinc dan log(zinc)", x = "Nilai", y = "Frekuensi") +
  theme_minimal()

ggplot() +
  geom_tile(data = meuse.grid, aes(x, y), fill = "grey92") +
  geom_point(data = st_drop_geometry(meuse_sf),
             aes(x, y, color = log_zinc, size = log_zinc)) +
  scale_color_viridis_c(option = "plasma") +
  coord_equal() +
  labs(title = "Sebaran Titik Pengamatan log(zinc)", x = "X", y = "Y") +
  theme_minimal()

Fungsi bantu dan pembagian fold untuk cross validation (fold yang sama dipakai oleh semua metode agar perbandingan adil):

rmse <- function(obs, pred) sqrt(mean((obs - pred)^2, na.rm = TRUE))

# fungsi untuk menggambar peta prediksi pada grid
peta_pred <- function(pred, judul, legend = "log(zinc)") {
  ggplot(data.frame(X = coords[, 1], Y = coords[, 2], z = pred),
         aes(X, Y, fill = z)) +
    geom_raster() +
    scale_fill_viridis_c(option = "plasma", name = legend) +
    coord_equal() +
    labs(title = judul, x = "X", y = "Y") +
    theme_minimal()
}

set.seed(2024)
K <- 10
folds <- sample(rep(1:K, length.out = nrow(meuse_sf)))
table(folds)
## folds
##  1  2  3  4  5  6  7  8  9 10 
## 16 16 16 16 16 15 15 15 15 15

Metode Deterministik

Metode deterministik menduga nilai pada lokasi baru langsung dari nilai dan jarak titik-titik pengamatan di sekitarnya, tanpa memodelkan struktur autokorelasi spasial dan tanpa memberikan ukuran ketidakpastian prediksi.

Closest Observation

Setiap titik grid diberi nilai dari satu pengamatan terdekat (nmax = 1). Hasilnya setara dengan poligon Thiessen/Voronoi.

pred_closest <- idw(log_zinc ~ 1, meuse_sf, grid_sf,
                    nmax = 1, idp = 0, debug.level = 0)

peta_pred(pred_closest$var1.pred, "Closest Observation")

Peta berbentuk blok-blok (poligon Thiessen) karena setiap sel hanya mengikuti satu titik terdekat, sehingga permukaan tidak halus dan sangat sensitif terhadap pengamatan tunggal.

Invers Distance Weight

IDW menduga nilai sebagai rata-rata terboboti dari seluruh pengamatan, dengan bobot berbanding terbalik terhadap jarak berpangkat idp:

\[\hat{z}(s_0) = \frac{\sum_{i=1}^{n} z(s_i)\, d_{0i}^{-p}}{\sum_{i=1}^{n} d_{0i}^{-p}}\]

Semakin besar p, semakin besar pengaruh titik yang dekat. Nilai p terbaik dicari dengan 10-fold cross validation.

idp_grid <- seq(0.5, 5, by = 0.5)

rmse_idp <- sapply(idp_grid, function(p) {
  pred <- rep(NA_real_, nrow(meuse_sf))
  for (i in 1:K) {
    tr <- meuse_sf[folds != i, ]
    te <- meuse_sf[folds == i, ]
    pred[folds == i] <- idw(log_zinc ~ 1, tr, te, idp = p, debug.level = 0)$var1.pred
  }
  rmse(meuse_sf$log_zinc, pred)
})

tab_idp <- data.frame(idp = idp_grid, RMSE = round(rmse_idp, 4))
idp_opt <- idp_grid[which.min(rmse_idp)]
kable(tab_idp, caption = "RMSE IDW untuk berbagai nilai idp (10-fold CV)")
RMSE IDW untuk berbagai nilai idp (10-fold CV)
idp RMSE
0.5 0.7033
1.0 0.6495
1.5 0.5822
2.0 0.5233
2.5 0.4853
3.0 0.4667
3.5 0.4610
4.0 0.4619
4.5 0.4656
5.0 0.4703 
cat("idp optimum:", idp_opt, "\n")
## idp optimum: 3.5
ggplot(tab_idp, aes(idp, RMSE)) +
  geom_line() + geom_point() +
  geom_point(data = tab_idp[tab_idp$idp == idp_opt, ], color = "red", size = 3) +
  labs(title = "Pemilihan idp pada IDW", x = "Pangkat jarak (idp)", y = "RMSE CV") +
  theme_minimal()

pred_idw <- idw(log_zinc ~ 1, meuse_sf, grid_sf, idp = idp_opt, debug.level = 0)

peta_pred(pred_idw$var1.pred, paste0("IDW (idp = ", idp_opt, ")"))

Nearest Neighbors

Nilai prediksi adalah rata-rata sederhana dari k pengamatan terdekat (idp = 0, semua tetangga berbobot sama). Jumlah tetangga k terbaik dicari dengan 10-fold cross validation.

k_grid <- 2:15

rmse_k <- sapply(k_grid, function(k) {
  pred <- rep(NA_real_, nrow(meuse_sf))
  for (i in 1:K) {
    tr <- meuse_sf[folds != i, ]
    te <- meuse_sf[folds == i, ]
    pred[folds == i] <- idw(log_zinc ~ 1, tr, te, nmax = k, idp = 0, debug.level = 0)$var1.pred
  }
  rmse(meuse_sf$log_zinc, pred)
})

tab_k <- data.frame(k = k_grid, RMSE = round(rmse_k, 4))
k_opt <- k_grid[which.min(rmse_k)]
kable(tab_k, caption = "RMSE NN untuk berbagai nilai k (10-fold CV)")
RMSE NN untuk berbagai nilai k (10-fold CV)
k RMSE
2 0.4825
3 0.4443
4 0.4818
5 0.4906
6 0.4860
7 0.4917
8 0.4972
9 0.5140
10 0.5182
11 0.5364
12 0.5422
13 0.5481
14 0.5572
15 0.5672 
cat("k optimum:", k_opt, "\n")
## k optimum: 3
ggplot(tab_k, aes(k, RMSE)) +
  geom_line() + geom_point() +
  geom_point(data = tab_k[tab_k$k == k_opt, ], color = "red", size = 3) +
  labs(title = "Pemilihan jumlah tetangga (k) pada NN", x = "k", y = "RMSE CV") +
  theme_minimal()

pred_nn <- idw(log_zinc ~ 1, meuse_sf, grid_sf,
               nmax = k_opt, idp = 0, debug.level = 0)

peta_pred(pred_nn$var1.pred, paste0("Nearest Neighbors (k = ", k_opt, ")"))

Pendekatan Ensembel

Pendekatan ensembel menggabungkan prediksi dari beberapa metode sekaligus dengan harapan kesalahan masing-masing metode saling mengimbangi. Ensembel pada tugas ini dibentuk dari tiga metode deterministik: Closest Observation, IDW, dan NN. Dua skema penggabungan digunakan:

  • Ensembel rata-rata (bobot sama): \(\hat{z}_{ens} = \frac{1}{3}(\hat{z}_{Closest} + \hat{z}_{IDW} + \hat{z}_{NN})\)
  • Ensembel terboboti: bobot proporsional terhadap \(1/\text{RMSE}\) tiap metode (dihitung pada bagian Evaluasi Model).
pred_ens <- (pred_closest$var1.pred + pred_idw$var1.pred + pred_nn$var1.pred) / 3

peta_pred(pred_ens, "Ensembel Rata-rata (Closest + IDW + NN)")

Berikut perbandingan peta keempat metode deterministik:

peta_det <- data.frame(
  X = coords[, 1], Y = coords[, 2],
  `Closest Observation` = pred_closest$var1.pred,
  IDW = pred_idw$var1.pred,
  `Nearest Neighbors` = pred_nn$var1.pred,
  `Ensembel (rata-rata)` = pred_ens,
  check.names = FALSE
) |>
  pivot_longer(-c(X, Y), names_to = "Metode", values_to = "log_zinc")

peta_det$Metode <- factor(peta_det$Metode,
  levels = c("Closest Observation", "IDW", "Nearest Neighbors", "Ensembel (rata-rata)"))

ggplot(peta_det, aes(X, Y, fill = log_zinc)) +
  geom_raster() +
  scale_fill_viridis_c(option = "plasma", name = "log(zinc)") +
  coord_equal() +
  facet_wrap(~ Metode, ncol = 2) +
  labs(title = "Peta Interpolasi Metode Deterministik") +
  theme_minimal() +
  theme(axis.text = element_text(size = 6))

Evaluasi Model

Evaluasi dilakukan dengan 10-fold cross validation: pada setiap iterasi model dibangun dari 9 fold (training) dan diuji pada 1 fold (testing). Hasil prediksi out-of-fold dikumpulkan lalu dihitung RMSE-nya.

cv_det <- function(data, folds, k_nn, idp) {
  res <- data.frame(obs = data$log_zinc, fold = folds,
                    closest = NA_real_, idw = NA_real_, nn = NA_real_)
  for (i in sort(unique(folds))) {
    tr  <- data[folds != i, ]
    te  <- data[folds == i, ]
    idx <- folds == i
    res$closest[idx] <- idw(log_zinc ~ 1, tr, te, nmax = 1,    idp = 0,   debug.level = 0)$var1.pred
    res$idw[idx]     <- idw(log_zinc ~ 1, tr, te,              idp = idp, debug.level = 0)$var1.pred
    res$nn[idx]      <- idw(log_zinc ~ 1, tr, te, nmax = k_nn, idp = 0,   debug.level = 0)$var1.pred
  }
  res
}

cv_res <- cv_det(meuse_sf, folds, k_nn = k_opt, idp = idp_opt)

# Ensembel rata-rata (Closest + IDW + NN)
cv_res$ens_mean <- (cv_res$closest + cv_res$idw + cv_res$nn) / 3
head(cv_res)

Bobot ensembel terboboti ditentukan dari kebalikan RMSE ketiga metode penyusun (semakin kecil RMSE, semakin besar bobotnya):

rmse_base <- c(Closest = rmse(cv_res$obs, cv_res$closest),
               IDW     = rmse(cv_res$obs, cv_res$idw),
               NN      = rmse(cv_res$obs, cv_res$nn))
w <- (1 / rmse_base) / sum(1 / rmse_base)
kable(data.frame(Metode = names(w), Bobot = round(w, 4)),
      caption = "Bobot Ensembel Terboboti (1/RMSE)")
Bobot Ensembel Terboboti (1/RMSE)
Metode Bobot
Closest Closest 0.2891
IDW IDW 0.3489
NN NN 0.3620 
cv_res$ens_wt <- w[["Closest"]] * cv_res$closest +
                 w[["IDW"]]     * cv_res$idw +
                 w[["NN"]]      * cv_res$nn

Catatan: bobot dihitung dari prediksi out-of-fold yang sama dengan yang dievaluasi, sehingga RMSE ensembel terboboti sedikit optimistik. Ensembel rata-rata bebas dari masalah ini.

ringkas <- function(obs, pred) {
  c(RMSE = rmse(obs, pred), MAE = mean(abs(obs - pred)), R2_CV = cor(obs, pred)^2)
}

hasil_det <- rbind(
  `Closest Observation`  = ringkas(cv_res$obs, cv_res$closest),
  `Inverse Distance Weight` = ringkas(cv_res$obs, cv_res$idw),
  `Nearest Neighbors`    = ringkas(cv_res$obs, cv_res$nn),
  `Ensembel (rata-rata)` = ringkas(cv_res$obs, cv_res$ens_mean),
  `Ensembel (terboboti)` = ringkas(cv_res$obs, cv_res$ens_wt)
) |> round(4)

kable(hasil_det, caption = "Evaluasi Metode Deterministik (10-fold CV, skala log(zinc))")
Evaluasi Metode Deterministik (10-fold CV, skala log(zinc))
RMSE MAE R2_CV
Closest Observation 0.5563 0.3890 0.4831
Inverse Distance Weight 0.4610 0.3549 0.5919
Nearest Neighbors 0.4443 0.3398 0.6207
Ensembel (rata-rata) 0.4546 0.3415 0.6010
Ensembel (terboboti) 0.4509 0.3404 0.6074 
cv_long_det <- cv_res |>
  select(obs, `Closest Observation` = closest, IDW = idw, `Nearest Neighbors` = nn,
         `Ensembel (rata-rata)` = ens_mean, `Ensembel (terboboti)` = ens_wt) |>
  pivot_longer(-obs, names_to = "Metode", values_to = "pred")

ggplot(cv_long_det, aes(obs, pred)) +
  geom_point(alpha = 0.6, color = "steelblue") +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
  facet_wrap(~ Metode, ncol = 3) +
  labs(title = "Observasi vs Prediksi - Metode Deterministik",
       x = "Observasi log(zinc)", y = "Prediksi log(zinc)") +
  theme_minimal()

Metode Stokastik

Metode stokastik memodelkan struktur autokorelasi spasial melalui semivariogram sehingga bobot pengamatan ditentukan oleh jarak sekaligus pola ketergantungan spasial data, dan prediksi dilengkapi dengan ukuran ketidakpastian (varians kriging). Metode yang digunakan adalah Ordinary Kriging.

Semivariogram

Model semivariogram (Spherical, Exponential, Gaussian) dicocokkan, kemudian dipilih model dengan sum of squares error (SSErr) terkecil. Fungsi yang sama dipakai ulang di dalam cross validation sehingga semivariogram diestimasi ulang pada setiap training set (tidak ada kebocoran data dari testing set).

fit_vgm <- function(d, models = c("Sph", "Exp", "Gau")) {
  v <- variogram(log_zinc ~ 1, d)
  best <- NULL
  best_sse <- Inf
  for (m in models) {
    f <- tryCatch(
      suppressWarnings(
        fit.variogram(v, vgm(psill = max(v$gamma) * 0.9, model = m,
                             range = max(v$dist) / 3, nugget = min(v$gamma)))
      ),
      error = function(e) NULL
    )
    if (!is.null(f) && !isTRUE(attr(f, "singular"))) {
      sse <- attr(f, "SSErr")
      if (sse < best_sse) { best <- f; best_sse <- sse }
    }
  }
  # cadangan jika semua model gagal
  if (is.null(best)) best <- vgm(psill = 0.6, model = "Sph", range = 900, nugget = 0.05)
  best
}

v_emp <- variogram(log_zinc ~ 1, meuse_sf)
v_fit <- fit_vgm(meuse_sf)
v_fit
plot(v_emp, v_fit, pch = 19, main = "Semivariogram log(zinc)")

Prediksi Kriging

pred_ok <- krige(log_zinc ~ 1, meuse_sf, grid_sf, model = v_fit, debug.level = 0)

peta_pred(pred_ok$var1.pred, "Ordinary Kriging")

ggplot(data.frame(X = coords[, 1], Y = coords[, 2], var = pred_ok$var1.var),
       aes(X, Y, fill = var)) +
  geom_raster() +
  scale_fill_viridis_c(option = "magma", name = "Varians") +
  coord_equal() +
  labs(title = "Varians Prediksi Ordinary Kriging") +
  theme_minimal()

Evaluasi Model Kriging dan Perbandingan Seluruh Metode

Kriging dievaluasi dengan fold yang sama seperti metode deterministik.

cv_res$ok <- NA_real_
for (i in 1:K) {
  tr  <- meuse_sf[folds != i, ]
  te  <- meuse_sf[folds == i, ]
  vf  <- fit_vgm(tr)
  cv_res$ok[folds == i] <- krige(log_zinc ~ 1, tr, te, model = vf, debug.level = 0)$var1.pred
}

Tabel perbandingan RMSE seluruh metode (diurutkan dari RMSE terkecil):

hasil_all <- rbind(
  `Closest Observation`     = ringkas(cv_res$obs, cv_res$closest),
  `Inverse Distance Weight` = ringkas(cv_res$obs, cv_res$idw),
  `Nearest Neighbors`       = ringkas(cv_res$obs, cv_res$nn),
  `Ensembel (rata-rata)`    = ringkas(cv_res$obs, cv_res$ens_mean),
  `Ensembel (terboboti)`    = ringkas(cv_res$obs, cv_res$ens_wt),
  `Kriging (Ordinary)`      = ringkas(cv_res$obs, cv_res$ok)
) |>
  as.data.frame() |>
  tibble::rownames_to_column("Metode") |>
  arrange(RMSE)

hasil_all[, -1] <- round(hasil_all[, -1], 4)
kable(hasil_all, caption = "Perbandingan Kinerja Seluruh Metode Interpolasi (10-fold CV, skala log(zinc))")
Perbandingan Kinerja Seluruh Metode Interpolasi (10-fold CV, skala log(zinc))
Metode RMSE MAE R2_CV
Kriging (Ordinary) 0.4046 0.3095 0.6851
Nearest Neighbors 0.4443 0.3398 0.6207
Ensembel (terboboti) 0.4509 0.3404 0.6074
Ensembel (rata-rata) 0.4546 0.3415 0.6010
Inverse Distance Weight 0.4610 0.3549 0.5919
Closest Observation 0.5563 0.3890 0.4831 
ggplot(hasil_all, aes(x = reorder(Metode, -RMSE), y = RMSE, fill = Metode)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = round(RMSE, 3)), hjust = -0.1) +
  coord_flip() +
  expand_limits(y = max(hasil_all$RMSE) * 1.12) +
  labs(title = "RMSE Cross Validation per Metode", x = NULL, y = "RMSE") +
  theme_minimal()

ggplot(data.frame(obs = cv_res$obs, pred = cv_res$ok), aes(obs, pred)) +
  geom_point(alpha = 0.6, color = "steelblue") +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Observasi vs Prediksi - Ordinary Kriging",
       x = "Observasi log(zinc)", y = "Prediksi log(zinc)") +
  theme_minimal()

Kesimpulan

  1. Berdasarkan 10-fold cross validation pada data meuse (peubah log(zinc)), metode dengan RMSE terkecil adalah Kriging (Ordinary) (RMSE = 0.4046), sedangkan metode dengan RMSE terbesar adalah Closest Observation (RMSE = 0.5563).
  2. Closest Observation hanya memakai satu tetangga sehingga sangat sensitif terhadap pengamatan tunggal dan menghasilkan peta berblok-blok.
  3. IDW (idp optimum = 3.5) dan NN (k optimum = 3) memanfaatkan lebih banyak tetangga sehingga permukaan lebih halus dan galat lebih kecil.
  4. Pendekatan Ensembel menggabungkan Closest Observation, IDW, dan NN; manfaat penggabungan dapat dilihat dengan membandingkan RMSE-nya terhadap metode penyusunnya pada tabel evaluasi.
  5. Kriging memanfaatkan struktur autokorelasi spasial lewat semivariogram (model terpilih: Sph, range = 897 m) dan memberikan ukuran ketidakpastian berupa varians prediksi, yang tidak dimiliki metode deterministik.

Referensi

  • Rusesta, R. R. Praktikum 10 : Interpolasi Titik. RPubs. https://rpubs.com/rositariarusesta/interpolasi_titik
  • Pebesma, E. J. (2004). Multivariable geostatistics in S: the gstat package. Computers & Geosciences, 30, 683-691.
  • Bivand, R. S., Pebesma, E., & Gomez-Rubio, V. (2013). Applied Spatial Data Analysis with R. Springer.