regspas 10

library(sp)
library(sf)

## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

library(gstat)
library(ggplot2)
library(tidyr)

data(meuse)
data(meuse.grid)

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

target_formula <- zinc ~ 1


print("Memproses Cross-Validation, mohon tunggu sebentar...")

## [1] "Memproses Cross-Validation, mohon tunggu sebentar..."

set.seed(33550336)

cv_closest <- krige.cv(target_formula, meuse_sf, nmax = 1)

cv_nn <- krige.cv(target_formula, meuse_sf, nmax = 5, set = list(idp = 0))

cv_idw <- krige.cv(target_formula, meuse_sf)

v <- variogram(target_formula, meuse_sf)
v_fit <- fit.variogram(v, vgm(psill = var(meuse_sf$zinc), model = "Sph", range = 900, nugget = min(v$gamma)))
cv_kriging <- krige.cv(target_formula, meuse_sf, model = v_fit)

ensemble_cv_pred <- (cv_idw$var1.pred + cv_kriging$var1.pred + cv_nn$var1.pred) / 3


calc_rmse <- function(observed, predicted) {
  sqrt(mean((observed - predicted)^2))
}
obs_aktual <- meuse_sf$zinc

tabel_rmse <- data.frame(
  Metode_Interpolasi = c("Closest Obs", "Nearest Neighbor", "IDW", "Kriging", "Ensembel"),
  RMSE = c(
    calc_rmse(obs_aktual, cv_closest$var1.pred),
    calc_rmse(obs_aktual, cv_nn$var1.pred),
    calc_rmse(obs_aktual, cv_idw$var1.pred),
    calc_rmse(obs_aktual, cv_kriging$var1.pred),
    calc_rmse(obs_aktual, ensemble_cv_pred)
  )
)

print("--- Hasil Perbandingan RMSE Cross Validation ---")

## [1] "--- Hasil Perbandingan RMSE Cross Validation ---"

print(tabel_rmse[order(tabel_rmse$RMSE), ])

##   Metode_Interpolasi     RMSE
## 4            Kriging 224.8003
## 5           Ensembel 247.1181
## 2   Nearest Neighbor 263.9874
## 3                IDW 278.2734
## 1        Closest Obs 309.9685

print("Memproses prediksi grid untuk plot...")

## [1] "Memproses prediksi grid untuk plot..."

pred_closest <- krige(target_formula, meuse_sf, meuse.grid_sf, nmax = 1)

## [inverse distance weighted interpolation]

pred_nn      <- krige(target_formula, meuse_sf, meuse.grid_sf, nmax = 5, set = list(idp = 0))

## [inverse distance weighted interpolation]

pred_idw     <- krige(target_formula, meuse_sf, meuse.grid_sf)

## [inverse distance weighted interpolation]

pred_kriging <- krige(target_formula, meuse_sf, meuse.grid_sf, model = v_fit)

## [using ordinary kriging]

pred_ensemble <- (pred_idw$var1.pred + pred_kriging$var1.pred + pred_nn$var1.pred) / 3


print("Men-generate plot...")

## [1] "Men-generate plot..."

grid_coords <- st_coordinates(meuse.grid_sf)

df_plot <- data.frame(
  x = grid_coords[, "X"],
  y = grid_coords[, "Y"],
  Closest_Obs = pred_closest$var1.pred,
  Nearest_Neighbor = pred_nn$var1.pred,
  IDW = pred_idw$var1.pred,
  Kriging = pred_kriging$var1.pred,
  Ensembel = pred_ensemble
)

metode_list <- c(
  "Closest_Obs" = "Closest Observation (1-NN)",
  "Nearest_Neighbor" = "Nearest Neighbor (5-NN Rata-rata)",
  "IDW" = "Inverse Distance Weighting (IDW)",
  "Kriging" = "Ordinary Kriging",
  "Ensembel" = "Ensembel (IDW + Kriging + 5-NN)"
)

for (kolom in names(metode_list)) {
  
  p <- ggplot(df_plot, aes(x = x, y = y, fill = .data[[kolom]])) +
    geom_tile() + 
    scale_fill_viridis_c(option = "magma", direction = -1) +
    theme_minimal() +
    coord_equal() + 
    labs(
      title = paste("Prediksi Spasial:", metode_list[[kolom]]),
      x = "Koordinat X",
      y = "Koordinat Y",
      fill = "Zinc (ppm)"
    ) +
    theme(
      plot.title = element_text(face = "bold", size = 14),
      legend.position = "right", 
      legend.title = element_text(face = "bold")
    )
  
  print(p)
}

Ordinary Kriging terbukti menjadi metode prediksi paling akurat dengan nilai Root Mean Square Error (RMSE) terendah sebesar 224.8003. Performa terbaik kedua dicapai oleh metode Ensembel (gabungan prediksi IDW, Kriging, dan 5-NN) yang memperoleh RMSE sebesar 247.1181, disusul oleh metode Nearest Neighbor dengan RMSE 263.9874 dan Inverse Distance Weighting (IDW) dengan RMSE 278.2734. Sebaliknya, tingkat eror tertinggi dihasilkan oleh metode Closest Observation (1-NN) dengan kinerja terendah yang mencapai RMSE 309.9685. Perbedaan performa kuantitatif ini juga sangat jelas terlihat pada hasil visualisasi peta spasialnya; metode Ordinary Kriging dan Ensembel mampu menampilkan pola gradasi transisi konsentrasi zinc yang jauh lebih mulus dan kontinu, sedangkan estimasi spasial pada metode Closest Observation dan Nearest Neighbor terlihat jauh lebih kaku dan terfragmentasi karena menyerupai blok-blok piksel poligon.