Tugas Prak 10
Annisa Marwah
2026-11-10
Library
library(sp)## Warning: package 'sp' was built under R version 4.5.3library(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 TRUElibrary(gstat)## Warning: package 'gstat' was built under R version 4.5.3library(automap)## Warning: package 'automap' was built under R version 4.5.3## Registered S3 method overwritten by 'stars':
## method from
## st_interpolate_aw.stars sflibrary(ggplot2)## Warning: package 'ggplot2' was built under R version 4.5.3library(dplyr)## Warning: package 'dplyr' was built under R version 4.5.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(viridis)## Warning: package 'viridis' was built under R version 4.5.3## Loading required package: viridisLitelibrary(patchwork)## Warning: package 'patchwork' was built under R version 4.5.3Data
data(meuse)
data(meuse.grid)
head(meuse)## x y cadmium copper lead zinc elev dist om ffreq soil lime
## 1 181072 333611 11.7 85 299 1022 7.909 0.00135803 13.6 1 1 1
## 2 181025 333558 8.6 81 277 1141 6.983 0.01222430 14.0 1 1 1
## 3 181165 333537 6.5 68 199 640 7.800 0.10302900 13.0 1 1 1
## 4 181298 333484 2.6 81 116 257 7.655 0.19009400 8.0 1 2 0
## 5 181307 333330 2.8 48 117 269 7.480 0.27709000 8.7 1 2 0
## 6 181390 333260 3.0 61 137 281 7.791 0.36406700 7.8 1 2 0
## landuse dist.m
## 1 Ah 50
## 2 Ah 30
## 3 Ah 150
## 4 Ga 270
## 5 Ah 380
## 6 Ga 470str(meuse)## 'data.frame': 155 obs. of 14 variables:
## $ x : num 181072 181025 181165 181298 181307 ...
## $ y : num 333611 333558 333537 333484 333330 ...
## $ cadmium: num 11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
## $ copper : num 85 81 68 81 48 61 31 29 37 24 ...
## $ lead : num 299 277 199 116 117 137 132 150 133 80 ...
## $ zinc : num 1022 1141 640 257 269 ...
## $ elev : num 7.91 6.98 7.8 7.66 7.48 ...
## $ dist : num 0.00136 0.01222 0.10303 0.19009 0.27709 ...
## $ om : num 13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
## $ ffreq : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ soil : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
## $ lime : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
## $ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
## $ dist.m : num 50 30 150 270 380 470 240 120 240 420 ...meuse_sf <- st_as_sf(
meuse,
coords = c("x", "y"),
crs = 28992
)
grid_sf <- st_as_sf(
meuse.grid,
coords = c("x", "y"),
crs = 28992
)
st_crs(meuse_sf)## Coordinate Reference System:
## User input: EPSG:28992
## wkt:
## PROJCRS["Amersfoort / RD New",
## BASEGEOGCRS["Amersfoort",
## DATUM["Amersfoort",
## ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4289]],
## CONVERSION["RD New",
## METHOD["Oblique Stereographic",
## ID["EPSG",9809]],
## PARAMETER["Latitude of natural origin",52.1561605555556,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",5.38763888888889,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9999079,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",155000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",463000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["easting (X)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing (Y)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["Netherlands - onshore, including Waddenzee, Dutch Wadden Islands and 12-mile offshore coastal zone."],
## BBOX[50.75,3.2,53.7,7.22]],
## ID["EPSG",28992]]ggplot() +
geom_sf(
data = meuse_sf,
aes(color = zinc),
size = 2
) +
scale_color_viridis_c() +
labs(
title = "Sebaran Konsentrasi Zinc",
color = "Zinc"
) +
theme_minimal()
Closest Observation
closest_obs <- function(train, newdata) {
nearest <- st_nearest_feature(
newdata,
train
)
pred <- train$zinc[nearest]
return(pred)
}Nearest Neighbors (NN)
nn_interp <- function(train, newdata) {
model_nn <- gstat(
formula = zinc ~ 1,
locations = train,
nmax = 5,
set = list(idp = 0)
)
pred <- predict(
model_nn,
newdata
)
return(pred$var1.pred)
}Inverse Distance Weighting (IDW)
idw_interp <- function(train, newdata) {
model_idw <- gstat(
formula = zinc ~ 1,
locations = train,
nmax = nrow(train),
set = list(idp = 1)
)
pred <- predict(
model_idw,
newdata
)
return(pred$var1.pred)
}Ensemble
# 1. Closest Observation
pred_closest <- closest_obs(
meuse_sf,
grid_sf
)
# 2. Nearest Neighbors
pred_nn <- nn_interp(
meuse_sf,
grid_sf
)## [inverse distance weighted interpolation]# 3. IDW
pred_idw <- idw_interp(
meuse_sf,
grid_sf
)## [inverse distance weighted interpolation]pred_ensemble <- (
pred_closest +
pred_nn +
pred_idw
) / 3Kriging
kriging_interp <- function(train, newdata) {
v <- variogram(
log(zinc) ~ 1,
data = train
)
initial_model <- vgm(
psill = var(log(train$zinc)) * 0.8,
model = "Sph",
range = 900,
nugget = var(log(train$zinc)) * 0.2
)
fitted_model <- tryCatch(
fit.variogram(v, initial_model),
error = function(e) initial_model
)
pred <- krige(
log(zinc) ~ 1,
locations = train,
newdata = newdata,
model = fitted_model
)
return(exp(pred$var1.pred))
}
pred_kriging <- kriging_interp(
meuse_sf,
grid_sf
)## [using ordinary kriging]Hasil Prediksi
hasil_prediksi <- grid_sf %>%
mutate(
Closest = pred_closest,
NN = pred_nn,
IDW = pred_idw,
Ensemble = pred_ensemble,
Kriging = pred_kriging
)
# Melihat hasil prediksi
head(hasil_prediksi)## Simple feature collection with 6 features and 10 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 181100 ymin: 333660 xmax: 181220 ymax: 333740
## Projected CRS: Amersfoort / RD New
## part.a part.b dist soil ffreq geometry Closest NN
## 1 1 0 0.0000000 1 1 POINT (181180 333740) 1022 681.2
## 2 1 0 0.0000000 1 1 POINT (181140 333700) 1022 681.2
## 3 1 0 0.0122243 1 1 POINT (181180 333700) 1022 681.2
## 4 1 0 0.0434678 1 1 POINT (181220 333700) 640 681.2
## 5 1 0 0.0000000 1 1 POINT (181100 333660) 1022 681.2
## 6 1 0 0.0122243 1 1 POINT (181140 333660) 1022 681.2
## IDW Ensemble Kriging
## 1 512.6391 738.6130 664.8873
## 2 531.5849 744.9283 751.7104
## 3 518.2480 740.4827 668.5832
## 4 506.7441 609.3147 594.4192
## 5 573.3646 758.8549 866.5249
## 6 544.8959 749.3653 761.6678Visualisasi hasil interpolasi
koordinat <- st_coordinates(hasil_prediksi)
peta_data <- cbind(
st_drop_geometry(hasil_prediksi),
x = koordinat[, 1],
y = koordinat[, 2]
)
buat_peta <- function(data, variabel, judul) {
ggplot() +
geom_raster(
data = data,
aes(
x = x,
y = y,
fill = .data[[variabel]]
)
) +
geom_point(
data = meuse,
aes(x = x, y = y),
color = "black",
size = 0.8,
alpha = 0.7
) +
scale_fill_viridis_c(
option = "viridis",
name = "Zinc\n(mg/kg)",
guide = guide_colorbar(
barheight = unit(5, "cm"),
barwidth = unit(0.5, "cm")
)
) +
labs(
title = judul,
x = "Koordinat X",
y = "Koordinat Y"
) +
coord_equal() +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(
size = 14,
face = "bold",
hjust = 0.5
),
axis.title = element_text(size = 11),
axis.text = element_text(size = 9),
legend.title = element_text(size = 10),
legend.text = element_text(size = 9),
panel.grid = element_line(
color = "grey90",
linewidth = 0.3
),
plot.margin = margin(10, 15, 10, 10)
)
}
# Membuat peta setiap metode
p1 <- buat_peta(
peta_data,
"Closest",
"Closest Observation"
)
p2 <- buat_peta(
peta_data,
"NN",
"Nearest Neighbors"
)
p3 <- buat_peta(
peta_data,
"IDW",
"Inverse Distance Weighting"
)
p4 <- buat_peta(
peta_data,
"Ensemble",
"Ensemble Interpolation"
)
p5 <- buat_peta(
peta_data,
"Kriging",
"Ordinary Kriging"
)
# Menampilkan satu per satu
print(p1)
print(p2)
print(p3)
print(p4)
print(p5)
Cross Validation dan perhitungan RMSE
RMSE <- function(actual, predicted) {
sqrt(mean((actual - predicted)^2, na.rm = TRUE))
}
set.seed(123)
K <- 5
fold_id <- sample(
rep(1:K, length.out = nrow(meuse_sf))
)
hasil_cv <- data.frame(
Fold = 1:K,
Closest = NA,
NN = NA,
IDW = NA,
Ensemble = NA,
Kriging = NA
)
pred_cv <- data.frame(
Aktual = meuse$zinc,
Closest = NA,
NN = NA,
IDW = NA,
Ensemble = NA,
Kriging = NA
)
for (k in 1:K) {
train <- meuse_sf[fold_id != k, ]
test <- meuse_sf[fold_id == k, ]
aktual <- test$zinc
p1 <- closest_obs(train, test)
p2 <- nn_interp(train, test)
p3 <- idw_interp(train, test)
p4 <- (p1 + p2 + p3) / 3
p5 <- kriging_interp(train, test)
hasil_cv$Closest[k] <- RMSE(aktual, p1)
hasil_cv$NN[k] <- RMSE(aktual, p2)
hasil_cv$IDW[k] <- RMSE(aktual, p3)
hasil_cv$Ensemble[k] <- RMSE(aktual, p4)
hasil_cv$Kriging[k] <- RMSE(aktual, p5)
idx <- which(fold_id == k)
pred_cv$Closest[idx] <- p1
pred_cv$NN[idx] <- p2
pred_cv$IDW[idx] <- p3
pred_cv$Ensemble[idx] <- p4
pred_cv$Kriging[idx] <- 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]print(hasil_cv)## Fold Closest NN IDW Ensemble Kriging
## 1 1 353.6444 329.2993 394.6411 324.9896 320.8955
## 2 2 239.4902 217.0730 310.8441 215.4481 183.4711
## 3 3 336.0191 250.4189 295.4537 262.4020 184.8995
## 4 4 398.2811 284.3271 326.3475 294.8776 215.7763
## 5 5 406.6861 244.0397 317.5324 277.8594 216.9282rata_rmse <- data.frame(
Metode = c(
"Closest Observation",
"Nearest Neighbors",
"IDW",
"Ensemble",
"Kriging"
),
RMSE = c(
mean(hasil_cv$Closest),
mean(hasil_cv$NN),
mean(hasil_cv$IDW),
mean(hasil_cv$Ensemble),
mean(hasil_cv$Kriging)
)
)
rata_rmse <- rata_rmse %>%
arrange(RMSE)
print(rata_rmse)## Metode RMSE
## 1 Kriging 224.3941
## 2 Nearest Neighbors 265.0316
## 3 Ensemble 275.1153
## 4 IDW 328.9638
## 5 Closest Observation 346.8242Visualisasi Perbandingan RMSE
ggplot(
rata_rmse,
aes(
x = reorder(Metode, RMSE),
y = RMSE
)
) +
geom_col(
fill = "steelblue",
width = 0.7
) +
geom_text(
aes(label = round(RMSE, 2)),
hjust = -0.2,
size = 4
) +
coord_flip() +
labs(
title = "Perbandingan RMSE Metode Interpolasi",
subtitle = "Evaluasi menggunakan 5-Fold Cross Validation",
x = "Metode Interpolasi",
y = "RMSE"
) +
expand_limits(
y = max(rata_rmse$RMSE) * 1.15
) +
theme_minimal()
Menentukan metode terbaik
Berdasarkan hasil evaluasi menggunakan 5-Fold Cross Validation, diperoleh nilai RMSE yang berbeda pada setiap metode interpolasi. Metode Ordinary Kriging menghasilkan nilai RMSE terkecil, yaitu 224,39, sehingga menjadi metode interpolasi terbaik dibandingkan metode lainnya. Selanjutnya, Nearest Neighbors menghasilkan RMSE sebesar 265,03, Ensemble sebesar 275,12, IDW sebesar 328,96, dan Closest Observation sebesar 346,82. Nilai RMSE yang lebih kecil menunjukkan tingkat kesalahan prediksi yang lebih rendah. Dengan demikian, Ordinary Kriging memiliki kemampuan prediksi konsentrasi zinc yang paling baik pada dataset Meuse dibandingkan keempat metode lainnya.