resp <- coop
resp$pred <- d$vble[st_nearest_feature(coop, d)]
pred <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
c1 <- tm_shape(pred) + lay_raster(0.6) + lay_judul("Closest Observation") 
c1

# Menerapkan IDW menggunakan semua lokasi titik
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, debug.level = 0)

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

pred <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
c2 <- tm_shape(pred) + lay_raster(0.6) + lay_judul("IDW")
c2

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

resp <- predict(res, coop, debug.level = 0)

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

pred <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
c3 <- tm_shape(pred) + lay_raster(0.6) + lay_judul("Nearest Neighbors")
c3

# Closest observation (Voronoi)
p1 <- d$vble[st_nearest_feature(coop, d)]

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

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

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

resp <- data.frame(x = st_coordinates(coop)[, 1],
                   y = st_coordinates(coop)[, 2], pred = p4)
resp <- st_as_sf(resp, coords = c("x", "y"), crs = st_crs(map))
pred <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
c4 <- tm_shape(pred) + lay_raster(0.6) + lay_judul("Ensembel (CLosest Observation + IDW + NN)")  
c4

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

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) # Ensemble

for (k in 1:5) {
  # Split data into test and train
  test  <- d[kf == k, ]
  train <- d[kf != k, ]

  # Closest observation (nilai titik terdekat = prediksi Voronoi)
  p1 <- train$vble[st_nearest_feature(test, train)]
  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, debug.level = 0)$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, debug.level = 0)$var1.pred
  rmse3[k] <- RMSE(test$vble, p3)

  # Ensemble (weights are inverse RMSE so lower RMSE higher weight)
  w <- 1 / c(rmse1[k], rmse2[k], rmse3[k])
  weights <- w / sum(w)
  p4 <- p1 * weights[1] + p2 * weights[2] + p3 * weights[3]
  rmse4[k] <- RMSE(test$vble, p4)
}

# RMSE obtained for each of the 5 splits
data.frame(closest.obs = rmse1, IDW = rmse2,
           nearest.neigh = rmse3, ensemble = rmse4)

# Average RMSE over the 5 splits
data.frame(closest.obs = mean(rmse1), IDW = mean(rmse2),
           nearest.neigh = mean(rmse3), ensemble = mean(rmse4))

library(sp)
library(gstat)
library(sf)
library(mapview)

# Peta interaktif memakai koordinat lintang-bujur bawaan jura
d_ll <- st_as_sf(dat, coords = c("long", "lat"), crs = 4326)
d_ll$vble <- d$vble
mapview(d_ll, zcol = "vble", map.types = "CartoDB.Voyager")

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

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

vgm()

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

# Model variogram awal
MODEL_VGM <- "Sph"   # model yang di-fit; coba juga "Exp" lalu bandingkan RMSE
plateau  <- mean(tail(v$gamma, 6))
vinitial <- vgm(psill = max(plateau - min(v$gamma), 0.1 * plateau), model = MODEL_VGM,
                range = max(v$dist) / 3, nugget = min(v$gamma))
vinitial

plot(v, vinitial, cex = 1.5)

# Menyesuaikan (fit) variogram dengan model MODEL_VGM
fv <- fit.variogram(object = v, model = vinitial)
fv

plot(v, fv, cex = 1.5)

library(automap)

# Autofit variogram
afvg <- autofitVariogram(vble ~ 1, as(d, "Spatial"))

# Plot hasil
plot(afvg)

library(ggplot2)

library(viridis)

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

# Melakukan prediksi pada grid
kpred <- predict(k, coop)

# Visualisasi prediksi
ggplot() + geom_sf(data = kpred, aes(color = var1.pred)) +
  geom_sf(data = d) +
  scale_color_viridis(name = ifelse(LOG_TRANS, paste0("log(", VAR, ")"), VAR)) +
  theme_bw()

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

# Peta raster kriging, dibandingkan dengan metode lain
resp <- kpred
resp$pred <- resp$var1.pred
pred <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
c5 <- tm_shape(pred) + lay_raster(0.6) + lay_judul("Kriging (OK)")

tmap_arrange(c1, c2, c3, c4, c5, ncol = 3)

rmse5 <- rep(NA, 5) # Ordinary kriging

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

  vt  <- variogram(vble ~ 1, train)
  plt <- mean(tail(vt$gamma, 6))
  init <- function(m) vgm(psill = max(plt - min(vt$gamma), 0.1 * plt), model = m,
                        range = max(vt$dist) / 3, nugget = min(vt$gamma))
  fvt <- tryCatch(fit.variogram(vt, init(MODEL_VGM)), error = function(e) NULL)
  if (is.null(fvt) || isTRUE(attr(fvt, "singular"))) fvt <- fit.variogram(vt, init("Exp"))

  kr <- gstat(formula = vble ~ 1, data = train, model = fvt)
  p5 <- predict(kr, test, debug.level = 0)$var1.pred
  rmse5[k] <- RMSE(test$vble, p5)
}

# RMSE per fold, semua metode
data.frame(closest.obs = rmse1, IDW = rmse2, nearest.neigh = rmse3,
           ensemble = rmse4, kriging = rmse5)

# Rata-rata RMSE, semua metode
data.frame(closest.obs = mean(rmse1), IDW = mean(rmse2),
           nearest.neigh = mean(rmse3), ensemble = mean(rmse4),
           kriging = mean(rmse5))