Data

# Memuat paket yang dibutuhkan
library(spData)

## Warning: package 'spData' was built under R version 4.5.3

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

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.2

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

## Loading required package: sp

## Warning: package 'sp' 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(sp)        
library(gstat)     

## Warning: package 'gstat' was built under R version 4.5.3

library(sf)        
library(mapview)   

## Warning: package 'mapview' was built under R version 4.5.3

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

Data yang akan di interpolasi adalah data zinc dari data meuse

Closest Observation

Metode ini menggunakan diagram Voronoi untuk membagi ruang menjadi poligon yang masing-masing berisi satu titik penghasil. Setiap titik dalam poligon lebih dekat ke titik penghasil daripada ke titik lainnya. Prediksi harga properti dilakukan berdasarkan nilai titik terdekat.

# Membuat batas wilayah berdasarkan grid
boundary <- as.polygons(
  ext(vect(meuse.grid)),
  crs = crs(vect(meuse.grid))
)
# Membuat diagram Voronoi 
v <- voronoi(x = vect(meuse), bnd = boundary)
plot(v)
# Menampilkan diagram Voronoi
points(vect(meuse), cex = 0.5)

zinc <- meuse$zinc
# Prediction
v <- st_as_sf(v)
v$zinc <- zinc
tm_shape(v) + tm_fill(col = "zinc", palette = "viridis", alpha = 0.6)

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the map variable `fill` namely 'palette' (rename to 'values') to fill.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.

# Mengambil nilai zinc berdasarkan poligon Voronoi
resp <- st_join(
  meuse.grid,
  v["zinc"],
  join = st_within
)

# Melihat hasil prediksi
head(resp)

## Simple feature collection with 6 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 181100 ymin: 333660 xmax: 181220 ymax: 333740
## Projected CRS: Amersfoort / RD New
##   part.a part.b      dist soil ffreq zinc              geometry
## 1      1      0 0.0000000    1     1   NA POINT (181180 333740)
## 2      1      0 0.0000000    1     1  403 POINT (181140 333700)
## 3      1      0 0.0122243    1     1  403 POINT (181180 333700)
## 4      1      0 0.0434678    1     1  258 POINT (181220 333700)
## 5      1      0 0.0000000    1     1  403 POINT (181100 333660)
## 6      1      0 0.0122243    1     1  403 POINT (181140 333660)

# Membuat kolom prediksi
resp$pred <- resp$zinc

# Mengubah data sf menjadi SpatVector
resp_vect <- vect(resp)

# Membuat raster template
r <- rast(
  vect(meuse.grid),
  ncols = 50,
  nrows = 50
)

# Melakukan rasterisasi
pred <- rasterize(
  resp_vect,
  r,
  field = "pred"
)

# Menampilkan hasil prediksi
names(pred) <- "zinc"
c1 <- tm_shape(pred) +
  tm_raster(
    palette = "viridis",
    alpha = 0.6
  ) 

## 

## ── 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`.

c1

Invers Distance Weight

Metode IDW (Inverse Distance Weighting) menggunakan bobot berdasarkan jarak antara lokasi yang belum teramati dan lokasi-lokasi yang sudah teramati. Semakin dekat jarak, semakin besar bobot yang diberikan pada nilai lokasi yang teramati.

library(gstat)

# Menerapkan IDW menggunakan semua lokasi
res1 <- gstat(formula = zinc ~ 1, locations = meuse, nmax = nrow(meuse), set = list(idp = 1))

# Mendapatkan prediksi zinc untuk lokasi yang belum teramati
resp1 <- predict(res1, meuse.grid)

## [inverse distance weighted interpolation]

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

# Mengambil koordinat grid
coords1 <- st_coordinates(meuse.grid)

# Menghitung jarak antar titik grid
dx <- min(diff(sort(unique(coords1[,1]))))
dy <- min(diff(sort(unique(coords1[,2]))))

# Membuat raster template
grid1 <- rast(
  xmin = min(coords1[,1]) - dx/2,
  xmax = max(coords1[,1]) + dx/2,
  ymin = min(coords1[,2]) - dy/2,
  ymax = max(coords1[,2]) + dy/2,
  ncols = length(unique(coords1[,1])),
  nrows = length(unique(coords1[,2])),
  crs = st_crs(meuse.grid)$wkt
)

# Mengubah hasil prediksi menjadi SpatVector
resp1_vect <- vect(resp1)

# Membuat raster prediksi IDW
pred1 <- rasterize(
  resp1_vect,
  grid1,
  field = "pred",
  fun = "mean"
)

c2 <- tm_shape(pred1) +
  tm_raster(
    alpha = 0.6,
    palette = "viridis"
  )

## 

## ── 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`.

c2

Nearest Neighboor

Metode Nearest Neighbors menghitung rata-rata harga dari k titik terdekat. Setiap tetangga dianggap memiliki bobot yang sama.

res2 <- gstat(formula = zinc ~ 1, locations = meuse, nmax = 5, set = list(idp = 0))

resp2 <- predict(res2, meuse.grid)

## [inverse distance weighted interpolation]

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

# Mengambil koordinat grid
coords2 <- st_coordinates(meuse.grid)

# Menghitung jarak antar titik grid
dx <- min(diff(sort(unique(coords2[,1]))))
dy <- min(diff(sort(unique(coords2[,2]))))

# Membuat raster template
grid2 <- rast(
  xmin = min(coords2[,1]) - dx/2,
  xmax = max(coords2[,1]) + dx/2,
  ymin = min(coords2[,2]) - dy/2,
  ymax = max(coords2[,2]) + dy/2,
  ncols = length(unique(coords2[,1])),
  nrows = length(unique(coords2[,2])),
  crs = st_crs(meuse.grid)$wkt
)

# Mengubah hasil prediksi menjadi SpatVector
resp2_vect <- vect(resp2)

# Membuat raster prediksi NN
pred2 <- rasterize(
  resp2_vect,
  grid2,
  field = "pred",
  fun = "mean"
)

c3 <- tm_shape(pred2) +
  tm_raster(
    alpha = 0.6,
    palette = "viridis"
  )

## 

## ── 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`.

c3

Ensembel

Pendekatan Ensembel menggabungkan hasil prediksi dari berbagai metode interpolasi (seperti Closest Observation, IDW, dan Nearest Neighbors) dengan memberikan bobot pada setiap metode. Bobot ini bisa berdasarkan kinerja masing-masing metode.

#Closestobservation(Voronoi)
v<-voronoi(x=vect(meuse),bnd= boundary)
v<-st_as_sf(v)
p1<-st_intersection(v,meuse.grid)$zinc

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

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

## [inverse distance weighted interpolation]

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

## [inverse distance weighted interpolation]

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

resp3<-data.frame(x=st_coordinates(meuse.grid)[,1], y=st_coordinates(meuse.grid)[,2], pred=p4)

resp3<-st_as_sf(resp3,coords= c("x","y"),crs= st_crs(boundary))
pred<-rasterize(resp3,grid1, field= "pred",fun= "mean")
c4<-tm_shape(pred)+tm_raster(alpha= 0.6,palette="viridis")

## 

## ── 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`.

c4

Kriging

Kriging adalah metode interpolasi geostatistik yang memperhitungkan autokorelasi spasial untuk membuat prediksi. Ordinary Kriging (OK) adalah salah satu varian Kriging yang paling umum digunakan, di mana nilai yang akan diprediksi adalah variabel z yang tidak diketahui, yang dianggap sebagai variabel acak.

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

# Variogram dengan plot standar
v <- variogram (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))

## Warning in fit.variogram(object = v, model = vgm(psill = 0.5, model = "Sph", :
## No convergence after 200 iterations: try different initial values?

fv

##   model     psill   range
## 1   Nug 0.1041661   0.000
## 2   Sph 0.5220512 897.703

plot(v,fv, cex= 1.5)

library(automap)

## Warning: package 'automap' was built under R version 4.5.3

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

# Plot hasil
plot(afvg)

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.5.3

library(viridis)
# Membuat objek gstat dengan model variogram hasil fitting (fv) untuk data zinc
k <- gstat(formula = 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
C5 <- ggplot() + geom_sf(data = kpred, aes(color = var1.pred)) +
  geom_sf(data = meuse) +
  scale_color_viridis(name = "zinc") + theme_bw()

C5

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

C6

Perbandingan RMSE Model dengan Cross Validation

Untuk menilai kinerja setiap metode, kita melakukan K-fold cross-validation dan menghitung 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.

set.seed(123)

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

# Pembagian data untuk K-fold
kf <- kfold(nrow(meuse), 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
rmse5 <- rep(NA, 5) # Kriging

for(k in 1:5){
#Splitdataintest andtrain
test<-meuse[kf== k,]
train<-meuse[kf !=k,]


#Closestobservation
v<-voronoi(x=vect(train), bnd= boundary)
v<-st_as_sf(v)
p1<-st_intersection(v,test)$zinc
rmse1[k] <-RMSE(test$zinc,p1)

#IDW
gs<-gstat(formula=zinc ~1,locations= train,
nmax=nrow(train), set= list(idp=1))
p2 <-predict(gs,test)$var1.pred
rmse2[k] <-RMSE(test$zinc,p2)

#Nearestneighbors
nn<-gstat(formula=zinc~1,locations= train,
nmax=5,set= list(idp=0))
p3 <-predict(nn,test)$var1.pred
rmse3[k] <-RMSE(test$zinc,p3)

#Ensemble(weightsareinverseRMSEsolowerRMSEhigherweight)
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$zinc,p4)

# Kriging
# Menghitung semivariogram empiris
vg <- variogram(zinc ~ 1, train)
# Menghitung jarak maksimum
max_dist <- max(dist(st_coordinates(train)))
# Menentukan model semivariogram awal
initial_model <- vgm(
  psill = var(train$zinc) * 0.8,
  model = "Sph",
  range = max_dist / 3,
  nugget = var(train$zinc) * 0.2
)
# Fitting model semivariogram
fit_vg <- fit.variogram(
  vg,
  model = initial_model,
  fit.ranges = FALSE
)
# Melihat model semivariogram
print(fit_vg)
# Prediksi Ordinary Kriging
krig <- krige(
  zinc ~ 1,
  locations = train,
  newdata = test,
  model = fit_vg
)
# Mengambil hasil prediksi
p5 <- krig$var1.pred
# Menghitung RMSE
rmse5[k] <- RMSE(test$zinc, p5)
}

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
##   model     psill    range
## 1   Nug  40037.02    0.000
## 2   Sph 176463.57 1475.064
## [using ordinary kriging]

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
##   model     psill    range
## 1   Nug  33137.43    0.000
## 2   Sph 180641.90 1480.255
## [using ordinary kriging]

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
##   model     psill    range
## 1   Nug  40263.92    0.000
## 2   Sph 117318.87 1480.255
## [using ordinary kriging]

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
##   model     psill    range
## 1   Nug  39416.49    0.000
## 2   Sph 155578.55 1398.325
## [using ordinary kriging]

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
##   model     psill    range
## 1   Nug  48677.78    0.000
## 2   Sph 150178.15 1480.255
## [using ordinary kriging]

#RMSEobtainedforeachofthe5splits
data.frame(closest.obs= rmse1,IDW= rmse2,
nearest.neigh= rmse3,ensemble= rmse4, kriging= rmse5)

##   closest.obs      IDW nearest.neigh ensemble  kriging
## 1    309.9009 282.3744      265.3339 250.9364 222.5966
## 2    385.8027 277.6231      275.7933 260.7226 243.6972
## 3    323.5046 438.7651      364.2145 337.4435 305.2135
## 4    237.6379 323.2493      255.5252 242.7105 210.6315
## 5    299.4034 329.9387      294.7488 281.7473 266.0208

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

##   closest.obs      IDW nearest.neigh ensemble  kriging
## 1    311.2499 330.3901      291.1231 274.7121 249.6319

Hasil perhitungan RMSE pada setiap metode interpolasi menunjukkan bahwa metode kriging memiliki nilai RMSE yang paling kecil. Artinya metode kriging merupakan metode dengan kesalahan interpolasi paling kecil diantara ke-5 metode yang digunakan sehingga metode kriging merupakan metode interpolasi terbaik.