ANALISIS INTERPOLASI SPASIAL PADA DATA MEUSE
Tri Ananda Marpaung|M0401241098
10 November 2026
## Warning: package 'sp' was built under R version 4.5.3## Warning: package 'gstat' was built under R version 4.5.3## 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## 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, union## Warning: package 'ggplot2' was built under R version 4.5.3## Warning: package 'FNN' was built under R version 4.5.3## Warning: package 'terra' was built under R version 4.5.3## terra 1.9.34## Warning: package 'whitebox' was built under R version 4.5.31 Memanggil Data
## 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 470## '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 ...## x y cadmium copper
## Min. :178605 Min. :329714 Min. : 0.200 Min. : 14.00
## 1st Qu.:179371 1st Qu.:330762 1st Qu.: 0.800 1st Qu.: 23.00
## Median :179991 Median :331633 Median : 2.100 Median : 31.00
## Mean :180005 Mean :331635 Mean : 3.246 Mean : 40.32
## 3rd Qu.:180630 3rd Qu.:332463 3rd Qu.: 3.850 3rd Qu.: 49.50
## Max. :181390 Max. :333611 Max. :18.100 Max. :128.00
##
## lead zinc elev dist
## Min. : 37.0 Min. : 113.0 Min. : 5.180 Min. :0.00000
## 1st Qu.: 72.5 1st Qu.: 198.0 1st Qu.: 7.546 1st Qu.:0.07569
## Median :123.0 Median : 326.0 Median : 8.180 Median :0.21184
## Mean :153.4 Mean : 469.7 Mean : 8.165 Mean :0.24002
## 3rd Qu.:207.0 3rd Qu.: 674.5 3rd Qu.: 8.955 3rd Qu.:0.36407
## Max. :654.0 Max. :1839.0 Max. :10.520 Max. :0.88039
##
## om ffreq soil lime landuse dist.m
## Min. : 1.000 1:84 1:97 0:111 W :50 Min. : 10.0
## 1st Qu.: 5.300 2:48 2:46 1: 44 Ah :39 1st Qu.: 80.0
## Median : 6.900 3:23 3:12 Am :22 Median : 270.0
## Mean : 7.478 Fw :10 Mean : 290.3
## 3rd Qu.: 9.000 Ab : 8 3rd Qu.: 450.0
## Max. :17.000 (Other):25 Max. :1000.0
## NA's :2 NA's : 1# Memilih variabel yang digunakan
data_meuse <- meuse %>%
select(x, y, zinc) %>%
na.omit()
# Jumlah observasi
nrow(data_meuse)## [1] 155## x y zinc
## 0 0 0# Visualisasi persebaran zinc
ggplot(data_meuse, aes(x = x, y = y, color = zinc)) +
geom_point(size = 3) +
scale_color_viridis_c() +
labs(
title = "Persebaran Konsentrasi Zinc pada Data Meuse",
x = "Koordinat X",
y = "Koordinat Y",
color = "Zinc"
) +
theme_minimal()
# 1. Closest Observation
pred_closest <- function(train, test) {
koordinat_train <- as.matrix(
train[, c("x", "y")]
)
koordinat_test <- as.matrix(
test[, c("x", "y")]
)
indeks <- FNN::get.knnx(
data = koordinat_train,
query = koordinat_test,
k = 1
)$nn.index[, 1]
return(train$zinc[indeks])
}# Membuat grid interpolasi
grid_interpolasi <- expand.grid(
x = seq(
min(data_meuse$x),
max(data_meuse$x),
length.out = 100
),
y = seq(
min(data_meuse$y),
max(data_meuse$y),
length.out = 100
)
)
# Melihat jumlah titik grid
dim(grid_interpolasi)## [1] 10000 2## x y
## 1 178605.0 329714
## 2 178633.1 329714
## 3 178661.3 329714
## 4 178689.4 329714
## 5 178717.5 329714
## 6 178745.7 329714# Fungsi interpolasi Closest Observation
pred_closest <- function(train, test) {
koordinat_train <- as.matrix(
train[, c("x", "y")]
)
koordinat_test <- as.matrix(
test[, c("x", "y")]
)
indeks <- FNN::get.knnx(
data = koordinat_train,
query = koordinat_test,
k = 1
)$nn.index[, 1]
prediksi <- train$zinc[indeks]
return(prediksi)
}# Melakukan interpolasi
hasil_closest <- pred_closest(
train = data_meuse,
test = grid_interpolasi
)
# Menggabungkan hasil prediksi dengan koordinat grid
grid_closest <- grid_interpolasi %>%
mutate(zinc_pred = hasil_closest)
# Melihat hasil interpolasi
head(grid_closest)## x y zinc_pred
## 1 178605.0 329714 783
## 2 178633.1 329714 783
## 3 178661.3 329714 783
## 4 178689.4 329714 783
## 5 178717.5 329714 783
## 6 178745.7 329714 783## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 113.0 199.0 403.0 600.7 933.0 1839.0ggplot() +
geom_raster(
data = grid_closest,
aes(
x = x,
y = y,
fill = zinc_pred
)
) +
geom_point(
data = data_meuse,
aes(x = x, y = y),
color = "black",
size = 0.7
) +
scale_fill_viridis_c(
option = "plasma",
name = "Prediksi Zinc"
) +
coord_equal() +
labs(
title = "Interpolasi Closest Observation",
subtitle = "Perkiraan konsentrasi zinc berdasarkan titik terdekat",
x = "Koordinat X",
y = "Koordinat Y"
) +
theme_minimal() +
theme(
axis.text.x = element_text(
angle = 45,
hjust = 1,
size = 8
),
axis.text.y = element_text(size = 8),
plot.title = element_text(
face = "bold",
size = 13
)
)
Berdasarkan hasil interpolasi menggunakan metode Closest Observation, diperoleh pola persebaran konsentrasi zinc yang membentuk wilayah-wilayah dengan batas yang relatif tegas. Metode ini menentukan nilai prediksi berdasarkan titik pengamatan yang memiliki jarak paling dekat dengan lokasi yang akan diestimasi. Setiap lokasi dalam wilayah tertentu akan memperoleh nilai zinc yang sama dari titik pengamatan terdekatnya.
Perbedaan warna pada peta menunjukkan variasi konsentrasi zinc secara spasial. Metode ini mampu menggambarkan persebaran nilai berdasarkan kedekatan lokasi, tetapi menghasilkan perubahan nilai yang tidak kontinu antarwilayah karena hanya menggunakan satu titik pengamatan terdekat. Titik
2 2. NATURAL NEIGHBOR
pred_nn <- function(train, test) {
folder <- tempfile("nn_")
dir.create(folder)
input_shp <- file.path(folder, "training.shp")
output_tif <- file.path(folder, "hasil_nn.tif")
# Membuat data spasial training
train_sf <- st_as_sf(
train,
coords = c("x", "y"),
crs = 28992,
remove = FALSE
)
st_write(
train_sf,
input_shp,
quiet = TRUE
)
# Interpolasi Natural Neighbor
whitebox::wbt_natural_neighbour_interpolation(
input = input_shp,
output = output_tif,
field = "zinc",
use_z = FALSE,
cell_size = 5,
clip = FALSE,
verbose_mode = FALSE
)
# Membaca raster
raster_nn <- terra::rast(output_tif)
# Menggunakan CRS raster yang dihasilkan
# Koordinat x-y tetap menggunakan angka aslinya
titik_uji <- terra::vect(
data.frame(
x = test$x,
y = test$y
),
geom = c("x", "y"),
crs = terra::crs(raster_nn)
)
# Ekstraksi prediksi
hasil <- terra::extract(
raster_nn,
titik_uji
)
prediksi <- hasil[[2]]
unlink(folder, recursive = TRUE)
return(prediksi)
}# Mengubah data menjadi objek sf
meuse_sf <- st_as_sf(
data_meuse,
coords = c("x", "y"),
crs = 28992,
remove = FALSE
)
# Melihat data spasial
head(meuse_sf)## Simple feature collection with 6 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 181025 ymin: 333260 xmax: 181390 ymax: 333611
## Projected CRS: Amersfoort / RD New
## x y zinc geometry
## 1 181072 333611 1022 POINT (181072 333611)
## 2 181025 333558 1141 POINT (181025 333558)
## 3 181165 333537 640 POINT (181165 333537)
## 4 181298 333484 257 POINT (181298 333484)
## 5 181307 333330 269 POINT (181307 333330)
## 6 181390 333260 281 POINT (181390 333260)## 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]]# Membuat folder sementara
folder_nn <- tempfile("natural_neighbor_")
dir.create(folder_nn)
# Menentukan lokasi file
input_nn <- file.path(folder_nn, "data_meuse.shp")
output_nn <- file.path(folder_nn, "hasil_nn.tif")
# Menyimpan data spasial sebagai shapefile
st_write(
meuse_sf,
input_nn,
quiet = TRUE
)
# Melakukan interpolasi Natural Neighbor
whitebox::wbt_natural_neighbour_interpolation(
input = input_nn,
output = output_nn,
field = "zinc",
use_z = FALSE,
cell_size = 5,
clip = FALSE,
verbose_mode = FALSE
)# Membaca hasil interpolasi
raster_nn <- terra::rast(output_nn)
# Melihat informasi raster
raster_nn## class : SpatRaster
## size : 780, 557, 1 (nrow, ncol, nlyr)
## resolution : 5, 5 (x, y)
## extent : 178605, 181390, 329711, 333611 (xmin, xmax, ymin, ymax)
## coord. ref. : unnamed
## source : hasil_nn.tif
## name : hasil_nn# Mengubah raster menjadi data frame
hasil_nn <- as.data.frame(
raster_nn,
xy = TRUE,
na.rm = TRUE
)
# Mengganti nama kolom
names(hasil_nn) <- c(
"x",
"y",
"zinc_pred"
)
# Melihat hasil
head(hasil_nn)## x y zinc_pred
## 61 178907.5 333608.5 711.0000
## 62 178912.5 333608.5 716.4384
## 63 178917.5 333608.5 720.7386
## 64 178922.5 333608.5 723.1122
## 65 178927.5 333608.5 724.5789
## 66 178932.5 333608.5 725.5641## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 113.1 211.8 473.2 601.7 902.9 1839.0ggplot() +
geom_raster(
data = hasil_nn,
aes(
x = x,
y = y,
fill = zinc_pred
)
) +
geom_point(
data = data_meuse,
aes(x = x, y = y),
color = "black",
size = 0.7
) +
scale_fill_viridis_c(
option = "plasma",
name = "Prediksi Zinc"
) +
coord_equal() +
labs(
title = "Interpolasi Natural Neighbor",
subtitle = "Perkiraan konsentrasi zinc berdasarkan titik tetangga alami",
x = "Koordinat X",
y = "Koordinat Y"
) +
theme_minimal() +
theme(
axis.text.x = element_text(
angle = 45,
hjust = 1,
size = 8
),
axis.text.y = element_text(size = 8),
plot.title = element_text(
face = "bold",
size = 13
)
)
Berdasarkan hasil interpolasi menggunakan metode Natural Neighbor,
diperoleh gambaran persebaran konsentrasi zinc yang lebih halus
dibandingkan metode Closest Observation. Metode Natural Neighbor
memperkirakan nilai pada lokasi yang tidak teramati dengan memanfaatkan
titik-titik pengamatan di sekitarnya yang memiliki hubungan spasial
berdasarkan pembagian wilayah Voronoi.
Nilai prediksi pada suatu lokasi diperoleh melalui pembobotan berdasarkan luas wilayah tetangga alami yang terbentuk setelah titik interpolasi dimasukkan ke dalam jaringan Voronoi. Dengan pendekatan tersebut, metode ini mempertimbangkan kontribusi beberapa titik pengamatan di sekitar lokasi yang diestimasi, bukan hanya satu titik terdekat. # 3. IDW
pred_idw <- function(train, test, power = 2) {
train_sp <- train
coordinates(train_sp) <- ~x+y
test_sp <- test
coordinates(test_sp) <- ~x+y
hasil <- gstat::idw(
zinc ~ 1,
locations = train_sp,
newdata = test_sp,
idp = power
)
return(hasil$var1.pred)
}# Menyiapkan data pengamatan
data_idw <- data_meuse %>%
select(x, y, zinc) %>%
filter(complete.cases(.)) %>%
distinct(x, y, .keep_all = TRUE)
# Memastikan data tersedia
head(data_idw)## x y zinc
## 1 181072 333611 1022
## 2 181025 333558 1141
## 3 181165 333537 640
## 4 181298 333484 257
## 5 181307 333330 269
## 6 181390 333260 281## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 113.0 198.0 326.0 469.7 674.5 1839.0# Membuat grid interpolasi IDW
grid_idw <- expand.grid(
x = seq(
min(data_idw$x),
max(data_idw$x),
length.out = 100
),
y = seq(
min(data_idw$y),
max(data_idw$y),
length.out = 100
)
)
# Menghitung prediksi IDW pada grid
prediksi_grid_idw <- pred_idw(
train = data_idw,
test = grid_idw,
power = 2
)## [inverse distance weighted interpolation]## x y zinc_pred
## 1 178605.0 329714 513.8386
## 2 178633.1 329714 525.7315
## 3 178661.3 329714 539.9492
## 4 178689.4 329714 555.9933
## 5 178717.5 329714 572.2124
## 6 178745.7 329714 585.2864## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 121.4 372.4 455.1 480.9 565.4 1681.8data_idw <- data_meuse %>%
select(x, y, zinc) %>%
filter(complete.cases(.)) %>%
distinct(x, y, .keep_all = TRUE)
head(data_idw)## x y zinc
## 1 181072 333611 1022
## 2 181025 333558 1141
## 3 181165 333537 640
## 4 181298 333484 257
## 5 181307 333330 269
## 6 181390 333260 281## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 113.0 198.0 326.0 469.7 674.5 1839.0## [inverse distance weighted interpolation]## [1] 793.8598## num 794ggplot() +
geom_raster(
data = grid_idw,
aes(
x = x,
y = y,
fill = zinc_pred
)
) +
geom_point(
data = data_idw,
aes(
x = x,
y = y
),
color = "black",
size = 0.8
) +
scale_fill_viridis_c(
option = "plasma",
name = "Prediksi Zinc"
) +
coord_equal() +
labs(
title = "Peta Interpolasi Inverse Distance Weighting",
subtitle = "IDW dengan parameter power = 2",
x = "Koordinat X",
y = "Koordinat Y"
) +
theme_minimal()
## [1] 10000## [1] 0## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 121.4 372.4 455.1 480.9 565.4 1681.8## [1] 121.4168 1681.7618Berdasarkan hasil interpolasi menggunakan metode Inverse Distance Weighting (IDW) dengan parameter power sebesar 2, diperoleh pola persebaran konsentrasi zinc yang dipengaruhi oleh jarak antara lokasi prediksi dengan titik-titik pengamatan. Metode ini memberikan bobot yang lebih besar kepada titik pengamatan yang memiliki jarak lebih dekat, sedangkan titik yang lebih jauh memperoleh bobot yang lebih kecil.
Penggunaan parameter power = 2 menunjukkan bahwa pengaruh jarak terhadap bobot bersifat kuadratik. Artinya, semakin jauh suatu titik pengamatan dari lokasi yang diprediksi, semakin kecil kontribusinya terhadap nilai interpolasi. Pola hasil interpolasi memperlihatkan perubahan konsentrasi zinc berdasarkan pengaruh nilai-nilai pengamatan di sekitarnya.
3 4. ORDINARY KRINGING’
pred_kriging <- function(train, test) {
train_sp <- train
coordinates(train_sp) <- ~x+y
test_sp <- test
coordinates(test_sp) <- ~x+y
# Variogram eksperimental
vg <- gstat::variogram(
zinc ~ 1,
train_sp
)
# Model variogram awal
model_awal <- gstat::vgm(
psill = var(train$zinc) * 0.8,
model = "Exp",
range = 500,
nugget = var(train$zinc) * 0.2
)
# Estimasi parameter variogram
model_fit <- gstat::fit.variogram(
vg,
model = model_awal
)
# Ordinary Kriging
hasil <- gstat::krige(
zinc ~ 1,
locations = train_sp,
newdata = test_sp,
model = model_fit
)
return(hasil$var1.pred)
}# Menyiapkan data pengamatan
data_ok <- data_meuse %>%
select(x, y, zinc) %>%
filter(complete.cases(.)) %>%
distinct(x, y, .keep_all = TRUE)
# Mengubah menjadi objek spasial
data_ok_sf <- st_as_sf(
data_ok,
coords = c("x", "y"),
crs = 28992
)
# Mengubah menjadi Spatial untuk gstat
data_ok_sp <- as(data_ok_sf, "Spatial")
# Memeriksa data
summary(data_ok$zinc)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 113.0 198.0 326.0 469.7 674.5 1839.0# Membuat semivariogram empiris
variogram_emp <- variogram(
zinc ~ 1,
data = data_ok_sp
)
# Melihat hasil semivariogram
print(variogram_emp)## np dist gamma dir.hor dir.ver id
## 1 57 79.29244 37362.96 0 0 var1
## 2 299 163.97367 72718.34 0 0 var1
## 3 419 267.36483 82655.53 0 0 var1
## 4 457 372.73542 111575.91 0 0 var1
## 5 547 478.47670 123080.69 0 0 var1
## 6 533 585.34058 147414.28 0 0 var1
## 7 574 693.14526 142891.51 0 0 var1
## 8 564 796.18365 153563.90 0 0 var1
## 9 589 903.14650 160217.69 0 0 var1
## 10 543 1011.29177 168867.36 0 0 var1
## 11 500 1117.86235 186528.20 0 0 var1
## 12 477 1221.32810 150320.76 0 0 var1
## 13 452 1329.16407 180399.66 0 0 var1
## 14 457 1437.25620 143139.84 0 0 var1
## 15 415 1543.20248 144112.31 0 0 var1# Visualisasi semivariogram empiris
plot(
variogram_emp,
main = "Experimental Semivariogram Zinc",
xlab = "Jarak (h)",
ylab = "Semivarian"
)
# Menghitung jarak maksimum antar titik
jarak_maks <- sqrt(
(max(data_ok$x) - min(data_ok$x))^2 +
(max(data_ok$y) - min(data_ok$y))^2
)
# Varians data zinc
var_zinc <- var(data_ok$zinc)
# Membuat model awal
model_sph <- vgm(
psill = 0.8 * var_zinc,
model = "Sph",
range = jarak_maks / 3,
nugget = 0.2 * var_zinc
)
model_exp <- vgm(
psill = 0.8 * var_zinc,
model = "Exp",
range = jarak_maks / 3,
nugget = 0.2 * var_zinc
)
model_gau <- vgm(
psill = 0.8 * var_zinc,
model = "Gau",
range = jarak_maks / 3,
nugget = 0.2 * var_zinc
)
# Fitting model semivariogram
fit_sph <- fit.variogram(
variogram_emp,
model_sph
)## Warning in fit.variogram(variogram_emp, model_sph): singular model in variogram
## fit## Warning in fit.variogram(variogram_emp, model_exp): No convergence after 200
## iterations: try different initial values?## Warning in fit.variogram(variogram_emp, model_gau): No convergence after 200
## iterations: try different initial values?## model psill range
## 1 Nug 30596.76 0.0000
## 2 Sph 80435.68 -475.8606## model psill range
## 1 Nug 15209.396 0.0000
## 2 Exp 6766.533 -781.4034## model psill range
## 1 Nug 27278.14 0.000
## 2 Gau 19648.48 -5378.415# Membuat tabel perbandingan
perbandingan_variogram <- data.frame(
Model = c("Spherical", "Exponential", "Gaussian"),
SSE = c(
attr(fit_sph, "SSErr"),
attr(fit_exp, "SSErr"),
attr(fit_gau, "SSErr")
)
)
print(perbandingan_variogram)## Model SSE
## 1 Spherical 83422683
## 2 Exponential 287526841
## 3 Gaussian 188456976
## model psill range
## 1 Nug 30596.76 0.0000
## 2 Sph 80435.68 -475.8606
## model psill range
## 1 Nug 15209.396 0.0000
## 2 Exp 6766.533 -781.4034
## model psill range
## 1 Nug 27278.14 0.000
## 2 Gau 19648.48 -5378.415# Menggabungkan model semivariogram
model_list <- list(
Spherical = fit_sph,
Exponential = fit_exp,
Gaussian = fit_gau
)
# Mengambil nilai SSE masing-masing model
sse_values <- sapply(
model_list,
function(x) attr(x, "SSErr")
)
# Melihat perbandingan SSE
print(sse_values)## Spherical Exponential Gaussian
## 83422683 287526841 188456976# Menentukan model dengan SSE terkecil
indeks_terbaik <- which.min(sse_values)
model_terbaik <- model_list[[indeks_terbaik]]
# Menampilkan model terpilih
print(model_terbaik)## model psill range
## 1 Nug 30596.76 0.0000
## 2 Sph 80435.68 -475.8606## [1] "Spherical"## model psill range
## 1 Nug 30596.76 0.0000
## 2 Sph 80435.68 -475.8606## model psill range
## 1 Nug 15209.396 0.0000
## 2 Exp 6766.533 -781.4034## model psill range
## 1 Nug 27278.14 0.000
## 2 Gau 19648.48 -5378.415## [1] 83422683## [1] 287526841## [1] 188456976# Menghitung varians zinc
var_zinc <- var(data_ok$zinc)
# Menghitung jarak maksimum
jarak_maks <- sqrt(
(max(data_ok$x) - min(data_ok$x))^2 +
(max(data_ok$y) - min(data_ok$y))^2
)
# Membuat model awal
model_sph <- vgm(
psill = 0.8 * var_zinc,
model = "Sph",
range = jarak_maks / 3,
nugget = 0.2 * var_zinc
)
model_exp <- vgm(
psill = 0.8 * var_zinc,
model = "Exp",
range = jarak_maks / 3,
nugget = 0.2 * var_zinc
)
model_gau <- vgm(
psill = 0.8 * var_zinc,
model = "Gau",
range = jarak_maks / 3,
nugget = 0.2 * var_zinc
)
# Fitting ulang
fit_sph <- fit.variogram(
variogram_emp,
model_sph,
fit.method = 6,
fit.sills = TRUE,
fit.ranges = TRUE
)## Warning in fit.variogram(variogram_emp, model_sph, fit.method = 6, fit.sills =
## TRUE, : singular model in variogram fitfit_exp <- fit.variogram(
variogram_emp,
model_exp,
fit.method = 6,
fit.sills = TRUE,
fit.ranges = TRUE
)## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial
## values?## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fitfit_gau <- fit.variogram(
variogram_emp,
model_gau,
fit.method = 6,
fit.sills = TRUE,
fit.ranges = TRUE
)## Warning in fit.variogram(variogram_emp, model_gau, fit.method = 6, fit.sills =
## TRUE, : No convergence after 200 iterations: try different initial values?## model psill range
## 1 Nug 36191.27 0.00000
## 2 Sph 91749.49 31.74033## model psill range
## 1 Nug 0 0.000
## 2 Exp 0 1596.623## model psill range
## 1 Nug 28178.99 0.000
## 2 Gau 29205.73 -5755.216# Memeriksa SSE
sse_values <- c(
Spherical = attr(fit_sph, "SSErr"),
Exponential = attr(fit_exp, "SSErr"),
Gaussian = attr(fit_gau, "SSErr")
)
print(sse_values)## Spherical Exponential Gaussian
## 24749929312 292221155929 188350669632# Menggabungkan model semivariogram
model_list <- list(
Spherical = fit_sph,
Exponential = fit_exp,
Gaussian = fit_gau
)
# Fungsi pemeriksaan model
cek_model <- function(m) {
parameter <- as.data.frame(m)
# Mengubah kolom parameter menjadi numerik
psill <- as.numeric(parameter$psill)
range <- as.numeric(parameter$range)
# Memeriksa validitas parameter
all(
is.finite(psill),
is.finite(range),
psill >= 0,
range >= 0
) && any(range > 0) && any(psill > 0)
}
# Memeriksa semua model
valid_model <- sapply(
model_list,
cek_model
)
print(valid_model)## Spherical Exponential Gaussian
## TRUE FALSE FALSE## Spherical Exponential Gaussian
## 24749929312 292221155929 188350669632# Menghilangkan model yang tidak valid
sse_valid <- sse_values
sse_valid[!valid_model] <- Inf
sse_valid[!is.finite(sse_valid)] <- Inf
if (all(is.infinite(sse_valid))) {
stop("Tidak ada model semivariogram yang valid.")
}
# Memilih model terbaik
indeks_terbaik <- which.min(sse_valid)
model_terbaik <- model_list[[indeks_terbaik]]
cat("Model terbaik:",
names(model_list)[indeks_terbaik])## Model terbaik: Spherical## model psill range
## 1 Nug 36191.27 0.00000
## 2 Sph 91749.49 31.74033## 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 470## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 155 obs. of 1 variable:
## .. ..$ zinc: num [1:155] 1022 1141 640 257 269 ...
## ..@ coords.nrs : num(0)
## ..@ coords : num [1:155, 1:2] 181072 181025 181165 181298 181307 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:155] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
## ..@ bbox : num [1:2, 1:2] 178605 329714 181390 333611
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel"| __truncated__
## .. .. ..$ comment: chr "PROJCRS[\"Amersfoort / RD New\",\n BASEGEOGCRS[\"Amersfoort\",\n DATUM[\"Amersfoort\",\n E"| __truncated__## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"# Membuat grid berdasarkan batas koordinat data meuse
grid_ok <- expand.grid(
x = seq(
min(data_ok_sp@coords[, 1]),
max(data_ok_sp@coords[, 1]),
length.out = 100
),
y = seq(
min(data_ok_sp@coords[, 2]),
max(data_ok_sp@coords[, 2]),
length.out = 100
)
)
# Mengubah menjadi objek spasial
coordinates(grid_ok) <- ~ x + y
# Menyamakan sistem koordinat
proj4string(grid_ok) <- proj4string(data_ok_sp)hasil_ok <- gstat::krige(
zinc ~ 1,
locations = data_ok_sp,
newdata = grid_ok,
model = model_terbaik
)## [using ordinary kriging]# Mengubah hasil Ordinary Kriging menjadi data frame
hasil_ok_df <- as.data.frame(hasil_ok)
# Memeriksa hasil prediksi
head(hasil_ok_df)## x y var1.pred var1.var
## 1 178605.0 329714 469.7161 128766.2
## 2 178633.1 329714 469.7161 128766.2
## 3 178661.3 329714 469.7161 128766.2
## 4 178689.4 329714 469.7161 128766.2
## 5 178717.5 329714 469.7161 128766.2
## 6 178745.7 329714 469.7161 128766.2## [1] 469.7161 469.7223 488.7706 484.9643 538.4518 485.1317 469.7163 468.8025
## [9] 452.9216 540.1496 517.7870 467.0418 469.0769 438.3114 469.7280 470.5687
## [17] 469.7388 469.7028 447.7693 469.6682 427.9878 398.1753 438.7568 450.8026
## [25] 400.1983 467.1277 458.3066 464.4134 326.9332 437.5739 460.7260 409.7605
## [33] 374.2829 447.8479 376.6909 465.3773 417.7589 449.4791 451.6606 470.3235
## [41] 452.5612 454.5274 447.8234 431.5022 477.9463 469.8983 469.0985 395.4327
## [49] 467.2390 483.2746 965.1986 470.9469 491.8568 476.1317 444.6848 468.8446
## [57] 421.0363 449.8032 410.8558 361.7605 472.6874 518.0805 486.9712 469.0140
## [65] 469.7153 405.2405 467.6730 462.1612 469.1595 468.3272 467.0821 311.3479
## [73] 466.0943 465.5011 461.7830 498.4125 471.4967 483.4237 477.2268 469.5262
## [81] 468.6332 467.8607 469.6710 467.6219 469.4097 449.5990 293.4518 473.5527
## [89] 471.0073 467.3995 469.5003 357.7680 433.4911 419.1579 435.0567 472.2278
## [97] 495.2459 484.1776 502.4974 479.0233 490.0677 471.6200 477.6529 472.7255
## [105] 470.3695 417.5736 433.2676 482.0024 475.5761 532.9659 487.2906 484.1642
## [113] 468.9620 469.7140 469.6846 493.4124 873.8019 470.5469 488.5803 452.7304
## [121] 458.0550 421.6870 430.7211 480.3227 470.0704 457.6607 324.4624 465.8531
## [129] 449.0954 455.1426 524.5815 486.7331 446.1149 412.5567 470.3015 482.4467
## [137] 566.9207 510.7881 495.5013 514.0838 470.4680 475.9572 470.2412 469.2508
## [145] 467.6910 439.3576 458.9159 706.0489 576.8506 496.4050 409.6619 335.2810
## [153] 440.0979 428.1922 457.6844 467.4675 469.4329 506.5660 390.2232 425.5403
## [161] 528.0543 722.5178 413.0892 460.1076 449.4612 469.2893 524.1911 711.7095
## [169] 439.4259 456.5891 468.4473 464.0532 535.6229 553.1558 446.8185 460.8839
## [177] 469.4923 385.4061 445.0362 364.8107 374.9111 451.4493 454.2753 546.2758
## [185] 471.8369 467.3887 273.5658 464.0163 350.9671 441.2812 485.0117 468.4581
## [193] 417.7991 456.6224 414.4701 347.4799 481.6158 676.7233 470.1931 439.4714
## [201] 461.9342 468.2669 405.6334 469.4439 438.0861 438.9920 461.7228 454.0539
## [209] 458.7244 400.3069 388.0297 467.0296 482.5369 531.9844 433.8080 394.7457
## [217] 469.7294 468.2714 433.0941 372.5934 416.7910 442.8480 468.7101 461.7962
## [225] 447.9697 514.1737 498.3597 462.1218 469.6984 414.5193 367.9773 457.3486
## [233] 465.0456 354.7535 469.6204 419.6992 370.5411 425.6981 465.6297 438.8462
## [241] 465.0272 506.6922 552.0428 469.5463 462.3072 452.6388 369.5786 428.1053
## [249] 467.7598 463.5223 449.2966 432.8412 464.9052 458.3371 478.6126 534.4116
## [257] 345.9843 410.0447 442.7493 358.0695 466.4588 470.9852 502.5943 458.2859
## [265] 348.7511 469.1136 349.3586 468.4203 469.4497 468.9293 487.8369 490.4513
## [273] 449.8448 466.7331 469.6696 469.4230 368.2054 457.7581 533.4103 538.9505
## [281] 465.8984 469.3428 473.9353 469.7441 517.3302 473.0649 433.7550 451.5832
## [289] 589.8359 566.8964 533.4793 479.5030 475.1091 643.4898 471.6951 462.7313
## [297] 448.2813 507.9443 480.8821 622.9259 470.8847 459.1377 441.4890 394.7986
## [305] 380.4410 609.6970 541.8764 489.1363 548.1567 597.0380 469.7921 531.7111
## [313] 632.8511 513.4379 489.1281 475.9975 692.4646 506.8107 559.0408 469.7188
## [321] 543.3150 479.5113 441.9819 469.1353 468.8887 446.2193 420.5070 395.5856
## [329] 473.3837 709.0435 475.2919 453.4076 469.7130 455.9475 410.8224 410.4004
## [337] 494.2299 469.8820 441.2565 892.7076 537.3255 474.8463 439.1893 468.6463
## [345] 436.5346 469.7261 468.8334 437.9445 551.8562 524.1258 446.6020 535.9910
## [353] 575.2975 468.8689 424.3067 446.1777 424.5095 439.7744 469.5829 713.6661
## [361] 554.1084 538.0366 510.7947 461.9039 454.2652 428.0200 468.4904 469.5080
## [369] 373.5312 459.9518 397.4629 423.6695 469.4826 333.5423 463.0264 485.6810
## [377] 366.3857 422.9691 473.9202 521.7957 470.0261 469.9653 412.4018 438.0772
## [385] 469.9442 469.8943 469.3683 465.6191 459.9855 469.8387 452.4192 440.4162
## [393] 424.5428 363.1901 551.2698 488.4508 403.7696 453.3174 468.9561 485.7410
## [401] 525.3522 460.2072 468.6412 469.2085 406.7656 466.1339 469.7055 466.2190
## [409] 464.9902 452.0924 469.6523 458.0424 455.5657 396.8052 432.8386 474.8937
## [417] 441.3661 463.6482 394.8730 441.5667 670.9819 479.2415 488.7382 470.3249
## [425] 466.7073 467.4449 355.2401 566.0383 708.5069 468.6414 469.0968 363.6720
## [433] 465.6666 448.9309 468.8193 468.4491 401.0437 468.3692 469.6884 451.7559
## [441] 397.7861 493.9194 471.4767 564.7502 471.5872 659.6306 470.5935 545.8293
## [449] 709.5163## 'data.frame': 10000 obs. of 4 variables:
## $ x : num 178605 178633 178661 178689 178718 ...
## $ y : num 329714 329714 329714 329714 329714 ...
## $ var1.pred: num 470 470 470 470 470 ...
## $ var1.var : num 128766 128766 128766 128766 128766 ...## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 273.6 469.7 469.7 469.7 469.7 965.2## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 65466 128766 128766 128517 128766 128766## [1] 449## [1] 451## x y var1.pred var1.var
## 1 178605.0 329714 469.7161 128766.2
## 2 178633.1 329714 469.7161 128766.2
## 3 178661.3 329714 469.7161 128766.2
## 4 178689.4 329714 469.7161 128766.2
## 5 178717.5 329714 469.7161 128766.2
## 6 178745.7 329714 469.7161 128766.2
## 7 178773.8 329714 469.7161 128766.2
## 8 178801.9 329714 469.7161 128766.2
## 9 178830.1 329714 469.7161 128766.2
## 10 178858.2 329714 469.7161 128766.2## model psill range
## 1 Nug 36191.27 0.00000
## 2 Sph 91749.49 31.74033## np dist gamma dir.hor dir.ver id
## 1 57 79.29244 37362.96 0 0 var1
## 2 299 163.97367 72718.34 0 0 var1
## 3 419 267.36483 82655.53 0 0 var1
## 4 457 372.73542 111575.91 0 0 var1
## 5 547 478.47670 123080.69 0 0 var1
## 6 533 585.34058 147414.28 0 0 var1
## 7 574 693.14526 142891.51 0 0 var1
## 8 564 796.18365 153563.90 0 0 var1
## 9 589 903.14650 160217.69 0 0 var1
## 10 543 1011.29177 168867.36 0 0 var1
## 11 500 1117.86235 186528.20 0 0 var1
## 12 477 1221.32810 150320.76 0 0 var1
## 13 452 1329.16407 180399.66 0 0 var1
## 14 457 1437.25620 143139.84 0 0 var1
## 15 415 1543.20248 144112.31 0 0 var1ggplot() +
geom_raster(
data = hasil_ok_df,
aes(x = x, y = y, fill = var1.pred)
) +
geom_point(
data = data_ok,
aes(x = x, y = y),
color = "black",
size = 0.8
) +
scale_fill_viridis_c(
option = "plasma",
name = "Prediksi Zinc"
) +
coord_equal() +
labs(
title = "Interpolasi Ordinary Kriging",
subtitle = "Prediksi konsentrasi zinc",
x = "Koordinat X",
y = "Koordinat Y"
) +
theme_minimal()
Berdasarkan hasil interpolasi menggunakan metode Ordinary Kriging,
diperoleh estimasi persebaran konsentrasi zinc dengan mempertimbangkan
hubungan ketergantungan spasial antar titik pengamatan. Berbeda dengan
metode interpolasi sebelumnya, Ordinary Kriging tidak hanya
mempertimbangkan jarak antar lokasi, tetapi juga memperhitungkan
struktur korelasi spasial melalui model semivariogram.
Pada penelitian ini, pemodelan semivariogram dilakukan dengan membandingkan tiga model, yaitu Spherical, Exponential, dan Gaussian. Pemilihan model semivariogram terbaik didasarkan pada nilai Sum of Squared Errors (SSE) terkecil. Model terpilih selanjutnya digunakan dalam proses Ordinary Kriging untuk memperoleh estimasi konsentrasi zinc pada lokasi yang tidak memiliki pengamatan.
Hasil interpolasi juga menghasilkan nilai varians prediksi (var1.var) yang menunjukkan tingkat ketidakpastian estimasi pada setiap lokasi. Dengan demikian, Ordinary Kriging tidak hanya memberikan nilai prediksi konsentrasi zinc, tetapi juga informasi mengenai ketidakpastian hasil interpolasi.
4 5. LEAVE-ONE-OUT CROSS VALIDATION
# 1. Menyiapkan data
data_cv <- data_meuse %>%
select(x, y, zinc) %>%
filter(complete.cases(.)) %>%
distinct(x, y, .keep_all = TRUE)
# Jumlah observasi
n_obs <- nrow(data_cv)
cat("Jumlah observasi:", n_obs, "\n")## Jumlah observasi: 155# 2. Menyiapkan penyimpanan hasil prediksi
hasil_cv <- data.frame(
Observasi = data_cv$zinc,
Closest = rep(NA_real_, n_obs),
NN = rep(NA_real_, n_obs),
IDW = rep(NA_real_, n_obs),
Kriging = rep(NA_real_, n_obs)
)
# 3. Melakukan LOOCV
for (i in seq_len(n_obs)) {
# Data training dan testing
data_train <- data_cv[-i, , drop = FALSE]
data_test <- data_cv[i, , drop = FALSE]
# A. Closest Observation
hasil_cv$Closest[i] <- tryCatch({
as.numeric(
pred_closest(data_train, data_test)
)
}, error = function(e) NA_real_)
# B. Natural Neighbor
hasil_cv$NN[i] <- tryCatch({
as.numeric(
pred_nn(data_train, data_test)
)
}, error = function(e) NA_real_)
# C. IDW
hasil_cv$IDW[i] <- tryCatch({
as.numeric(
pred_idw(
train = data_train,
test = data_test,
power = 2
)
)
}, error = function(e) NA_real_)
# D. Ordinary Kriging
hasil_cv$Kriging[i] <- tryCatch({
as.numeric(
pred_kriging(data_train, data_test)
)
}, error = function(e) NA_real_)
# Menampilkan progres
if (i %% 10 == 0 || i == n_obs) {
cat(
"Observasi ke-", i,
"dari", n_obs, "\n"
)
}
}## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 10 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 20 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 30 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 40 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 50 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 60 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 70 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 80 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 90 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 100 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 110 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 120 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 130 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 140 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 150 dari 155
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Observasi ke- 155 dari 155# 4. Menghitung prediksi Ensemble
hasil_cv$Ensemble <- rowMeans(
hasil_cv[, c("IDW", "NN", "Kriging")],
na.rm = FALSE
)
# 5. Melihat hasil prediksi
head(hasil_cv)## Observasi Closest NN IDW Kriging Ensemble
## 1 1022 1141 NA 793.8598 922.9235 NA
## 2 1141 1022 1009.2551 727.3774 880.7553 872.4626
## 3 640 1022 668.3533 605.2803 659.5549 644.3961
## 4 257 640 494.7849 477.0037 471.4921 481.0936
## 5 269 281 282.3858 368.8423 265.9119 305.7133
## 6 281 269 NA 347.1991 262.4199 NA## 'data.frame': 155 obs. of 6 variables:
## $ Observasi: num 1022 1141 640 257 269 ...
## $ Closest : num 1141 1022 1022 640 281 ...
## $ NN : num NA 1009 668 495 282 ...
## $ IDW : num 794 727 605 477 369 ...
## $ Kriging : num 923 881 660 471 266 ...
## $ Ensemble : num NA 872 644 481 306 ...# 6. Menghitung jumlah prediksi valid
jumlah_valid <- data.frame(
Metode = c(
"Closest Observation",
"Natural Neighbor",
"IDW",
"Ordinary Kriging",
"Ensemble"
),
Jumlah_Valid = sapply(
hasil_cv[, c(
"Closest",
"NN",
"IDW",
"Kriging",
"Ensemble"
)],
function(x) sum(is.finite(x))
)
)
jumlah_valid$Persentase <- round(
jumlah_valid$Jumlah_Valid / n_obs * 100,
2
)
print(jumlah_valid)## Metode Jumlah_Valid Persentase
## Closest Closest Observation 155 100.00
## NN Natural Neighbor 151 97.42
## IDW IDW 155 100.00
## Kriging Ordinary Kriging 155 100.00
## Ensemble Ensemble 151 97.425 6. PERBANDINGAN HASIL EVALUASI
# Menggunakan observasi valid yang sama
# untuk seluruh metode
data_lengkap <- hasil_cv[
complete.cases(hasil_cv),
]
# Memastikan semua nilai prediksi bersifat finite
data_lengkap <- data_lengkap[
apply(
data_lengkap,
1,
function(x) all(is.finite(x))
),
]
# Memeriksa jumlah observasi yang digunakan
cat(
"Jumlah data evaluasi:",
nrow(data_lengkap),
"\n"
)## Jumlah data evaluasi: 151# Daftar metode
metode <- c(
"Closest",
"NN",
"IDW",
"Kriging",
"Ensemble"
)
# Menyiapkan tabel evaluasi
hasil_evaluasi <- data.frame(
Metode = c(
"Closest Observation",
"Natural Neighbor",
"IDW",
"Ordinary Kriging",
"Ensemble"
),
RMSE = NA_real_,
MAE = NA_real_,
ME = NA_real_
)
# Menghitung ukuran kesalahan
for (i in seq_along(metode)) {
prediksi <- data_lengkap[[metode[i]]]
error <- data_lengkap$Observasi - prediksi
hasil_evaluasi$RMSE[i] <- sqrt(
mean(error^2)
)
hasil_evaluasi$MAE[i] <- mean(
abs(error)
)
hasil_evaluasi$ME[i] <- mean(error)
}
# Mengurutkan berdasarkan RMSE terkecil
hasil_evaluasi <- hasil_evaluasi[
order(hasil_evaluasi$RMSE),
]
# Membulatkan hasil
hasil_evaluasi[, -1] <- round(
hasil_evaluasi[, -1],
4
)
# Menampilkan tabel
print(hasil_evaluasi)## Metode RMSE MAE ME
## 2 Natural Neighbor 216.1357 139.7691 -16.1101
## 4 Ordinary Kriging 229.9060 151.9544 1.3869
## 5 Ensemble 230.9172 159.7996 -5.5144
## 3 IDW 280.6564 205.9729 -1.8201
## 1 Closest Observation 313.6804 194.2384 8.5430# Menentukan metode terbaik
cat(
"\nMetode terbaik berdasarkan RMSE:",
hasil_evaluasi$Metode[1],
"\n"
)##
## Metode terbaik berdasarkan RMSE: Natural Neighbor## RMSE terkecil: 216.13576 7. VISUALISASI PERBANDINGAN RMSE
library(ggplot2)
ggplot(
hasil_evaluasi,
aes(
x = reorder(Metode, RMSE),
y = RMSE
)
) +
geom_col(
fill = "steelblue",
width = 0.7
) +
geom_text(
aes(label = round(RMSE, 2)),
hjust = -0.15,
size = 4
) +
coord_flip() +
expand_limits(
y = max(hasil_evaluasi$RMSE) * 1.15
) +
labs(
title = "Perbandingan RMSE Metode Interpolasi",
subtitle = "Leave-One-Out Cross Validation",
x = "Metode Interpolasi",
y = "RMSE"
) +
theme_minimal()
Berdasarkan grafik perbandingan nilai Root Mean Square Error (RMSE)
menggunakan metode Leave-One-Out Cross Validation (LOOCV), terlihat
bahwa setiap metode interpolasi menghasilkan tingkat kesalahan prediksi
yang berbeda. Metode Closest Observation memiliki nilai RMSE tertinggi
sebesar 313,68, diikuti oleh IDW sebesar 280,66. Sementara itu, metode
Ensemble menghasilkan RMSE sebesar 230,92, dan Ordinary Kriging sebesar
229,91.
Metode Natural Neighbor memperoleh nilai RMSE terendah, yaitu 216,14, yang menunjukkan bahwa metode tersebut memiliki tingkat kesalahan prediksi paling kecil dibandingkan metode lainnya. Dengan demikian, berdasarkan hasil evaluasi LOOCV, Natural Neighbor merupakan metode interpolasi dengan kinerja prediksi terbaik, sedangkan Closest Observation memiliki tingkat kesalahan prediksi tertinggi pada dataset yang digunakan
7 Kesimpulan
Berdasarkan hasil penerapan lima metode interpolasi pada dataset Meuse, diperoleh bahwa setiap metode memiliki karakteristik dan tingkat akurasi yang berbeda. Evaluasi menggunakan Leave-One-Out Cross Validation (LOOCV) menunjukkan bahwa metode Natural Neighbor memiliki nilai RMSE terendah sebesar 216,14, sehingga menjadi metode dengan kinerja prediksi terbaik dibandingkan metode lainnya. Sementara itu, metode Closest Observation menghasilkan nilai RMSE tertinggi sebesar 313,68. Hasil ini menunjukkan bahwa metode Natural Neighbor memberikan estimasi konsentrasi zinc yang lebih akurat berdasarkan kriteria RMSE pada data yang dianalisis.