# ============================================================
# Tugas Interpolasi Spasial Data Meuse
# Metode: Closest Obs, NN, IDW, Kriging, Ensemble + Cross Validation
# ============================================================

# 1. Memuat Library yang Dibutuhkan
library(sp)
## Warning: package 'sp' was built under R version 4.5.3
library(sf)
## Warning: package 'sf' was built under R version 4.5.3
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
library(terra)
## Warning: package 'terra' was built under R version 4.5.3
## terra 1.9.34
library(tmap)
## Warning: package 'tmap' was built under R version 4.5.3
library(viridis)
## Warning: package 'viridis' was built under R version 4.5.3
## Loading required package: viridisLite
## Warning: package 'viridisLite' was built under R version 4.5.2
library(dismo)
## Warning: package 'dismo' was built under R version 4.5.3
## Loading required package: raster
## Warning: package 'raster' was built under R version 4.5.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:terra':
## 
##     intersect, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(gstat)
## Warning: package 'gstat' was built under R version 4.5.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.3
# ------------------------------------------------------------------------------
# 2. Persiapan Data (Meuse & Grid)
# ------------------------------------------------------------------------------
data(meuse)
data(meuse.grid)

# Mengonversi dataset meuse & meuse.grid menjadi objek 'sf' (CRS EPSG 28992)
d <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)
meuse_grid_sf <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)

# Membuat target variabel (log(zinc) atau zinc)
d$vble <- d$zinc 

# Membuat batas wilayah (outer boundary/map) dari convex hull data grid meuse
map <- st_union(st_geometry(meuse_grid_sf)) %>% st_convex_hull() %>% st_sf()

# Membuat grid raster untuk visualisasi
grid <- rast(ext(map), nrows = 100, ncols = 100, crs = crs(map))

# Mengambil titik pusat sel raster untuk lokasi prediksi
xy <- xyFromCell(grid, 1:ncell(grid))
coop <- st_as_sf(as.data.frame(xy), coords = c("x", "y"), crs = st_crs(map))
coop <- st_filter(coop, map) # Menyaring lokasi prediksi dalam batas map

# ------------------------------------------------------------------------------
# 3. Model Variogram untuk Kriging
# ------------------------------------------------------------------------------
v <- variogram(vble ~ 1, data = d)
fv <- fit.variogram(v, model = vgm(psill = 100000, model = "Sph", range = 900, nugget = 20000))

# Visualisasi Variogram
plot(v, fv, main = "Model Semivariogram (Zinc)")

# ------------------------------------------------------------------------------
# 4. Interpolasi Grid & Visualisasi Peta (5 Metode)
# ------------------------------------------------------------------------------

# A. Closest Observation (Voronoi Diagram)
v_vor <- voronoi(x = vect(d), bnd = map)
v_vor <- st_as_sf(v_vor)
resp1 <- st_intersection(v_vor, coop)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
resp1$pred <- resp1$vble
pred1 <- rasterize(resp1, grid, field = "pred", fun = "mean")
c1 <- tm_shape(pred1) + tm_raster(alpha = 0.6, palette = "viridis", title = "Closest")
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## map variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
## [v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
## map variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
# B. Nearest Neighbors (nmax = 5, idp = 0)
nn_mod <- gstat(formula = vble ~ 1, locations = d, nmax = 5, set = list(idp = 0))
resp2 <- predict(nn_mod, coop, debug.level = 0)
resp2$pred <- resp2$var1.pred
pred2 <- rasterize(resp2, grid, field = "pred", fun = "mean")
c2 <- tm_shape(pred2) + tm_raster(alpha = 0.6, palette = "viridis", title = "NN (k=5)")
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
## [v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
## map variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
# C. Inverse Distance Weighting (IDW)
idw_mod <- gstat(formula = vble ~ 1, locations = d, nmax = nrow(d), set = list(idp = 2))
resp3 <- predict(idw_mod, coop, debug.level = 0)
resp3$pred <- resp3$var1.pred
pred3 <- rasterize(resp3, grid, field = "pred", fun = "mean")
c3 <- tm_shape(pred3) + tm_raster(alpha = 0.6, palette = "viridis", title = "IDW")
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
## [v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
## map variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
# D. Ordinary Kriging
kr_mod <- gstat(formula = vble ~ 1, data = d, model = fv)
resp4 <- predict(kr_mod, coop, debug.level = 0)
resp4$pred <- resp4$var1.pred
pred4 <- rasterize(resp4, grid, field = "pred", fun = "mean")
c4 <- tm_shape(pred4) + tm_raster(alpha = 0.6, palette = "viridis", title = "Kriging")
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
## [v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
## map variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
# E. Ensemble (Rata-rata 3 Metode: Closest, IDW, Kriging)
p_ens <- (resp1$pred + resp3$pred + resp4$pred) / 3
resp5 <- data.frame(st_coordinates(coop), pred = p_ens)
resp5 <- st_as_sf(resp5, coords = c("X", "Y"), crs = st_crs(map))
pred5 <- rasterize(resp5, grid, field = "pred", fun = "mean")
c5 <- tm_shape(pred5) + tm_raster(alpha = 0.6, palette = "viridis", title = "Ensemble")
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
## [v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
## map variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
# Tampilkan Gabungan Peta Prediksi
tmap_arrange(c1, c2, c3, c4, c5, ncol = 3)

# ------------------------------------------------------------------------------
# 5. Evaluasi Model (K-Fold Cross Validation & RMSE)
# ------------------------------------------------------------------------------
set.seed(123)

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

# Pembagian data K-Fold (5 folds)
k_fold <- 5
kf <- kfold(nrow(d), k = k_fold)

# Wadah menyimpan nilai RMSE
rmse_closest <- rep(NA, k_fold)
rmse_nn      <- rep(NA, k_fold)
rmse_idw     <- rep(NA, k_fold)
rmse_kriging <- rep(NA, k_fold)
rmse_ens     <- rep(NA, k_fold)

for(k in 1:k_fold) {
  # Split data train & test
  test  <- d[kf == k, ]
  train <- d[kf != k, ]
  
  # 1. Closest Observation
  v_tr <- voronoi(x = vect(train), bnd = map)
  v_tr <- st_as_sf(v_tr)
  p1   <- st_intersection(v_tr, test)$vble
  rmse_closest[k] <- RMSE(test$vble, p1)
  
  # 2. Nearest Neighbors (k=5)
  nn_fit <- gstat(formula = vble ~ 1, locations = train, nmax = 5, set = list(idp = 0))
  p2     <- predict(nn_fit, test, debug.level = 0)$var1.pred
  rmse_nn[k] <- RMSE(test$vble, p2)
  
  # 3. IDW
  idw_fit <- gstat(formula = vble ~ 1, locations = train, nmax = nrow(train), set = list(idp = 2))
  p3      <- predict(idw_fit, test, debug.level = 0)$var1.pred
  rmse_idw[k] <- RMSE(test$vble, p3)
  
  # 4. Ordinary Kriging
  v_k <- variogram(vble ~ 1, data = train)
  fv_k <- fit.variogram(v_k, model = vgm(psill = 100000, model = "Sph", range = 900, nugget = 20000))
  kr_fit <- gstat(formula = vble ~ 1, data = train, model = fv_k)
  p4     <- predict(kr_fit, test, debug.level = 0)$var1.pred
  rmse_kriging[k] <- RMSE(test$vble, p4)
  
  # 5. Ensemble (Menggunakan bobot Inverse-RMSE dari 3 metode: Closest, IDW, Kriging)
  w <- 1 / c(rmse_closest[k], rmse_idw[k], rmse_kriging[k])
  weights <- w / sum(w)
  p5 <- p1 * weights[1] + p3 * weights[2] + p4 * weights[3]
  rmse_ens[k] <- RMSE(test$vble, p5)
}
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
# ------------------------------------------------------------------------------
# 6. Output Hasil Perbandingan RMSE
# ------------------------------------------------------------------------------

# Tabel Nilai RMSE untuk setiap Fold
rmse_folds <- data.frame(
  Fold = 1:5,
  Closest.Obs = rmse_closest,
  NN.k5 = rmse_nn,
  IDW = rmse_idw,
  Kriging = rmse_kriging,
  Ensemble = rmse_ens
)

# Rata-rata RMSE (Average RMSE)
rmse_summary <- data.frame(
  Closest.Obs = mean(rmse_closest),
  NN.k5       = mean(rmse_nn),
  IDW         = mean(rmse_idw),
  Kriging     = mean(rmse_kriging),
  Ensemble    = mean(rmse_ens)
)

cat("\n--- HASIL RMSE TIAP FOLD ---\n")
## 
## --- HASIL RMSE TIAP FOLD ---
print(rmse_folds)
##   Fold Closest.Obs    NN.k5      IDW  Kriging Ensemble
## 1    1    309.9009 265.3339 253.9304 214.2559 235.1864
## 2    2    385.8027 275.7933 249.1851 246.3962 261.7940
## 3    3    323.5046 364.2145 371.0162 300.6180 306.0539
## 4    4    237.6379 255.5252 259.6547 205.4405 217.3980
## 5    5    299.4034 294.7488 285.9440 251.4817 258.7884
cat("\n--- RATA-RATA RMSE (EVALUASI AKHIR) ---\n")
## 
## --- RATA-RATA RMSE (EVALUASI AKHIR) ---
print(rmse_summary)
##   Closest.Obs    NN.k5      IDW  Kriging Ensemble
## 1    311.2499 291.1231 283.9461 243.6385 255.8441

HASIL DAN PEMBAHASAN

  1. Pemodelan Semivariogram (Ordinary Kriging) Berdasarkan pemodelan semivariogram eksperimental untuk variabel Zinc, digunakan model fiting variogram tipe Spherical. Model ini memiliki nilai nugget sekitar 25000 yang mencerminkan variasi pada skala sangat kecil serta kesalahan pengukuran. Nilai sill berada pada kisaran 160000 yang menunjukkan batas variansi spasial maksimum, serta nilai range pada kisaran 800 hingga 900 meter yang menandakan jarak jangkauan ketergantungan spasial antar sampel.

  2. Visualisasi Peta Hasil Interpolasi Spasial Berdasarkan peta hasil prediksi kelima metode interpolasi, terlihat perbedaan karakteristik visual yang cukup jelas. Metode Closest Observation menghasilkan pola poligon Voronoi dengan batas wilayah yang sangat tegas karena nilai prediksi hanya ditentukan oleh satu titik sampel terdekat. Metode Nearest Neighbors (k=5) menghasilkan tampilan yang sedikit lebih halus daripada Closest Observation, namun batas antar wilayah masih terlihat kaku. Metode Inverse Distance Weighting (IDW) memperlihatkan efek bulls eye di sekitar titik-titik sampel bernilai tinggi. Sementara itu, metode Ordinary Kriging menghasilkan permukaan interpolasi yang paling mulus (smooth) karena memperhitungkan autokorelasi spasial melalui variogram. Metode Ensemble berhasil menyeimbangkan hasil estimasi visual dari metode-metode penyusunnya.

  3. Evaluasi Kinerja Model Menggunakan 5-Fold Cross Validation Pengujian kinerja model dilakukan menggunakan 5-Fold Cross Validation dengan indikator nilai Root Mean Square Error (RMSE). Dari hasil pengujian, metode Closest Observation memperoleh rata-rata RMSE sebesar 311.25. Metode Nearest Neighbors (k=5) menghasilkan rata-rata RMSE sebesar 291.12, sedangkan metode IDW memperoleh rata-rata RMSE sebesar 283.95. Penggabungan metode melalui Ensemble menghasilkan rata-rata RMSE sebesar 255.84. Adapun metode Ordinary Kriging memperoleh nilai rata-rata RMSE terkecil yaitu sebesar 243.64.

KESIMPULAN

Berdasarkan evaluasi 5-Fold Cross Validation, metode Ordinary Kriging merupakan metode interpolasi terbaik untuk data kadar Zinc pada dataset Meuse karena menghasilkan nilai rata-rata RMSE paling rendah, yaitu sebesar 243.64. Metode Ensemble berada di urutan kedua terbaik dengan nilai RMSE 255.84. Sebaliknya, metode Closest Observation memiliki tingkat kesalahan tertinggi dengan RMSE sebesar 311.25. Hal ini menunjukkan bahwa pemodelan ketergantungan spasial menggunakan fungsi variogram pada Ordinary Kriging mampu memberikan hasil estimasi lokasi yang paling akurat dibandingkan metode deterministik lainnya.