Tugas Praktikum Regspas Pertemuan 10

Persiapan Data

# Memuat pustaka yang diperlukan
library(gstat)
## Warning: package 'gstat' was built under R version 4.5.3
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.50
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.3
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.2
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.2
## 
## 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(automap)
## Warning: package 'automap' was built under R version 4.5.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.2
# Memuat dataset Jura
data(jura)

# Mengambil data pengamatan
d <- prediction.dat

# Mengubah data menjadi objek sf
d <- st_as_sf(d, coords = c("Xloc", "Yloc"))

# Menambahkan variabel yang akan diinterpolasi
d$vble <- d$Cd

Membuat peta dan grid prediksi

# Membuat batas wilayah berdasarkan lokasi pengamatan
map <- st_union(d) %>%
  st_convex_hull() %>%
  st_sf()

# Menampilkan peta lokasi pengamatan
tmap_mode("plot")
## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`
tm_shape(map) +
  tm_polygons(alpha = 0.3) +
  tm_shape(d) +
  tm_dots("vble", palette = "viridis")
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.

# Membuat grid raster 100x100 sel
grid <- rast(map, nrows = 100, ncols = 100)

# Mengambil koordinat titik pusat setiap sel grid
xy <- xyFromCell(grid, 1:ncell(grid))

# Membuat objek sf untuk lokasi prediksi
coop <- st_as_sf(
  as.data.frame(xy),
  coords = c("x", "y"),
  crs = st_crs(map)
)

# Menyaring lokasi prediksi yang berada dalam batas wilayah
coop <- st_filter(coop, map)

# Menampilkan lokasi prediksi
qtm(coop)

Closest Observation

# Mengubah batas wilayah menjadi SpatVector
map_vect <- terra::vect(map)

# Membuat diagram Voronoi
v <- terra::voronoi(
  x = terra::vect(d),
  bnd = map_vect
)

# Menampilkan diagram Voronoi
terra::plot(v)

# Prediction
v <- st_as_sf(v)
tm_shape(v) + tm_fill(col = "vble", 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`.

# Mengubah objek Voronoi menjadi sf
v <- sf::st_as_sf(v)

# Menggabungkan poligon Voronoi dengan titik grid
resp <- sf::st_intersection(v, coop)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
# Menambahkan nilai prediksi
resp$pred <- resp$vble

# Mengubah hasil intersection menjadi SpatVector
resp_vect <- terra::vect(resp)

# Melakukan rasterisasi
pred <- terra::rasterize(
  resp_vect,
  grid,
  field = "pred",
  fun = "mean"
)

# Menampilkan hasil interpolasi Closest Observation
c1 <- 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`.
c1

Hasil interpolasi menunjukkan bahwa sebagian besar wilayah memiliki nilai Cd relatif rendah, sedangkan nilai Cd yang lebih tinggi hanya ditemukan pada beberapa lokasi tertentu. Hasil rasterisasi mempertahankan pola persebaran nilai tersebut dalam bentuk grid.

Invers Distance Weight

# 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)
## [inverse distance weighted interpolation]
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")
## 
## ── 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

Hasil interpolasi menunjukkan bahwa sebagian besar wilayah memiliki nilai Cd pada kisaran 1.0–1.4. Nilai Cd yang lebih tinggi, yaitu 1.5–3.0, cenderung ditemukan di bagian utara wilayah, sedangkan nilai yang lebih rendah tersebar pada beberapa lokasi tertentu. Pola ini menunjukkan adanya variasi nilai Cd secara spasial dengan konsentrasi nilai yang relatif lebih tinggi di bagian utara.

Nearest Neighbors

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

resp <- predict(res, coop)
## [inverse distance weighted interpolation]
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")
## 
## ── 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

Hasil interpolasi menunjukkan bahwa sebagian besar wilayah memiliki nilai Cd pada kisaran 1.0–1.9. Sementara itu, nilai Cd yang lebih tinggi tersebar pada beberapa lokasi, terutama di bagian utara dan tengah wilayah. Hal ini menunjukkan adanya variasi nilai Cd secara spasial pada wilayah pengamatan.

Pendekatan Ensembel

#Closestobservation(Voronoi)
v<-voronoi(x=vect(d),bnd= map)
v<-st_as_sf(v)
p1<-st_intersection(v,coop)$vble
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
#IDW
gs<-gstat(formula=vble~1,locations= d,nmax= nrow(d),
set= list(idp=1))
p2<-predict(gs,coop)$var1.pred
## [inverse distance weighted interpolation]
#Nearestneighbors
nn<-gstat(formula=vble~1,locations= d,nmax= 5,
set= list(idp=0))
p3<-predict(nn,coop)$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]
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")
## 
## ── 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

Hasil interpolasi menunjukkan bahwa sebagian besar wilayah memiliki nilai Cd pada kisaran 0.5–1.9. Nilai Cd yang lebih tinggi, yaitu 2.0–4.0, cenderung terkonsentrasi di bagian utara dan beberapa lokasi di bagian tengah wilayah. Sementara itu, nilai Cd yang lebih rendah tersebar di beberapa bagian wilayah pengamatan.

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

Kriging

# Memuat paket yang dibutuhkan
library(sp)
library(gstat)
library(sf)
library(mapview)
## Warning: package 'mapview' was built under R version 4.5.3
# Memuat dataset jura
data(jura)

# Mengonversi dataset jura menjadi objek sf
jura.pred <- st_as_sf(jura.pred, coords = c("long", "lat"), crs = 4326)
jura.grid <- st_as_sf(jura.grid, coords = c("long", "lat"), crs = 4326)

# Menampilkan peta interaktif
mapview(jura.pred, zcol = "Cd", map.types = "CartoDB.Voyager")
# Menampilkan peta jura.grid
mapview(jura.grid, map.types = "CartoDB.Voyager")
# Variogram dengan cloud plot
vc <- variogram(log(vble) ~ 1, d, cloud = TRUE)
plot(vc)

Berdasarkan grafik variogram cloud, sebagian besar pasangan titik memiliki nilai semivarian yang rendah pada jarak dekat, menunjukkan bahwa lokasi yang berdekatan cenderung memiliki nilai Cd yang mirip. Seiring bertambahnya jarak antarobservasi, sebaran nilai semivarian semakin luas dan beberapa pasangan titik menunjukkan nilai semivarian yang tinggi. Hal ini mengindikasikan bahwa kemiripan nilai Cd cenderung menurun ketika jarak antar lokasi semakin jauh, sehingga terdapat indikasi ketergantungan spasial pada data.

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

Berdasarkan grafik variogram standar, nilai semivarian meningkat seiring bertambahnya jarak antarobservasi hingga mencapai nilai tertinggi sekitar jarak 1.2. Setelah itu, semivarian cenderung stabil dan sedikit menurun. Hal ini menunjukkan adanya ketergantungan spasial pada data Cd, dengan kemiripan nilai antar lokasi yang semakin berkurang seiring bertambahnya jarak.

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

Visualisasi menunjukkan berbagai bentuk model variogram teoretis yang dapat digunakan dalam analisis geostatistika, seperti Nugget, Exponential, Spherical, Gaussian, dan model lainnya. Setiap model memiliki karakteristik pola semivarian yang berbeda sehingga dapat dipilih untuk menyesuaikan variogram empiris berdasarkan data pengamatan.

# Menyesuaikan variogram dengan model Spherical
fv <- fit.variogram(
  object = v,
  model = vgm(
    psill = 0.5,
    model = "Sph",
    range = 1.2,
    nugget = 0.1
  )
)

# Menampilkan hasil fitting
fv
##   model     psill   range
## 1   Nug 0.2103211 0.00000
## 2   Sph 0.3191918 1.04677

Berdasarkan hasil fitting model variogram Spherical, diperoleh nilai nugget sebesar 0.2103, partial sill sebesar 0.3192, dan range sebesar 1.0468. Nilai sill total diperoleh dari penjumlahan nugget dan partial sill, yaitu sebesar 0.5295. Hal ini menunjukkan bahwa ketergantungan spasial pada data Cd berlangsung hingga jarak sekitar 1.0468. Setelah melewati jarak tersebut, semivarian cenderung mencapai kondisi stabil.

plot(v,fv, cex= 1.5)

Berdasarkan grafik, model variogram Spherical cukup baik dalam mengikuti pola variogram empiris. Nilai semivarian meningkat seiring bertambahnya jarak hingga mencapai kondisi stabil (sill) sekitar 0.5295 pada jarak (range) 1.0468. Hal ini menunjukkan bahwa ketergantungan spasial antar lokasi semakin melemah setelah jarak tersebut.

library(ggplot2)
library(viridis)

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

# Melakukan prediksi variabel pada grid menggunakan model Kriging
kpred <- predict(k, coop)
## [using ordinary kriging]
# Visualisasi prediksi log(Cd) dengan warna berdasarkan var1.pred
ggplot() +
  geom_sf(data = kpred, aes(color = var1.pred)) +
  geom_sf(data = d) +
  scale_color_viridis(name = "log(Cd)") +
  theme_bw()

Berdasarkan hasil visualisasi Ordinary Kriging, nilai prediksi Cd menunjukkan variasi spasial di wilayah pengamatan. Area berwarna kuning hingga hijau menunjukkan nilai Cd yang relatif tinggi, sedangkan area berwarna biru hingga ungu menunjukkan nilai yang lebih rendah. Nilai Cd yang tinggi terlihat terkonsentrasi pada beberapa lokasi di bagian utara dan tengah wilayah, sementara nilai yang lebih rendah tersebar di beberapa bagian lainnya. Titik hitam menunjukkan lokasi pengamatan aktual yang digunakan dalam proses interpolasi.

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

Berdasarkan peta varians prediksi Kriging, sebagian besar wilayah memiliki nilai varians sekitar 0.25–0.35 yang menunjukkan tingkat ketidakpastian prediksi relatif rendah. Sementara itu, varians yang lebih tinggi (hijau–kuning) terlihat pada beberapa bagian tepi wilayah, terutama di sisi timur. Hal ini menunjukkan bahwa ketidakpastian hasil prediksi cenderung lebih tinggi pada lokasi tersebut.

Evaluasi Model

Membandingkan RMSE dengan menggunakan Cross Validation

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
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){

  # Split data into test and train
  test <- d[kf == k,]
  train <- d[kf != k,]

  # Closest observation
  v <- voronoi(x = vect(train), bnd = map)
  v <- st_as_sf(v)
  p1 <- st_intersection(v, test)$vble
  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)$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)$var1.pred
  rmse3[k] <- RMSE(test$vble, p3)

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

  # Kriging
  vg_train <- variogram(log(vble) ~ 1, data = train)

  fv_train <- fit.variogram(
    vg_train,
    model = vgm(psill = 0.5, model = "Sph",
                range = 1.2, nugget = 0.1)
  )

  kpred_cv <- krige(
    log(vble) ~ 1,
    locations = train,
    newdata = test,
    model = fv_train
  )

  p5 <- exp(kpred_cv$var1.pred)
  rmse5[k] <- RMSE(test$vble, p5)
}
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
# Menampilkan RMSE setiap fold
data.frame(
  closest.obs = rmse1,
  IDW = rmse2,
  nearest.neigh = rmse3,
  ensemble = rmse4,
  kriging = rmse5
)
##   closest.obs       IDW nearest.neigh  ensemble   kriging
## 1   1.1044389 0.7560668     0.8515053 0.8142814 0.8271888
## 2   1.0748607 0.6628500     0.6750161 0.6934629 0.6423065
## 3   0.9867737 0.6956564     0.7882901 0.6580956 0.7006500
## 4   1.0648430 0.7954476     0.9217603 0.8465482 0.8354231
## 5   1.0433718 0.9310081     1.0405151 0.9464067 1.0210067
# Menghitung rata-rata RMSE setiap metode
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    1.054858 0.7682058     0.8554174 0.791759 0.805315

Berdasarkan hasil 5-fold Cross Validation, metode IDW menghasilkan nilai RMSE rata-rata terkecil sebesar 0.7682, sehingga memiliki tingkat akurasi prediksi terbaik dibandingkan metode lainnya. Metode Ensemble memperoleh RMSE sebesar 0.7918, diikuti Kriging sebesar 0.8053, Nearest Neighbors sebesar 0.8554, dan Closest Observation sebesar 1.0549.

Dengan demikian, IDW merupakan metode interpolasi terbaik untuk memprediksi nilai Cd pada dataset Jura karena menghasilkan kesalahan prediksi paling rendah. Sebaliknya, Closest Observation memiliki performa paling rendah karena menghasilkan nilai RMSE terbesar.