Tugas Pertemuan 9: Interpolasi Spasial Data Coalash

Metode yang digunakan sesuai ketentuan tugas: Closest Observation, Nearest Neighbors (NN), Inverse Distance Weighting (IDW), Ensembel (minimal 3 metode), dan Kriging, lalu RMSE masing-masing dibandingkan dengan K-fold Cross Validation.

Data yang Digunakan

Data coalash berupa data.frame dengan koordinat lokal x, y dan peubah coalash. Karena koordinatnya lokal (bukan lintang/bujur), objek sf tidak diberi CRS. Hal ini tidak masalah karena interpolasi hanya membutuhkan jarak antartitik. Peubah coalash disalin ke kolom vble agar sama dengan penamaan pada materi.

# Memuat pustaka yang diperlukan
library(sp)
library(gstat)
library(sf)
library(terra)
library(tmap)
library(viridis)
library(dismo)
library(dplyr)
library(ggplot2)
library(automap)

# Memuat data coalash
data(coalash)
str(coalash)

## 'data.frame':    208 obs. of  3 variables:
##  $ x      : int  1 1 1 2 2 2 2 2 2 2 ...
##  $ y      : int  14 15 16 8 10 11 12 14 15 16 ...
##  $ coalash: num  10.21 9.92 11.17 10.01 11.15 ...

# Mengubah menjadi objek sf (tanpa CRS)
d <- st_as_sf(coalash, coords = c("x", "y"))

# Menghapus titik dengan koordinat kembar (mencegah matriks singular pada Kriging)
d <- d[!duplicated(st_coordinates(d)), ]

# Menambahkan kolom vble untuk peubah yang diinterpolasi
d$vble <- d$coalash

nrow(d)

## [1] 208

summary(d$vble)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   7.000   8.960   9.785   9.779  10.568  17.610

# Membuat batas area: convex hull dari seluruh titik observasi
map <- st_convex_hull(st_union(d)) %>% st_sf()

# Menampilkan peta kadar abu
tmap_mode("plot")
tm_shape(map) + tm_polygons(alpha = 0.3) + tm_shape(d) +
  tm_dots("vble", palette = "viridis")

map adalah batas area (convex hull titik observasi) yang digunakan untuk membatasi area prediksi.

Grid Lokasi Prediksi

Data coalash tidak memiliki grid prediksi bawaan, sehingga dibuat grid raster 100 x 100 sel yang menutupi area observasi. Titik pusat tiap sel menjadi lokasi prediksi.

# Membuat grid raster 100x100 sel
grid <- rast(vect(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 area
coop <- st_filter(coop, map)

# Menampilkan peta lokasi prediksi
qtm(coop)

Fungsi Bantu

Beberapa fungsi bantu dibuat agar langkah yang sama dapat dipakai ulang pada prediksi grid dan pada cross validation.

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

# Closest observation (Voronoi): nilai diambil dari poligon Voronoi yang memuat lokasi.
# st_join() menjaga urutan lokasi prediksi; baris ganda (titik tepat di batas poligon) dibuang.
pred_voronoi <- function(train, new) {
  v <- st_as_sf(voronoi(x = vect(train), bnd = vect(map)))
  new_g <- st_sf(geometry = st_geometry(new))
  r <- st_join(new_g, v["vble"])
  r <- r[!duplicated(st_coordinates(r)), ]
  r$vble
}

# IDW dan Nearest Neighbors (via gstat)
pred_gstat <- function(train, new, nmax, idp) {
  g <- gstat(formula = vble ~ 1, locations = train,
             nmax = nmax, set = list(idp = idp))
  predict(g, new, debug.level = 0)$var1.pred
}

# Model variogram awal berdasarkan skala data
# (skala koordinat coalash berbeda dengan data berkoordinat meter, jadi tidak memakai angka tetap)
vgm_awal <- function(train) {
  xy <- st_coordinates(train)
  diag_len <- sqrt(diff(range(xy[, 1]))^2 + diff(range(xy[, 2]))^2)
  vgm(psill = 0.8 * var(train$vble), model = "Sph",
      range = diag_len / 3, nugget = 0.1 * var(train$vble))
}

# Fit variogram; jika hasil singular, pakai model awal
fit_vgm <- function(train) {
  v    <- variogram(vble ~ 1, train)
  init <- vgm_awal(train)
  fv   <- suppressWarnings(fit.variogram(v, model = init))
  if (isTRUE(attr(fv, "singular"))) fv <- init
  fv
}

# Ordinary Kriging
pred_ok <- function(train, new) {
  fv <- fit_vgm(train)
  krige(vble ~ 1, train, new, model = fv, debug.level = 0)$var1.pred
}

# Mengubah vektor prediksi pada lokasi coop menjadi raster
buat_raster <- function(p) {
  resp <- coop
  resp$pred <- p
  rasterize(vect(resp), grid, field = "pred", fun = "mean")
}

Metode Deterministik

Closest Observation

Diagram Voronoi membagi ruang menjadi poligon yang masing-masing memuat satu titik observasi. Setiap lokasi di dalam poligon lebih dekat ke titik pembentuknya daripada ke titik lain, sehingga prediksinya adalah nilai titik terdekat.

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

# Prediksi
v <- st_as_sf(v)
tm_shape(v) + tm_fill(col = "vble", palette = "viridis", alpha = 0.6)

Setiap poligon Voronoi mewakili wilayah prediksi yang dipengaruhi oleh titik terdekatnya.

p1 <- pred_voronoi(d, coop)
c1 <- tm_shape(buat_raster(p1)) + tm_raster(alpha = 0.6, palette = "viridis")
c1

Inverse Distance Weighting (IDW)

Metode IDW memberi bobot berdasarkan jarak antara lokasi yang belum teramati dan lokasi yang sudah teramati. Semakin dekat jaraknya, semakin besar bobotnya. Digunakan semua titik observasi (nmax = nrow(d)) dan power idp = 2 (nilai ini dapat diubah).

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

# 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
p2 <- resp$pred

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

Nearest Neighbors

Metode ini menghitung rata-rata nilai dari k = 5 titik terdekat. Setiap tetangga dianggap memiliki bobot yang sama (idp = 0).

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
p3 <- resp$pred

pred <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
c3 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
c3

Metode Stokastik: Kriging

Ordinary Kriging (OK) memperhitungkan autokorelasi spasial dan mengasumsikan rata-rata populasi tidak diketahui tetapi konstan di seluruh wilayah. Langkahnya: membangun variogram empiris, men-fit model semivariogram, lalu memprediksi lokasi baru.

Keterangan parameter semivariogram:

Nugget: semivarian saat jarak mendekati nol (variasi mikrospasial dan kesalahan pengukuran).
Sill: nilai maksimum (plateau) semivarian, yaitu total varians spasial.
Range: jarak saat semivarian mencapai sill, yaitu batas jangkauan ketergantungan spasial.

Visualisasi interaktif dengan mapview pada materi tidak dipakai di sini karena koordinat coalash bersifat lokal sehingga tidak dapat ditumpangkan pada peta dasar.

Variogram

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

# Daftar model variogram yang tersedia
show.vgms(par.strip.text = list(cex = 0.75))

Fitting Model Semivariogram

# Model variogram awal (psill, model, range, nugget dihitung dari skala data)
vinitial <- vgm_awal(d)
vinitial

##   model     psill    range
## 1   Nug 0.1629284 0.000000
## 2   Sph 1.3034271 8.875685

# Memvisualisasikan variogram dengan model awal
plot(v, vinitial, cex = 1.5)

# Menyesuaikan variogram dengan model "Sph"
fv <- fit.variogram(object = v, model = vinitial)
if (isTRUE(attr(fv, "singular"))) fv <- vinitial   # pengaman jika hasil singular
fv

##   model     psill    range
## 1   Nug 1.0838788 0.000000
## 2   Sph 0.4837233 8.390055

plot(v, fv, cex = 1.5)

Interpretasi: nugget = 1.0839, partial sill = 0.4837, dan range = 8.39 (dalam satuan koordinat). Artinya ketergantungan spasial kadar abu berlangsung hingga jarak sekitar 8.39; di luar jarak itu observasi relatif tidak berkorelasi.

Sebagai pembanding, autofitVariogram() dari paket automap memilih model terbaik secara otomatis:

# Autofit variogram untuk data kadar abu
afvg <- autofitVariogram(vble ~ 1, as_Spatial(d))

# Plot hasil
plot(afvg)

Prediksi Kriging

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

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

## [using ordinary kriging]

p5 <- kpred$var1.pred

# Visualisasi prediksi dengan warna berdasarkan var1.pred
ggplot() + geom_sf(data = kpred, aes(color = var1.pred)) +
  geom_sf(data = d) +
  scale_color_viridis(name = "Kadar abu") + coord_sf(datum = NA) + theme_bw()

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

Pendekatan Ensembel

Pendekatan ensembel menggabungkan prediksi dari beberapa metode dengan bobot tertentu. Pada tugas ini dibuat tiga ensembel, masing-masing menggabungkan 3 metode:

Ensembel 1 (bobot sama): Closest Observation + IDW + NN, bobot masing-masing 1/3.
Ensembel 2 (bobot invers RMSE): Closest Observation + IDW + NN, metode dengan RMSE lebih kecil diberi bobot lebih besar. Pada cross validation, bobot dihitung dari leave-one-out di dalam data latih agar data uji tidak ikut menentukan bobot.
Ensembel 3 (bobot sama, melibatkan Kriging): IDW + NN + Kriging, bobot masing-masing 1/3.

# Ensembel 1 (bobot sama)
weights <- c(1/3, 1/3, 1/3)
p4 <- p1 * weights[1] + p2 * weights[2] + p3 * weights[3]
c4 <- tm_shape(buat_raster(p4)) + tm_raster(alpha = 0.6, palette = "viridis")

# Ensembel 3 (IDW + NN + Kriging, bobot sama)
p6 <- (p2 + p3 + p5) / 3
c6 <- tm_shape(buat_raster(p6)) + tm_raster(alpha = 0.6, palette = "viridis")

# Peta Kriging dalam bentuk raster
c5 <- tm_shape(buat_raster(p5)) + tm_raster(alpha = 0.6, palette = "viridis")

Gabungan seluruh peta dalam satu frame (urutan: Closest Observation, IDW, NN, Ensembel 1, Kriging, Ensembel 3):

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

Evaluasi Model

Untuk menilai kinerja tiap metode digunakan 5-fold cross validation dan RMSE. Seluruh proses (termasuk fitting variogram pada Kriging) dilakukan hanya dengan data latih pada setiap lipatan.

set.seed(123)

# Pembagian data untuk K-fold
kf <- kfold(nrow(d), k = 5)

# RMSE leave-one-out di dalam data latih (untuk bobot Ensembel 2)
loo_rmse <- function(train) {
  n <- nrow(train)
  loo <- sapply(1:n, function(i) {
    c(closest = pred_gstat(train[-i, ], train[i, ], 1,     1),
      idw     = pred_gstat(train[-i, ], train[i, ], n - 1, 2),
      nn      = pred_gstat(train[-i, ], train[i, ], 5,     0))
  })
  apply(loo, 1, function(p) RMSE(train$vble, p))
}

# Vektor untuk menyimpan nilai RMSE tiap metode
rmse1 <- rep(NA, 5) # Closest observation
rmse2 <- rep(NA, 5) # IDW
rmse3 <- rep(NA, 5) # Nearest neighbors
rmse4 <- rep(NA, 5) # Ensembel 1 (bobot sama: Closest + IDW + NN)
rmse5 <- rep(NA, 5) # Kriging
rmse6 <- rep(NA, 5) # Ensembel 2 (bobot invers RMSE: Closest + IDW + NN)
rmse7 <- rep(NA, 5) # Ensembel 3 (bobot sama: IDW + NN + Kriging)

pred_cv <- vector("list", 5)

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

  # Closest observation
  p1 <- pred_voronoi(train, test)
  rmse1[k] <- RMSE(test$vble, p1)

  # IDW
  gs <- gstat(formula = vble ~ 1, locations = train,
              nmax = nrow(train), set = list(idp = 2))
  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)

  # Kriging
  p5 <- pred_ok(train, test)
  rmse5[k] <- RMSE(test$vble, p5)

  # Ensembel 1 (bobot sama)
  p4 <- (p1 + p2 + p3) / 3
  rmse4[k] <- RMSE(test$vble, p4)

  # Ensembel 2 (bobot invers RMSE dari LOO pada data latih: RMSE kecil, bobot besar)
  w <- 1 / loo_rmse(train)[c("closest", "idw", "nn")]
  weights <- w / sum(w)
  p6 <- p1 * weights[1] + p2 * weights[2] + p3 * weights[3]
  rmse6[k] <- RMSE(test$vble, p6)

  # Ensembel 3 (IDW + NN + Kriging, bobot sama)
  p7 <- (p2 + p3 + p5) / 3
  rmse7[k] <- RMSE(test$vble, p7)

  pred_cv[[k]] <- data.frame(obs = test$vble, Closest = p1, IDW = p2, NN = p3,
                             Ensembel1 = p4, Kriging = p5,
                             Ensembel2 = p6, Ensembel3 = p7)
}

# RMSE yang diperoleh pada masing-masing dari 5 split
data.frame(closest.obs = rmse1, IDW = rmse2, nearest.neigh = rmse3,
           ensemble1 = rmse4, kriging = rmse5,
           ensemble2 = rmse6, ensemble3 = rmse7)

##   closest.obs       IDW nearest.neigh ensemble1   kriging ensemble2 ensemble3
## 1    1.338597 0.8966235     0.9922960 0.9666776 0.8723656  0.949733 0.8983698
## 2    1.503519 1.0835911     1.2366343 1.2042684 1.2166928  1.186589 1.1642250
## 3    1.294278 1.0091334     0.9548066 1.0176335 0.9883335  1.003835 0.9690469
## 4    1.546807 1.0591014     1.2108443 1.1734685 1.0668098  1.152860 1.0971449
## 5    1.575616 1.5043077     1.5445686 1.4977107 1.4543084  1.497754 1.4922017

# Rata-rata RMSE dari 5 split
data.frame(closest.obs = mean(rmse1), IDW = mean(rmse2),
           nearest.neigh = mean(rmse3), ensemble1 = mean(rmse4),
           kriging = mean(rmse5), ensemble2 = mean(rmse6),
           ensemble3 = mean(rmse7))

##   closest.obs      IDW nearest.neigh ensemble1  kriging ensemble2 ensemble3
## 1    1.451763 1.110551       1.18783  1.171952 1.119702  1.158154  1.124198

Perbandingan RMSE

rmse_rata <- data.frame(
  Metode = c("Closest Observation", "IDW", "Nearest Neighbors",
             "Ensembel 1 (C+IDW+NN)", "Kriging",
             "Ensembel 2 (C+IDW+NN, bobot invers RMSE)",
             "Ensembel 3 (IDW+NN+Kriging)"),
  RMSE = c(mean(rmse1), mean(rmse2), mean(rmse3), mean(rmse4),
           mean(rmse5), mean(rmse6), mean(rmse7))
) %>% arrange(RMSE)

rmse_rata

##                                     Metode     RMSE
## 1                                      IDW 1.110551
## 2                                  Kriging 1.119702
## 3              Ensembel 3 (IDW+NN+Kriging) 1.124198
## 4 Ensembel 2 (C+IDW+NN, bobot invers RMSE) 1.158154
## 5                    Ensembel 1 (C+IDW+NN) 1.171952
## 6                        Nearest Neighbors 1.187830
## 7                      Closest Observation 1.451763

ggplot(rmse_rata, aes(x = reorder(Metode, RMSE), y = RMSE, fill = RMSE)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = round(RMSE, 3)), hjust = -0.1, size = 3.5) +
  coord_flip(ylim = c(0, max(rmse_rata$RMSE) * 1.15)) +
  scale_fill_viridis(direction = -1) +
  theme_bw() +
  labs(x = "Metode", y = "Rata-rata RMSE (5-fold CV)",
       title = "Perbandingan RMSE antarmetode")

Kesimpulan

Berdasarkan rata-rata RMSE dari 5-fold cross validation, metode dengan kinerja terbaik adalah IDW (RMSE = 1.1106), sedangkan yang terburuk adalah Closest Observation (RMSE = 1.4518).