Tugas Praktikum 9 Regspas

1 Interpolasi Spasial

Interpolasi spasial merupakan metode yang digunakan untuk memperkirakan nilai peubah spasial pada lokasi yang tidak diketahui berdasarkan data yang tersedia pada lokasi tetangga yang diketahui (Cressie, 1993).

Pada tugas ini digunakan data meuse dari paket sp yang berisi konsentrasi logam berat di tanah pada dataran banjir Sungai Meuse, Belanda. Peubah yang diinterpolasi adalah kadar zinc (zinc, ppm). Metode yang dibandingkan:

1. Closest Observation (Voronoi)
2. Inverse Distance Weighting (IDW)
3. Nearest Neighbors (NN)
4. Ensembel (Closest Observation + IDW + NN)
5. Ordinary Kriging

Seluruh metode dibandingkan menggunakan RMSE dari K-fold Cross Validation.

2 Data yang digunakan

# Memuat pustaka yang diperlukan
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.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(sp)
library(gstat)

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

library(mapview)

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

# dataset meuse dan meuse.grid dari R
data(meuse)
data(meuse.grid)

# Mengonversi menjadi objek 'sf' dengan CRS 28992 (Amersfoort / RD New)
meuse <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)
meuse.grid <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)

Peta interaktif data zinc dan grid lokasi prediksi:

mapview(meuse, zcol = "zinc", map.types = "CartoDB.Voyager")

mapview(meuse.grid, map.types = "CartoDB.Voyager")

# Memuat data titik pengamatan
d <- meuse

# Menambahkan kolom peubah yang akan diinterpolasi (kadar zinc)
d$vble <- d$zinc

# Membuat batas wilayah studi dari meuse.grid (setiap titik grid berukuran 40 m x 40 m)
map <- meuse.grid %>%
  st_buffer(dist = 20, endCapStyle = "SQUARE") %>%
  st_union() %>%
  st_sf()

# Membuat batas persegi panjang (bounding box) untuk diagram Voronoi
bb <- st_as_sfc(st_bbox(st_buffer(st_union(st_geometry(d)), 500))) %>% st_sf()

# Menampilkan peta kadar zinc
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`.

Untuk memprediksi kadar zinc, dibuat grid raster yang mencakup seluruh area dan titik pusat grid digunakan sebagai lokasi prediksi.

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

# Mengambil koordinat titik pusat dari setiap sel dalam 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 studi
coop <- st_filter(coop, map)

# Penanda urutan lokasi prediksi (agar urutan hasil Voronoi sama dengan urutan coop)
coop$id <- 1:nrow(coop)

# Menampilkan peta lokasi prediksi
qtm(coop)

3 Metode Deterministik

3.1 Closest Observation

Metode ini menggunakan diagram Voronoi untuk membagi ruang menjadi poligon yang masing-masing berisi satu titik pengamatan. Prediksi pada suatu lokasi adalah nilai titik pengamatan terdekat.

# Membuat diagram Voronoi dari lokasi pengamatan
v <- voronoi(x = vect(d), bnd = bb)
plot(v)
# Menampilkan titik pengamatan
points(vect(d), cex = 0.5)

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

resp <- st_intersection(v, coop)

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

resp$pred <- resp$vble
pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c1 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis") +
  tm_layout(main.title = "Closest Observation", main.title.size = 0.9)

## 

## ── 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`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

c1

3.2 Inverse Distance Weighting

Metode IDW memberikan bobot berdasarkan jarak antara lokasi yang belum teramati dan lokasi yang sudah teramati. Semakin dekat jarak, semakin besar bobotnya.

# Menerapkan IDW menggunakan semua lokasi pengamatan
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)

## [inverse distance weighted interpolation]

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

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

## 

## ── 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`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

c2

3.3 Nearest Neighbors

Metode Nearest Neighbors menghitung rata-rata nilai dari k titik terdekat (di sini k = 5). Setiap tetangga memiliki bobot yang sama.

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") +
  tm_layout(main.title = "Nearest Neighbors", main.title.size = 0.9)

## 

## ── 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`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

c3

3.4 Pendekatan Ensembel

Pendekatan ensembel menggabungkan hasil prediksi dari tiga metode (Closest Observation, IDW, dan Nearest Neighbors) dengan bobot tertentu.

# Closest observation (Voronoi)
v <- voronoi(x = vect(d), bnd = bb)
v <- st_as_sf(v)
tmp <- st_intersection(v, coop)

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

p1 <- tmp$vble[match(coop$id, tmp$id)]

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

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

## [inverse distance weighted interpolation]

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

## 

## ── 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`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

c4

4 Metode Stokastik

4.1 Kriging

Kriging adalah metode interpolasi geostatistik yang memperhitungkan autokorelasi spasial untuk membuat prediksi. Ordinary Kriging (OK) mengasumsikan nilai rata-rata populasi tidak diketahui, tetapi konstan di seluruh wilayah. Pada bagian ini variabel yang dimodelkan adalah log(zinc).

Keterangan:

• Nugget: Nilai semivarian ketika jarak antarobservasi mendekati nol, yang merepresentasikan variasi mikrospasial dan kesalahan pengukuran.
• Sill: Nilai maksimum atau plateau semivarian yang merepresentasikan total varians spasial.
• Range: Jarak ketika semivarian mencapai sill, yang menunjukkan batas jangkauan ketergantungan spasial.

Secara umum:

• Semivarian rendah pada jarak dekat → indikasi bahwa lokasi berdekatan memiliki nilai yang mirip.
• Semivarian meningkat seiring jarak → indikasi adanya struktur ketergantungan spasial.
• Semivarian mencapai sill → variasi spasial relatif stabil setelah jarak tertentu.

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

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

vgm()

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

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(log(zinc) ~ 1, meuse)

# Plot hasil
plot(afvg)

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

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

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

5 Evaluasi Model

Untuk menilai kinerja setiap metode dilakukan K-fold cross-validation (K = 5) dan dihitung RMSE (Root Mean Square Error).

set.seed(123)

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

# Fungsi untuk fit model semivariogram log(zinc) (nilai awal sama dengan bagian Kriging)
fit_vgm <- function(dat) {
  ve <- variogram(log(zinc) ~ 1, dat)
  init <- vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1)
  f <- tryCatch(fit.variogram(ve, init), error = function(e) init)
  if (isTRUE(attr(f, "singular"))) f <- init
  f
}

# 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
rmse5 <- rep(NA, 5) # Kriging

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

  # Closest observation
  v <- voronoi(x = vect(train), bnd = bb)
  v <- st_as_sf(v)
  tmp <- st_intersection(v, test)
  p1 <- tmp$vble[match(test$id, tmp$id)]
  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 (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)

  # Ordinary Kriging pada log(zinc), lalu di-back-transform ke skala zinc
  fit_k <- fit_vgm(train)
  p5 <- exp(krige(log(zinc) ~ 1, locations = train, newdata = test, model = fit_k)$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]

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

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

5.1 Perbandingan RMSE

hasil <- data.frame(
  Metode = c("Closest Observation", "IDW", "Nearest Neighbors", "Ensemble", "Kriging"),
  RMSE = c(mean(rmse1), mean(rmse2), mean(rmse3), mean(rmse4), mean(rmse5))
)
hasil <- hasil[order(hasil$RMSE), ]
knitr::kable(hasil, digits = 2, row.names = FALSE,
             caption = "Rata-rata RMSE 5-fold Cross Validation (urut dari terkecil)")

Rata-rata RMSE 5-fold Cross Validation (urut dari terkecil)

Metode	RMSE
Kriging	238.30
Ensemble	274.71
Nearest Neighbors	291.12
Closest Observation	311.25
IDW	330.39

barplot(setNames(hasil$RMSE, hasil$Metode), col = viridis(nrow(hasil)),
        ylab = "RMSE", main = "Perbandingan RMSE (5-fold CV)", las = 2, cex.names = 0.8)

5.2 Kesimpulan

Berdasarkan rata-rata RMSE dari 5-fold cross validation, metode dengan kinerja terbaik (RMSE terkecil) adalah Kriging dengan RMSE sebesar 238.3, sehingga menjadi metode interpolasi terbaik untuk memprediksi kadar zinc pada data meuse. sedangkan metode dengan RMSE terbesar adalah IDW dengan RMSE sebesar 330.39.