Konsep Dasar

Interpolasi spasial mengasumsikan bahwa nilai variabel yang diukur (misalnya, suhu, ketinggian, kelembaban tanah, dll.) di suatu lokasi akan serupa dengan nilai di lokasi-lokasi tetangga. Dengan kata lain, ada korelasi spasial, di mana lokasi yang lebih dekat cenderung memiliki nilai yang lebih mirip dibandingkan lokasi yang lebih jauh.

Data yang digunakan

Dalam bagian ini, seluruh metode interpolasi menggunakan data meuse dari paket sp (Pebesma & Bivand), yang berisi 155 titik pengamatan konsentrasi logam berat pada tanah lapisan atas (topsoil) di dataran banjir Sungai Meuse, Belanda. Variabel yang diinterpolasi adalah konsentrasi seng (zinc, ppm). Karena sebaran zinc menjulur ke kanan, digunakan transformasi log(zinc) pada seluruh metode, sama seperti pada bagian Kriging, sehingga RMSE antarmetode dapat dibandingkan secara adil. Lokasi prediksi adalah meuse.grid (grid 40 m x 40 m, 3.103 titik).

Langkah-langkah:

1. Memuat data meuse dan meuse.grid, lalu mengubahnya menjadi objek sf (CRS 28992, RD New).
2. Membuat variabel vble = log(zinc).

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

data(meuse)
data(meuse.grid)

# Data pengamatan (lokasi sampel)
d <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)
d$vble <- log(d$zinc)

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

# Batas wilayah studi (gabungan sel grid 40 m), dipakai untuk plot Voronoi
map <- st_union(st_buffer(coop, dist = 20, endCapStyle = "SQUARE"))

ggplot() +
  geom_sf(data = map, fill = "grey90", color = NA) +
  geom_sf(data = d, aes(color = vble), size = 1.8) +
  scale_color_viridis_c(name = "log(zinc)") + theme_bw()

Lokasi prediksi (coop) adalah titik pusat sel meuse.grid:

ggplot() + geom_sf(data = coop, size = 0.2) + theme_bw()

RMSE <- function(observed, predicted) {
  sqrt(mean((observed - predicted)^2))
}

# Closest observation: nilai titik terdekat (setara dengan sel Voronoi)
pred_closest <- function(train, newdata) {
  idx <- st_nearest_feature(newdata, train)
  train$vble[idx]
}

# IDW dengan semua titik latih, idp = 1
pred_idw <- function(train, newdata, idp = 1) {
  g <- gstat(formula = vble ~ 1, locations = train,
             nmax = nrow(train), set = list(idp = idp))
  predict(g, newdata, debug.level = 0)$var1.pred
}

# Nearest neighbors: rata-rata k tetangga terdekat (idp = 0 -> bobot sama)
pred_nn <- function(train, newdata, k = 5) {
  g <- gstat(formula = vble ~ 1, locations = train,
             nmax = k, set = list(idp = 0))
  predict(g, newdata, debug.level = 0)$var1.pred
}

# Ordinary kriging: variogram di-fit ulang pada data latih
pred_ok <- function(train, newdata) {
  vg <- variogram(vble ~ 1, train)
  fv <- suppressWarnings(
    fit.variogram(vg, vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1))
  )
  g <- gstat(formula = vble ~ 1, data = train, model = fv)
  predict(g, newdata, debug.level = 0)$var1.pred
}

Metode Deteministik

# Membuat diagram Voronoi dari lokasi pengamatan
vgeom <- st_collection_extract(st_voronoi(st_union(st_geometry(d))), "POLYGON")
vor <- st_sf(geometry = vgeom)
vor <- st_join(vor, d["vble"])
vor <- suppressWarnings(st_intersection(vor, map))

ggplot() +
  geom_sf(data = vor, aes(fill = vble), alpha = 0.7, color = "white", linewidth = 0.1) +
  geom_sf(data = d, size = 0.5) +
  scale_fill_viridis_c(name = "log(zinc)") + theme_bw()

p1 <- pred_closest(d, coop)

# Menerapkan IDW menggunakan semua lokasi pengamatan, lalu prediksi ke grid
p2 <- pred_idw(d, coop, idp = 1)

p3 <- pred_nn(d, coop, k = 5)

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

buat <- function(p, nama) {
  x <- coop
  x$pred <- p
  x$metode <- nama
  x[, c("metode", "pred")]
}

hasil <- rbind(
  buat(p1, "Closest observation"),
  buat(p2, "IDW"),
  buat(p3, "Nearest neighbors"),
  buat(p4, "Ensemble (bobot sama)")
)
hasil$metode <- factor(hasil$metode, levels = unique(hasil$metode))

ggplot(hasil) +
  geom_sf(aes(color = pred), size = 0.5) +
  facet_wrap(~ metode, ncol = 2) +
  scale_color_viridis_c(name = "log(zinc)") + theme_bw()

Metode Stokastik

# Memuat paket yang dibutuhkan
library(sp)        
library(gstat)     
library(sf)        
library(mapview)   

# Memuat dataset meuse dan meuse.grid yang sudah tersedia di R
data(meuse)        
data(meuse.grid)   

# Mengonversi dataset meuse menjadi objek 'sf' (simple feature) dengan koordinat 'x' dan 'y', dan CRS 28992 (proyeksi UTM)
meuse <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)
meuse.grid <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)

# Menampilkan peta interaktif 
mapview(meuse, zcol = "zinc", map.types = "CartoDB.Voyager")

mapview(meuse.grid, map.types = "CartoDB.Voyager")

# Variogram dengan cloud plot
vc <- variogram(log(zinc) ~ 1, meuse, cloud = TRUE)
plot(vc)

# Variogram dengan plot standar
v <- variogram(log(zinc) ~ 1, data = meuse)
plot(v)

vgm()

##    short                                      long
## 1    Nug                              Nug (nugget)
## 2    Exp                         Exp (exponential)
## 3    Sph                           Sph (spherical)
## 4    Gau                            Gau (gaussian)
## 5    Exc        Exclass (Exponential class/stable)
## 6    Mat                              Mat (Matern)
## 7    Ste Mat (Matern, M. Stein's parameterization)
## 8    Cir                            Cir (circular)
## 9    Lin                              Lin (linear)
## 10   Bes                              Bes (bessel)
## 11   Pen                      Pen (pentaspherical)
## 12   Per                            Per (periodic)
## 13   Wav                                Wav (wave)
## 14   Hol                                Hol (hole)
## 15   Log                         Log (logarithmic)
## 16   Pow                               Pow (power)
## 17   Spl                              Spl (spline)
## 18   Leg                            Leg (Legendre)
## 19   Err                   Err (Measurement error)
## 20   Int                           Int (Intercept)

# Menampilkan daftar model variogram
show.vgms(par.strip.text = list(cex = 0.75))

# 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 (vinitial) dan batas cutoff 1000
plot(v, vinitial, cutoff = 1000, cex = 1.5)

# Menyesuaikan variogram dengan model "Sph" menggunakan data v dan parameter
fv<-fit.variogram(object=v,model= 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(v,fv, cex= 1.5)

library(automap)
# Autofit variogram untuk data log-zinc
afvg <- autofitVariogram(log(zinc) ~ 1, meuse)

# Plot hasil
plot(afvg)

library(ggplot2)
library(viridis)
# Membuat objek gstat dengan model variogram hasil fitting (fv) untuk data zinc
k <- gstat(formula = log(zinc) ~ 1, data = meuse, model = fv)

# Melakukan prediksi variabel pada grid menggunakan model gstat yang sudah dibuat
kpred <- predict(k, meuse.grid)

## [using ordinary kriging]

# Visualisasi prediksi log(zinc) dengan warna berdasarkan var1.pred
ggplot() + geom_sf(data = kpred, aes(color = var1.pred)) +
  geom_sf(data = meuse) +
  scale_color_viridis(name = "log(zinc)") + theme_bw()

# Visualisasi prediksi varians dengan warna berdasarkan var1.var
ggplot() + geom_sf(data = kpred, aes(color = var1.var)) +
  geom_sf(data = meuse) +
  scale_color_viridis(name = "variance") + theme_bw()

Evaluasi Model (Cross Validation)

Untuk menilai kinerja setiap metode, dilakukan K-fold cross-validation (K = 5) dan dihitung RMSE (Root Mean Square Error). Cross-validation digunakan untuk memastikan bahwa model tidak hanya cocok dengan data latih tetapi juga dapat memprediksi data baru dengan baik. Metode yang dibandingkan:

1. Closest observation
2. IDW (idp = 1)
3. Nearest neighbors (k = 5)
4. Ensembel bobot sama (rata-rata Closest, IDW, NN)
5. Ensembel bobot 1/RMSE (Closest, IDW, NN)
6. Ordinary Kriging (variogram Sph di-fit ulang pada setiap data latih)

Langkah-langkah:

1. Membagi data menjadi K bagian dengan dismo::kfold().
2. Pada setiap fold, seluruh model (termasuk variogram Kriging) dilatih hanya pada data latih, lalu memprediksi data uji.
3. Bobot ensembel 1/RMSE dihitung dari cross validation di dalam data latih (inner CV), bukan dari data uji, supaya tidak terjadi kebocoran data (data leakage).
4. RMSE dihitung pada skala log(zinc) (skala model) dan pada skala asli (ppm, setelah exp()).

set.seed(123)
K <- 5

# Bobot ensembel = 1/RMSE hasil CV internal pada data latih
bobot_invrmse <- function(train, k_inner = 3) {
  fi <- kfold(nrow(train), k = k_inner)
  e <- matrix(NA, k_inner, 3)
  for (j in 1:k_inner) {
    tr <- train[fi != j, ]
    te <- train[fi == j, ]
    e[j, ] <- c(RMSE(te$vble, pred_closest(tr, te)),
                RMSE(te$vble, pred_idw(tr, te)),
                RMSE(te$vble, pred_nn(tr, te)))
  }
  w <- 1 / colMeans(e)
  w / sum(w)
}

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

metode <- c("Closest observation", "IDW", "Nearest neighbors",
            "Ensemble (bobot sama)", "Ensemble (bobot 1/RMSE)", "Kriging (OK)")
rmse_log  <- matrix(NA, K, length(metode),
                    dimnames = list(paste("Fold", 1:K), metode))
rmse_zinc <- rmse_log

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

  q1 <- pred_closest(train, test)
  q2 <- pred_idw(train, test)
  q3 <- pred_nn(train, test)
  q4 <- (q1 + q2 + q3) / 3

  w  <- bobot_invrmse(train)
  q5 <- w[1] * q1 + w[2] * q2 + w[3] * q3

  q6 <- pred_ok(train, test)

  P <- cbind(q1, q2, q3, q4, q5, q6)
  rmse_log[f, ]  <- apply(P, 2, function(p) RMSE(test$vble, p))
  rmse_zinc[f, ] <- apply(P, 2, function(p) RMSE(exp(test$vble), exp(p)))
}

RMSE pada skala log(zinc) untuk setiap fold dan rata-ratanya:

tab_log <- rbind(rmse_log, "Rata-rata" = colMeans(rmse_log))
kable(round(tab_log, 4), caption = "RMSE 5-fold CV, skala log(zinc)")

RMSE pada skala asli (ppm):

tab_zinc <- rbind(rmse_zinc, "Rata-rata" = colMeans(rmse_zinc))
kable(round(tab_zinc, 2), caption = "RMSE 5-fold CV, skala zinc (ppm, back-transform exp)")

rata <- data.frame(metode = factor(metode, levels = metode),
                   rmse = colMeans(rmse_log))

ggplot(rata, aes(x = metode, y = rmse, fill = metode)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = round(rmse, 3)), vjust = -0.4) +
  scale_fill_viridis_d() +
  labs(x = NULL, y = "Rata-rata RMSE (log(zinc))") +
  theme_bw() + theme(axis.text.x = element_text(angle = 25, hjust = 1))

Kesimpulan

Berdasarkan rata-rata RMSE 5-fold cross validation, metode dengan galat terkecil adalah Kriging (OK) (RMSE = 0.4312 pada skala log(zinc)), sedangkan metode dengan galat terbesar adalah IDW (RMSE = 0.6358).

Catatan Tambahan

1. Data. Seluruh metode (Closest observation, IDW, NN, Ensembel, Kriging) kini memakai data meuse sesuai ketentuan tugas. Data properties (Athena) tidak lagi dipakai.
2. Variabel respons sama. Semua metode memakai log(zinc). Pada versi sebelumnya Kriging memakai log(zinc), sedangkan metode lain memakai variabel berbeda sehingga RMSE tidak sebanding.
3. Kriging masuk perbandingan RMSE. Variogram di-fit ulang pada data latih di setiap fold (bukan memakai variogram dari seluruh data), agar tidak ada kebocoran informasi dari data uji.
4. Ensembel tanpa kebocoran data. Pada versi sebelumnya bobot 1/RMSE dihitung dari RMSE pada data uji yang sama dengan yang dievaluasi sehingga RMSE ensembel terlalu optimistis. Di sini bobot dihitung dari inner CV pada data latih. Ensembel bobot sama juga disertakan sebagai pembanding.
5. Urutan prediksi Voronoi. Pada versi sebelumnya st_intersection(v, coop)$vble mengembalikan baris berurutan menurut poligon Voronoi, bukan menurut urutan coop/test, sehingga dapat tidak sejajar dengan p2, p3, dan nilai observasi saat dihitung RMSE atau digabung dalam ensembel. Di sini st_nearest_feature() dipakai (hasil identik dengan sel Voronoi) sehingga urutan selalu sejajar. Diagram Voronoi tetap ditampilkan untuk visualisasi.
6. Grid prediksi. Untuk meuse, lokasi prediksi memakai meuse.grid langsung, sehingga tidak perlu membuat grid raster 100x100 sendiri. Peta memakai ggplot2 agar konsisten dengan bagian Kriging dan tidak bergantung pada versi tmap (argumen palette sudah usang di tmap v4).
7. Pesan R. Teks seperti [inverse distance weighted interpolation] dan [using ordinary kriging] adalah pesan keluaran R, bukan kode. Di sini disembunyikan dengan debug.level = 0 dan message = FALSE.
8. Nearest Neighbors. Metode ini memakai gstat dengan idp = 0 (bobot jarak dinolkan, semua tetangga berbobot sama) dan nmax = 5, sehingga hasilnya rata-rata 5 tetangga terdekat. Pesan “inverse distance weighted” pada keluaran R muncul karena fungsi yang dipakai sama dengan IDW.
9. Pembagian fold. set.seed(123) dipakai agar hasil dapat direproduksi. RMSE bergantung pada pembagian fold, sehingga angka dapat sedikit berbeda jika seed diubah.
10. Back-transform. RMSE skala ppm dihitung dengan exp() sederhana (tanpa koreksi bias log-normal), sehingga dipakai hanya untuk gambaran besaran galat. Perbandingan utama memakai RMSE skala log(zinc).
11. Penulisan. Judul “Metode Deteministik” pada dokumen asli sebaiknya “Deterministik” (dibiarkan agar teks asli tidak berubah).

Sumber Data dan Referensi

• Data meuse (paket sp): https://search.r-project.org/CRAN/refmans/sp/html/meuse.html
• Paket sp di CRAN: https://cran.r-project.org/package=sp
• Tutorial gstat dengan data meuse (Pebesma): https://search.r-project.org/CRAN/packages/gstat/vignettes/gstat.pdf
• Cressie, N. (1993). Statistics for Spatial Data. Wiley.