Pendahuluan

Menerapkan interpolasi spasial pada data meuse (kadar logam berat di tanah dataran banjir Sungai Meuse, Belanda) dengan metode Closest Observation, Nearest Neighbors (NN), IDW, Ensembel, dan Kriging. Kinerja kelima metode dibandingkan dengan RMSE hasil K-fold Cross Validation.

Variabel yang diinterpolasi adalah log(zinc), sama seperti pada contoh kriging di praktikum. Semua RMSE berada pada skala log(zinc). Ensembel menggabungkan empat metode: Closest Observation, IDW, NN, dan Kriging.

Data yang Digunakan

# Memuat pustaka yang diperlukan
library(sp)
library(sf)
library(terra)
library(tmap)
library(viridis)
library(dismo)
library(dplyr)
library(gstat)
library(ggplot2)
library(knitr)

# Memuat data meuse (titik pengamatan) dan meuse.grid (area prediksi)
data(meuse)
data(meuse.grid)

# Mengubah menjadi objek sf (CRS 28992 = Amersfoort / RD New)
d <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)
grid_pts <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)

# Menambahkan kolom variabel yang diinterpolasi: log(zinc)
d$vble <- log(d$zinc)

# Membuat batas wilayah studi dari sel grid 40 m x 40 m
# (ditambah buffer kecil di sekitar titik agar semua titik sampel berada di dalam batas)
sel <- st_buffer(grid_pts, dist = 20, endCapStyle = "SQUARE")
map <- st_union(c(st_geometry(sel), st_geometry(st_buffer(d, 20)))) %>% st_sf()

# Menampilkan peta log(zinc)
tmap_mode("plot")
tm_shape(map) + 
  tm_polygons(alpha = 0.3) + 
  tm_shape(d) + 
  tm_dots("vble", palette = "viridis")

Berikut grid raster yang digunakan untuk prediksi. Titik pusat tiap sel menjadi lokasi prediksi.

# Membuat grid raster 100x100 sel untuk menutupi wilayah studi
grid <- rast(map, nrows = 100, ncols = 100)

# Mengambil koordinat titik pusat dari setiap sel dalam grid
xy <- xyFromCell(grid, 1:ncell(grid))

# Membuat objek sf untuk lokasi prediksi
coop <- st_as_sf(as.data.frame(xy), coords = c("x", "y"), crs = st_crs(map))

# Menyaring lokasi prediksi yang berada dalam batas wilayah studi
coop <- st_filter(coop, map)

# Menampilkan peta lokasi prediksi
qtm(coop)

Metode Deterministik

Closest Observation

# Membuat diagram Voronoi dari lokasi sampel
v <- terra::voronoi(x = vect(d), bnd = vect(map))
plot(v)
points(vect(d), cex = 0.5)

# Prediction
v <- st_as_sf(v)
tm_shape(v) + 
  tm_fill(col = "vble", palette = "viridis", alpha = 0.6)

resp <- st_intersection(v, coop)
resp$pred <- resp$vble
pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c1 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis", title = "Closest Obs")
c1

Inverse Distance Weighting (IDW)

# Menerapkan IDW menggunakan semua lokasi sampel
res <- gstat(formula = vble ~ 1, locations = d, nmax = nrow(d), set = list(idp = 1))

# Mendapatkan prediksi untuk lokasi yang belum teramati
resp <- predict(res, coop)

## [inverse distance weighted interpolation]

resp$x <- st_coordinates(resp)[,1]
resp$y <- st_coordinates(resp)[,2]
resp$pred <- resp$var1.pred

# Membuat peta prediksi dengan metode IDW
pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c2 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis", title = "IDW")
c2

Nearest Neighbors

res <- gstat(formula = vble ~ 1, locations = d, nmax = 5, set = list(idp = 0))
resp <- predict(res, coop)

## [inverse distance weighted interpolation]

resp$x <- st_coordinates(resp)[,1]
resp$y <- st_coordinates(resp)[,2]
resp$pred <- resp$var1.pred

pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c3 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis", title = "NN")
c3

Metode Stokastik: Kriging

Variogram

# Variogram dengan cloud plot
vc <- variogram(vble ~ 1, d, cloud = TRUE)
plot(vc)

# Variogram dengan plot standar
vg <- variogram(vble ~ 1, data = d)
plot(vg)

Fitting Model Variogram

# Menentukan model variogram awal dengan parameter psill, model, range, dan nugget
vinitial <- vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1)

# Memvisualisasikan variogram dengan model awal
plot(vg, vinitial, cutoff = 1000, cex = 1.5)

# Menyesuaikan model variogram dengan data
fv <- fit.variogram(object = vg, model = vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1))
fv

plot(vg, fv, cex = 1.5)

Prediksi Kriging

# Membuat objek gstat dengan model variogram hasil fitting
k <- gstat(formula = vble ~ 1, locations = d, model = fv)

# Prediksi pada seluruh lokasi grid
resp <- predict(k, coop)

## [using ordinary kriging]

resp$x <- st_coordinates(resp)[,1]
resp$y <- st_coordinates(resp)[,2]
resp$pred <- resp$var1.pred

pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c4 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis", title = "Kriging")
c4

# Visualisasi varians prediksi
ggplot() + 
  geom_sf(data = resp, aes(color = var1.var)) + 
  geom_sf(data = d) + 
  scale_color_viridis(name = "variance") + 
  theme_bw()

Pendekatan Ensembel

Ensembel menggabungkan prediksi dari empat metode: Closest Observation, IDW, Nearest Neighbors, dan Kriging, dengan bobot yang sama.

# Closest observation (Voronoi)
# (st_join menjaga urutan lokasi prediksi sama dengan urutan coop)
v <- terra::voronoi(x = vect(d), bnd = vect(map))
v <- st_as_sf(v)
v$pred <- v$vble
p1 <- st_join(coop, v["pred"])$pred

# IDW
gs <- gstat(formula = vble ~ 1, locations = d, nmax = nrow(d), set = list(idp = 1))
p2 <- predict(gs, coop)$var1.pred

## [inverse distance weighted interpolation]

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

## [inverse distance weighted interpolation]

# Kriging
kk <- gstat(formula = vble ~ 1, locations = d, model = fv)
p4 <- predict(kk, coop)$var1.pred

## [using ordinary kriging]

# Ensemble (equal weights)
weights <- c(1/4, 1/4, 1/4, 1/4)
p5 <- p1 * weights[1] + p2 * weights[2] + p3 * weights[3] + p4 * weights[4]

resp <- data.frame(x = st_coordinates(coop)[,1], y = st_coordinates(coop)[,2], pred = p5)
resp <- st_as_sf(resp, coords = c("x", "y"), crs = st_crs(map))
pred <- rasterize(resp, grid, field = "pred", fun = "mean")

c5 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis", title = "Ensembel")
c5

# Gabungkan semua plot dalam satu frame
tmap_arrange(c1, c2, c3, c4, c5, ncol = 3)

Evaluasi Model: Cross Validation

Dilakukan \(K\)-fold cross validation dengan \(K = 5\). Pada setiap lipatan, model (termasuk variogram kriging) dibangun ulang hanya dengan data latih, lalu RMSE dihitung pada data uji.

Ensembel dievaluasi dengan dua skema bobot:

Ensembel (bobot sama): bobot masing-masing metode \(1/4\).
Ensembel (bobot 1/RMSE): bobot proporsional dengan kebalikan RMSE. Bobot ini dihitung dari RMSE pada data uji yang sama, sehingga hasilnya cenderung sedikit optimistis.

set.seed(123)

# Fungsi untuk menghitung RMSE
RMSE <- function(observed, predicted) {
  sqrt(mean((observed - predicted)^2))
}

# Pembagian data untuk K-fold
kf <- kfold(nrow(d), k = 5)

# Vectors to store the RMSE values obtained with each method
rmse1 <- rep(NA, 5) # Closest observation
rmse2 <- rep(NA, 5) # IDW
rmse3 <- rep(NA, 5) # Nearest neighbors
rmse4 <- rep(NA, 5) # Kriging
rmse5 <- rep(NA, 5) # Ensemble (equal weights)
rmse6 <- rep(NA, 5) # Ensemble (weights = inverse RMSE)

for (k in 1:5) {
  # Split data in test and train
  test <- d[kf == k, ]
  train <- d[kf != k, ]
  
  # Closest observation
  v <- terra::voronoi(x = vect(train), bnd = vect(map))
  v <- st_as_sf(v)
  v$pred <- v$vble
  p1 <- st_join(test, v["pred"])$pred
  rmse1[k] <- RMSE(test$vble, p1)
  
  # IDW
  gs <- gstat(formula = vble ~ 1, locations = train, nmax = nrow(train), set = list(idp = 1))
  p2 <- predict(gs, test)$var1.pred
  rmse2[k] <- RMSE(test$vble, p2)
  
  # Nearest neighbors
  nn <- gstat(formula = vble ~ 1, locations = train, nmax = 5, set = list(idp = 0))
  p3 <- predict(nn, test)$var1.pred
  rmse3[k] <- RMSE(test$vble, p3)
  
  # Kriging (variogram di-fit ulang menggunakan data latih)
  vtr <- variogram(vble ~ 1, data = train)
  fvt <- suppressWarnings(
    fit.variogram(vtr, model = vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1)))
  kr <- gstat(formula = vble ~ 1, locations = train, model = fvt)
  p4 <- predict(kr, test)$var1.pred
  rmse4[k] <- RMSE(test$vble, p4)
  
  # Ensemble (equal weights)
  p5 <- (p1 + p2 + p3 + p4) / 4
  rmse5[k] <- RMSE(test$vble, p5)
  
  # Ensemble (weights are inverse RMSE so lower RMSE higher weight)
  w <- 1 / c(rmse1[k], rmse2[k], rmse3[k], rmse4[k])
  weights <- w / sum(w)
  p6 <- p1 * weights[1] + p2 * weights[2] + p3 * weights[3] + p4 * weights[4]
  rmse6[k] <- RMSE(test$vble, p6)
}

## [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]

# RMSE obtained for each of the 5 splits
rmse_fold <- data.frame(closest.obs = rmse1, IDW = rmse2, nearest.neigh = rmse3,
                        kriging = rmse4, ensemble.equal = rmse5, ensemble.invRMSE = rmse6)
kable(rmse_fold, digits = 4, caption = "RMSE tiap lipatan (skala log(zinc))")

RMSE tiap lipatan (skala log(zinc))
closest.obs	IDW	nearest.neigh	kriging	ensemble.equal	ensemble.invRMSE
0.6327	0.5496	0.4958	0.4281	0.4508	0.4458
0.5987	0.5828	0.4943	0.4183	0.4492	0.4439
0.6328	0.7843	0.6516	0.5055	0.5995	0.5864
0.4874	0.6333	0.5316	0.4198	0.4633	0.4542
0.5593	0.6293	0.4791	0.3840	0.4463	0.4347

# Average RMSE over the 5 splits
rmse_mean <- data.frame(closest.obs = mean(rmse1), IDW = mean(rmse2), nearest.neigh = mean(rmse3),
                        kriging = mean(rmse4), ensemble.equal = mean(rmse5),
                        ensemble.invRMSE = mean(rmse6))
kable(rmse_mean, digits = 4, caption = "Rata-rata RMSE 5-fold CV")

Rata-rata RMSE 5-fold CV
closest.obs	IDW	nearest.neigh	kriging	ensemble.equal	ensemble.invRMSE
0.5822	0.6358	0.5305	0.4312	0.4818	0.473

ringkas <- data.frame(
  Metode = c("Closest Obs", "IDW", "NN", "Kriging",
             "Ensembel (bobot sama)", "Ensembel (bobot 1/RMSE)"),
  `Rata-rata RMSE` = as.numeric(rmse_mean[1, ]),
  `Simpangan Baku` = apply(rmse_fold, 2, sd),
  check.names = FALSE
)
ringkas <- ringkas[order(ringkas$`Rata-rata RMSE`), ]
ringkas$Peringkat <- 1:nrow(ringkas)
rownames(ringkas) <- NULL
kable(ringkas, digits = 4, caption = "Ringkasan RMSE urut dari terbaik")

Ringkasan RMSE urut dari terbaik
Metode	Rata-rata RMSE	Simpangan Baku	Peringkat
Kriging	0.4312	0.0449	1
Ensembel (bobot 1/RMSE)	0.4730	0.0638	2
Ensembel (bobot sama)	0.4818	0.0661	3
NN	0.5305	0.0704	4
Closest Obs	0.5822	0.0610	5
IDW	0.6358	0.0899	6

ggplot(ringkas, aes(x = reorder(Metode, `Rata-rata RMSE`), y = `Rata-rata RMSE`, fill = Metode)) +
  geom_col(show.legend = FALSE) +
  geom_errorbar(aes(ymin = `Rata-rata RMSE` - `Simpangan Baku`,
                    ymax = `Rata-rata RMSE` + `Simpangan Baku`), width = 0.2) +
  scale_fill_viridis_d() +
  coord_flip() +
  theme_bw() +
  labs(x = NULL, y = "Rata-rata RMSE (log zinc)",
       title = "Perbandingan RMSE 5-fold CV (batang galat = simpangan baku antar lipatan)")

Interpretasi

Metode terbaik. Berdasarkan rata-rata RMSE 5-fold CV, metode dengan galat prediksi terkecil adalah Kriging (RMSE = 0.4312), sedangkan yang terbesar adalah IDW (RMSE = 0.6358). Selisihnya sekitar 32.2% dari RMSE terbesar.
Closest Observation (RMSE = 0.5822). Peta prediksinya berupa blok-blok poligon Voronoi dengan batas tegas, karena setiap lokasi hanya meminjam nilai satu titik terdekat. Metode ini tidak meredam variasi lokal atau pencilan, sehingga cenderung kurang akurat dibanding metode yang merata-ratakan banyak titik.
NN (RMSE = 0.5305) dan IDW (RMSE = 0.6358). Keduanya menghasilkan permukaan yang lebih halus. NN memberi bobot sama pada 5 tetangga terdekat, sedangkan IDW memberi bobot lebih besar pada titik yang lebih dekat. Kelemahannya, tidak ada ukuran ketidakpastian dan hasilnya bergantung pada pilihan parameter (\(k\) pada NN dan pangkat jarak pada IDW).
Kriging (RMSE = 0.4312). Kriging memakai struktur autokorelasi spasial dari variogram (model spherical), sehingga bobot tiap titik mengikuti pola ketergantungan spasial data, bukan hanya jarak. Kriging juga memberi varians prediksi, yang rendah di sekitar titik sampel dan membesar di area yang jauh dari sampel.
Ensembel (bobot sama: 0.4818, bobot 1/RMSE: 0.473). Penggabungan empat metode meredam galat masing-masing, sehingga prediksinya lebih stabil. Namun ensembel tidak selalu mengalahkan metode tunggal terbaik, terutama jika ada metode yang jauh lebih buruk, karena metode itu ikut menarik rata-rata. Bobot 1/RMSE memberi peran lebih besar pada metode yang lebih akurat, tetapi bobotnya dihitung dari data uji yang sama sehingga hasilnya sedikit optimistis.
Catatan. Hasil bergantung pada pembagian lipatan acak (set.seed(123)) dan data yang hanya 155 titik, sehingga selisih RMSE yang kecil antar metode (lihat simpangan baku antar lipatan) belum tentu bermakna. Pada data meuse, zinc berasosiasi kuat dengan jarak ke sungai, sehingga kriging dengan kovariat (universal kriging) berpotensi memberi galat lebih kecil, tetapi itu di luar cakupan tugas ini.

Kesimpulan

Dari metode yang dibandingkan lewat 5-fold cross validation, Kriging memberikan kinerja prediksi terbaik untuk log(zinc) pada data meuse, dengan rata-rata RMSE 0.4312.

Tugas Praktikum Interpolasi REGSPAS - Data Meuse

Lara Tria Novita

2026-11-09