library(sp)        # data 'meuse'
library(sf)        # objek spasial (sf)
library(terra)     # raster & voronoi
library(gstat)     # IDW, NN, variogram, kriging
library(dismo)     # kfold
library(automap)   # autofit variogram
library(dplyr)
library(tidyr)
library(ggplot2)
library(viridis)
library(tmap)
library(patchwork) # menyusun beberapa peta

# Memuat data meuse
data(meuse)
data(meuse.riv)

# Konversi ke sf (koordinat x,y; CRS RD New = 28992)
d <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)

# Peubah yang diinterpolasi
d$vble <- d$zinc

# Convex hull dari titik pengamatan
map <- st_sf(geometry = st_convex_hull(st_union(d)))

# Konversi data sungai Meuse
riv_sf <- st_sf(geometry = st_sfc(st_linestring(meuse.riv), crs = 28992))

summary(d$zinc)

library(sf)
data(meuse.grid)   # reset ke data.frame asli
meuse.grid <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)

# Fungsi bantu: titik prediksi -> raster -> peta ggplot
to_rast <- function(sf_obj, col) {
  rasterize(vect(sf_obj), grid, field = col, fun = "mean")
}
peta_rast <- function(r, judul, nama = "Zinc", lim = NULL) {
  df <- as.data.frame(r, xy = TRUE); names(df)[3] <- "val"
  ggplot(df, aes(x, y, fill = val)) + geom_raster() + coord_equal() +
    scale_fill_viridis(name = nama, limits = lim) +
    theme_bw() + ggtitle(judul) +
    theme(axis.title = element_blank(), axis.text = element_blank(),
          axis.ticks = element_blank())
}
lim_zinc <- range(d$zinc)   # skala warna sama untuk semua peta

bb <- st_bbox(d)   # kotak pembatas titik data

ggplot() +
  geom_sf(data = map, fill = "grey95", alpha = 0.5) +
  geom_sf(data = riv_sf, colour = "#9fc5e8", linewidth = 1) +
  geom_sf(data = d, aes(colour = vble), size = 2) +
  coord_sf(xlim = c(bb["xmin"], bb["xmax"]),
           ylim = c(bb["ymin"], bb["ymax"]),
           expand = TRUE) +
  scale_colour_viridis(name = "Zinc (ppm)") +
  theme_bw() + ggtitle("Sebaran Kadar Zinc di Data meuse")

grid <- rast(map, nrows = 100, ncols = 100)
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)

qtm(coop)

# Membuat diagram Voronoi dari lokasi properti
v <- voronoi(x = vect(d), bnd = map)
plot(v)
# Menampilkan diagram Voronoi dengan harga properti
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")

# Tampilkan peta
c1

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

resp <- predict(res, coop)

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")
c3

library(gstat)

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

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

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

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

# Tampilkan peta
c2

#Closestobservation(Voronoi)
v<-voronoi(x=vect(d),bnd= map)
v<-st_as_sf(v)
p1<-st_intersection(v,coop)$vble

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

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

#Ensemble(equalweights)
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(resp,grid, field= "pred",fun= "mean")
c4<-tm_shape(pred)+tm_raster(alpha= 0.6,palette="viridis")
c4

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