Data

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

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

# Memuat data meuse
library(sp)
data(meuse)
data(meuse.grid)

# Mengubah data meuse menjadi objek sf (CRS 28992)
d <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)

# Menambahkan kolom log(zinc) ke dalam dataset
d$vble <- log(d$zinc)

# Membuat batas wilayah dari sebaran titik meuse
map <- st_buffer(st_convex_hull(st_union(d)), 100) %>% st_sf()

# Menampilkan peta log(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`.

# Membuat grid raster 100x100 sel untuk menutupi wilayah
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
coop <- st_filter(coop, map)

# Menampilkan peta lokasi prediksi
qtm(coop)

Metode Deteministik

# Membuat diagram Voronoi dari lokasi titik
v <- voronoi(x = vect(d), bnd = map)
plot(v)
# Menampilkan diagram Voronoi dengan titik observasi
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")

##

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

#IDW

library(gstat)

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

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

##

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

Pendekatan Ensembel

#Closestobservation(Voronoi)
v<-voronoi(x=vect(d),bnd= map)
v<-st_as_sf(v)
p1<-st_join(coop,v)$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]

#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

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

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 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<-d[kf== k,]
train<-d[kf !=k,]


#Closestobservation
v<-voronoi(x=vect(train), bnd= map)
v<-st_as_sf(v)
p1<-st_join(test,v)$vble.y
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)

#Nearestneighbors
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(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$vble,p4)

#Kriging
vtr<-variogram(vble~1,data=train)
ftr<-fit.variogram(vtr,model=vgm(psill=0.5,model="Sph",range=900,nugget=0.1))
kr<-gstat(formula=vble~1,data=train,model=ftr)
p5<-predict(kr,test)$var1.pred
rmse5[k] <-RMSE(test$vble,p5)
}

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]

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

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

Metode Stokastik

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

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

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

Tugas Praktikum Regresi Spasial

Dahayu Salsabila

10 November 2026

Data

Metode Deteministik

Nearest Neighbors

Pendekatan Ensembel

Cross Validation

Metode Stokastik