Regresi Spasial Pertemuan 9

# Memuat pustaka yang diperlukan
library(sp)        # untuk data meuse
## Warning: package 'sp' was built under R version 4.5.2
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.2
## terra 1.8.93
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
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
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
# Memuat data meuse
data(meuse)
data(meuse.area)

# Memuat data titik observasi (setara 'properties')
d <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)

# Menambahkan kolom variabel yang diprediksi (setara harga per m2)
d$vble <- d$zinc

# Membuat peta dari batas wilayah studi (setara batas administratif Athena)
map <- st_sfc(st_polygon(list(as.matrix(meuse.area))), crs = 28992) %>%
  st_sf()

# Menampilkan peta kadar zinc
tmap_mode("plot")
## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`
## This message is displayed once per session.
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`.

data(meuse.grid)
grid <- rast(meuse.grid[, c("x", "y", "dist")], type = "xyz", crs = "EPSG:28992")

Closest Observation

v <- terra::voronoi(x = vect(d), bnd = vect(map))
plot(v)
points(vect(d), cex = 0.5)

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 visual 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(st_make_valid(v), map)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
resp$pred <- resp$vble
pred <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
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
## visual variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
c1

# Inverse Distance Weight

library(gstat)
## Warning: package 'gstat' was built under R version 4.5.3
data(meuse.grid)
coop <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)

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

# 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
pred <- rasterize(vect(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
## visual 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 Neighbours

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(vect(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
## visual variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
c3

Pendekatan Ensemble

# Closest observation (Voronoi)
v <- terra::voronoi(x = vect(d), bnd = vect(map))
v <- st_make_valid(st_as_sf(v))
p1 <- st_join(coop, v["vble"])$vble

# 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]
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(vect(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
## visual variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
c4

tmap_arrange(
  c1 + tm_title("Closest Observation"),
  c2 + tm_title("IDW"),
  c3 + tm_title("Nearest Neighbours"),
  c4 + tm_title("Ensemble"),
  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 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")
# 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()
##    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))
fv
##   model      psill    range
## 1   Nug 0.05066017   0.0000
## 2   Sph 0.59060556 897.0044
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)
## 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 = 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()

Evaluasi Model

library(sf)
library(gstat)

set.seed(123)
k_fold <- 5
folds <- sample(rep(1:k_fold, length.out = nrow(d)))

# Tabel penampung prediksi
cv <- data.frame(obs = d$vble,
                 voronoi = NA, idw = NA, knn = NA, kriging = NA)

for (i in 1:k_fold) {
  train <- d[folds != i, ]
  test  <- d[folds == i, ]

  # 1. Closest observation (Voronoi)
  nearest <- st_nearest_feature(test, train)
  cv$voronoi[folds == i] <- train$vble[nearest]

  # 2. IDW
  gs <- gstat(formula = vble ~ 1, locations = train,
              nmax = nrow(train), set = list(idp = 1))
  cv$idw[folds == i] <- predict(gs, test, debug.level = 0)$var1.pred

  # 3. Nearest neighbours (k = 5)
  nn <- gstat(formula = vble ~ 1, locations = train,
              nmax = 5, set = list(idp = 0))
  cv$knn[folds == i] <- predict(nn, test, debug.level = 0)$var1.pred

  # 4. Kriging (log-zinc, lalu back-transform)
  kr <- gstat(formula = log(vble) ~ 1, data = train, model = fv)
  cv$kriging[folds == i] <- exp(predict(kr, test, debug.level = 0)$var1.pred)
}

# 5. Ensemble (rata-rata Voronoi, IDW, NN)
cv$ensemble <- (cv$voronoi + cv$idw + cv$knn) / 3
rmse <- function(obs, pred) sqrt(mean((obs - pred)^2))

hasil <- sapply(cv[, -1], function(p) rmse(cv$obs, p))
hasil <- sort(round(hasil, 2))
hasil
##  kriging      knn ensemble      idw  voronoi 
##   229.40   267.83   277.51   330.75   351.95
barplot(hasil, col = viridis::viridis(length(hasil)),
        ylab = "RMSE (zinc, ppm)", main = "Perbandingan RMSE (5-fold CV)",
        las = 2)

Berdasarkan 5-fold cross validation pada data meuse (variabel zinc), kriging menghasilkan akurasi terbaik dengan RMSE 229,40 ppm. Nilai ini sekitar 14% lebih rendah dari Nearest Neighbours (267,83) dan 35% lebih rendah dari Closest Observation/Voronoi (351,95).