Siapkan library dan data
library(sp)
## Warning: package 'sp' was built under R version 4.5.3
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(gstat)
## Warning: package 'gstat' 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:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
## Warning: package 'ggplot2' 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
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
data(meuse)
data(meuse.grid)
df <- meuse
grid <- meuse.grid
head(grid)
str(df)
## '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 ...
Data Prep
df <- sf::st_as_sf(df, coords = c("x", "y"), crs = 28992, remove = FALSE)
grid <- sf::st_as_sf(grid, coords = c("x", "y"), crs = 28992, remove = FALSE)
sf::st_geometry(df)
## Geometry set for 155 features
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 178605 ymin: 329714 xmax: 181390 ymax: 333611
## Projected CRS: Amersfoort / RD New
## First 5 geometries:
## POINT (181072 333611)
## POINT (181025 333558)
## POINT (181165 333537)
## POINT (181298 333484)
## POINT (181307 333330)
sf::st_geometry(grid)
## Geometry set for 3103 features
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 178460 ymin: 329620 xmax: 181540 ymax: 333740
## Projected CRS: Amersfoort / RD New
## First 5 geometries:
## POINT (181180 333740)
## POINT (181140 333700)
## POINT (181180 333700)
## POINT (181220 333700)
## POINT (181100 333660)
AED
ggplot(df) +
geom_sf(aes(color = zinc), size = 2.5) +
scale_color_viridis_c() +
labs(
title = "Persebaran Konsentrasi Zinc",
color = "Zinc"
) +
theme_minimal()
ggplot(grid) +
geom_sf(color = "blue", size = 0.5) +
labs(
title = "Lokasi Grid Prediksi"
) +
theme_minimal()
Metode interpolasi
RMSE <- function(observed, predicted) {sqrt(mean((observed - predicted)^2))}
set.seed(123)
k <- 4
kf <- kfold(nrow(df), k = k)
rmse1 <- rep(NA, k)
rmse2 <- rep(NA, k)
rmse3 <- rep(NA, k)
rmse4 <- rep(NA, k)
for (i in 1:k) {
test <- df[kf == i, ]
train <- df[kf != i, ]
# 1. Closest Observation
nearest <- st_nearest_feature(test, train)
p1 <- train$zinc[nearest]
rmse1[i] <- RMSE(test$zinc, p1)
# 2. INVERSE DISTANCE WEIGHTING (IDW)
idw_model <- gstat(
formula = zinc ~ 1,
locations = train,
nmax = nrow(train),
set = list(idp = 1)
)
p2 <- predict(idw_model, test)$var1.pred
rmse2[i] <- RMSE(test$zinc, p2)
# 3. NEAREST NEIGHBORS (NN)
nn_model <- gstat(
formula = zinc ~ 1,
locations = train,
nmax = 5,
set = list(idp = 0)
)
p3 <- predict(nn_model, test)$var1.pred
rmse3[i] <- RMSE(test$zinc, p3)
# 4. ENSEMBLE
p4 <- (p1 + p2 + p3) / 3
rmse4[i] <- RMSE(test$zinc, p4)
}
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
Evaluasi RMSE
hasil_cv <- data.frame(
Fold = 1:k,
Closest_Observation = rmse1,
IDW = rmse2,
Nearest_Neighbors = rmse3,
Ensemble = rmse4
)
print(hasil_cv)
## Fold Closest_Observation IDW Nearest_Neighbors Ensemble
## 1 1 393.7434 282.7544 232.2025 248.7509
## 2 2 364.5879 316.9761 290.8121 274.1471
## 3 3 271.5254 417.0272 333.2966 317.7759
## 4 4 287.9924 304.7345 280.6150 266.7027
hasil_akhir <- data.frame(
Metode = c(
"Closest Observation",
"IDW",
"Nearest Neighbors",
"Ensemble"
),
Mean_RMSE = c(
mean(rmse1),
mean(rmse2),
mean(rmse3),
mean(rmse4)
)
)
hasil_akhir <- hasil_akhir %>%
arrange(Mean_RMSE)
print(hasil_akhir)
## Metode Mean_RMSE
## 1 Ensemble 276.8441
## 2 Nearest Neighbors 284.2316
## 3 Closest Observation 329.4623
## 4 IDW 330.3730
Interpolasi
# A. CLOSEST OBSERVATION
nearest_full <- st_nearest_feature(grid, df)
grid$pred_closest <- df$zinc[nearest_full]
# B. IDW
idw_full <- gstat(formula = zinc ~ 1, locations = df, nmax = nrow(df), set = list(idp = 1))
grid$pred_idw <- predict(idw_full, grid)$var1.pred
## [inverse distance weighted interpolation]
# C. NEAREST NEIGHBORS
nn_full <- gstat(formula = zinc ~ 1, locations = df,nmax = 5, set = list(idp = 0))
grid$pred_nn <- predict(nn_full, grid)$var1.pred
## [inverse distance weighted interpolation]
# D. ENSEMBLE
grid$pred_ensemble <- (grid$pred_closest +grid$pred_idw +grid$pred_nn) / 3
Visualisasi
plot_interpolasi <- function(data, variabel, judul) {
ggplot() +
geom_sf(
data = data,
aes(color = .data[[variabel]]),
size = 1
) +
geom_sf(
data = df,
color = "black",
size = 0.5
) +
scale_color_viridis_c() +
labs(
title = judul,
color = "Prediksi Zinc"
) +
theme_minimal()
}
p1 <- plot_interpolasi(
grid,
"pred_closest",
"Closest Observation"
)
p2 <- plot_interpolasi(
grid,
"pred_idw",
"Inverse Distance Weighting"
)
p3 <- plot_interpolasi(
grid,
"pred_nn",
"Nearest Neighbors"
)
p4 <- plot_interpolasi(
grid,
"pred_ensemble",
"Ensemble"
)
print(p1)
print(p2)
print(p3)
print(p4)
