Tugas Regresi Spasial Praktikum 10
Fitria Angelita Delvi
2026-11-10
Persiapan Data
# Membaca dataset
data <- read.csv("C:\\Users\\Fitria Angelita D\\Downloads\\curah_hujan_jabar_2024-12-04_458_stasiun.csv")
# Melihat 6 baris pertama
head(data)## Date Station_Id Rainfall_station District Sub_District Latitude
## 1 2024-12-04 32010102a Ciriung Kab. Bogor Cibinong -6.460417
## 2 2024-12-04 32010603a BP3K Jonggol Kab. Bogor Jonggol -6.467889
## 3 2024-12-04 32010606a Dayeuh/Sukanagara Kab. Bogor Jonggol -6.509972
## 4 2024-12-04 32010701a Klapanunggal Kab. Bogor Cileungsi -6.414222
## 5 2024-12-04 32010704a Situ Tungilis Kab. Bogor Cileungsi -6.419806
## 6 2024-12-04 32010802a Cariu Kab. Bogor Cariu -6.504222
## Longitude Elev...m.. Missing_Data.... Rainfall_mm
## 1 106.8564 146 0.7755474 0
## 2 107.0659 137 0.0000000 2
## 3 107.0200 169 0.0000000 15
## 4 106.9621 90 0.0000000 11
## 5 107.0313 100 0.0000000 5
## 6 107.1329 128 0.0000000 2# Melihat struktur data
str(data)## 'data.frame': 458 obs. of 10 variables:
## $ Date : chr "2024-12-04" "2024-12-04" "2024-12-04" "2024-12-04" ...
## $ Station_Id : chr "32010102a" "32010603a" "32010606a" "32010701a" ...
## $ Rainfall_station: chr "Ciriung" "BP3K Jonggol" "Dayeuh/Sukanagara" "Klapanunggal" ...
## $ District : chr "Kab. Bogor" "Kab. Bogor" "Kab. Bogor" "Kab. Bogor" ...
## $ Sub_District : chr "Cibinong" "Jonggol" "Jonggol" "Cileungsi" ...
## $ Latitude : num -6.46 -6.47 -6.51 -6.41 -6.42 ...
## $ Longitude : num 107 107 107 107 107 ...
## $ Elev...m.. : num 146 137 169 90 100 128 115 122 88 144 ...
## $ Missing_Data....: num 0.776 0 0 0 0 ...
## $ Rainfall_mm : num 0 2 15 11 5 2 5 9.2 5 19 ...# Melihat dimensi data
dim(data)## [1] 458 10# Ringkasan statistik
summary(data)## Date Station_Id Rainfall_station District
## Length:458 Length:458 Length:458 Length:458
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## Sub_District Latitude Longitude Elev...m..
## Length:458 Min. :-7.765 Min. :106.4 Min. : 0.0
## Class :character 1st Qu.:-7.039 1st Qu.:107.1 1st Qu.: 36.0
## Mode :character Median :-6.776 Median :107.6 Median : 197.0
## Mean :-6.793 Mean :107.6 Mean : 347.7
## 3rd Qu.:-6.492 3rd Qu.:108.2 3rd Qu.: 579.0
## Max. :-6.022 Max. :108.8 Max. :1880.0
## NA's :9
## Missing_Data.... Rainfall_mm
## Min. : 0.0000 Min. : 0.00
## 1st Qu.: 0.0000 1st Qu.: 3.00
## Median : 0.5018 Median : 14.00
## Mean : 2.1936 Mean : 39.84
## 3rd Qu.: 0.5931 3rd Qu.: 46.88
## Max. :59.6715 Max. :357.00
## Berdasarkan hasil pembacaan dataset, diperoleh sebanyak 458 observasi dan 10 variabel yang mencakup informasi stasiun pengamatan, lokasi geografis, serta curah hujan di Provinsi Jawa Barat pada 4 Desember 2024. Variabel utama yang digunakan dalam analisis interpolasi adalah Latitude, Longitude, dan Rainfall_mm. Hasil statistik deskriptif menunjukkan bahwa curah hujan memiliki nilai minimum sebesar 0 mm, maksimum 357 mm, rata-rata 39,84 mm, dan median 14 mm. Nilai rata-rata yang lebih besar dibandingkan median mengindikasikan bahwa distribusi curah hujan cenderung menceng ke kanan. Dengan demikian, dataset memiliki variasi nilai curah hujan yang cukup besar antarstasiun pengamatan dan dapat digunakan untuk analisis interpolasi spasial.
Periksa Kualitas Data
# 1. Memeriksa missing value setiap variabel
colSums(is.na(data))## Date Station_Id Rainfall_station District
## 0 0 0 0
## Sub_District Latitude Longitude Elev...m..
## 0 0 0 9
## Missing_Data.... Rainfall_mm
## 0 0# 2. Memeriksa missing value variabel utama
sum(is.na(data$Latitude))## [1] 0sum(is.na(data$Longitude))## [1] 0sum(is.na(data$Rainfall_mm))## [1] 0# 3. Memeriksa koordinat duplikat
sum(duplicated(data[, c("Latitude", "Longitude")]))## [1] 0# 4. Memeriksa jumlah stasiun unik
length(unique(data$Station_Id))## [1] 458# 5. Memeriksa ringkasan variabel utama
summary(data[, c(
"Latitude",
"Longitude",
"Rainfall_mm"
)])## Latitude Longitude Rainfall_mm
## Min. :-7.765 Min. :106.4 Min. : 0.00
## 1st Qu.:-7.039 1st Qu.:107.1 1st Qu.: 3.00
## Median :-6.776 Median :107.6 Median : 14.00
## Mean :-6.793 Mean :107.6 Mean : 39.84
## 3rd Qu.:-6.492 3rd Qu.:108.2 3rd Qu.: 46.88
## Max. :-6.022 Max. :108.8 Max. :357.00Hasil pemeriksaan kualitas data menunjukkan bahwa seluruh 458 observasi memiliki nilai Latitude, Longitude, dan Rainfall_mm yang lengkap tanpa missing value. Selain itu, tidak ditemukan koordinat stasiun yang duplikat dan seluruh observasi memiliki identitas stasiun yang unik. Terdapat sembilan nilai kosong pada variabel elevasi, tetapi hal tersebut tidak memengaruhi analisis karena elevasi tidak digunakan sebagai variabel utama interpolasi. Oleh karena itu, dataset dinilai telah memenuhi kebutuhan awal untuk analisis interpolasi spasial tanpa memerlukan penghapusan observasi.
Eksplorasi Data
Histogram Curah Hujan
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.5.2# Histogram curah hujan
ggplot(data, aes(x = Rainfall_mm)) +
geom_histogram(
bins = 30,
fill = "steelblue",
color = "white"
) +
labs(
title = "Distribusi Curah Hujan Jawa Barat",
x = "Curah Hujan (mm)",
y = "Frekuensi"
) +
theme_minimal()
Berdasarkan histogram, distribusi curah hujan menunjukkan pola menceng ke kanan (positively skewed). Sebagian besar stasiun pengamatan memiliki curah hujan relatif rendah, sedangkan frekuensi pengamatan dengan curah hujan tinggi semakin sedikit. Hal ini menunjukkan bahwa distribusi curah hujan tidak simetris dan terdapat beberapa lokasi dengan nilai curah hujan yang jauh lebih tinggi dibandingkan mayoritas pengamatan.
Boxplot Curah Hujan
# Boxplot curah hujan
ggplot(data, aes(y = Rainfall_mm)) +
geom_boxplot(
fill = "skyblue",
color = "black"
) +
labs(
title = "Boxplot Curah Hujan Jawa Barat",
y = "Curah Hujan (mm)"
) +
theme_minimal()
Berdasarkan boxplot, terlihat adanya sejumlah pencilan (outlier) pada bagian atas distribusi, yang menunjukkan beberapa stasiun memiliki nilai curah hujan jauh lebih tinggi dibandingkan sebagian besar stasiun lainnya. Kondisi ini mengindikasikan adanya variasi curah hujan yang cukup besar antar lokasi pengamatan. Pencilan tersebut tetap dipertahankan karena belum terdapat bukti bahwa nilai-nilai tersebut merupakan kesalahan pengukuran.
Visualisasi Persebaran Spasial
# Scatter plot spasial
ggplot(
data,
aes(
x = Longitude,
y = Latitude,
color = Rainfall_mm
)
) +
geom_point(size = 3) +
scale_color_viridis_c(
name = "Curah Hujan (mm)"
) +
labs(
title = "Persebaran Spasial Curah Hujan Jawa Barat",
x = "Longitude",
y = "Latitude"
) +
theme_minimal()
Berdasarkan scatter plot, titik-titik pengamatan tersebar di berbagai wilayah Jawa Barat dengan nilai curah hujan yang bervariasi. Beberapa lokasi di bagian barat daya menunjukkan konsentrasi curah hujan yang relatif tinggi, sedangkan sebagian besar lokasi lainnya menunjukkan nilai yang lebih rendah. Perbedaan nilai antar lokasi tersebut memperlihatkan adanya variasi spasial curah hujan yang menjadi dasar untuk melakukan interpolasi guna memperkirakan nilai curah hujan pada lokasi yang tidak memiliki stasiun pengamatan.
Statistik Deskriptif
# Statistik deskriptif
mean(data$Rainfall_mm)## [1] 39.83581median(data$Rainfall_mm)## [1] 14sd(data$Rainfall_mm)## [1] 64.10964var(data$Rainfall_mm)## [1] 4110.047min(data$Rainfall_mm)## [1] 0max(data$Rainfall_mm)## [1] 357Berdasarkan hasil perhitungan statistik deskriptif, diperoleh rata-rata curah hujan sebesar 39,84 mm, median 14 mm, standar deviasi 64,11 mm, dan varians 4.110,047. Nilai minimum curah hujan adalah 0 mm, sedangkan nilai maksimum mencapai 357 mm. Rata-rata yang lebih besar dibandingkan median menunjukkan bahwa distribusi data cenderung menceng ke kanan. Sementara itu, standar deviasi yang relatif besar menunjukkan adanya keragaman nilai curah hujan yang cukup tinggi antarstasiun pengamatan.
Persiapan Data Spasial
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 TRUElibrary(gstat)## Warning: package 'gstat' was built under R version 4.5.3library(terra)## Warning: package 'terra' was built under R version 4.5.3## terra 1.9.11library(ggplot2)# Mengubah data menjadi objek spasial
data_sf <- st_as_sf(
data,
coords = c("Longitude", "Latitude"),
crs = 4326,
remove = FALSE
)
# Melihat objek spasial
data_sf## Simple feature collection with 458 features and 10 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 106.4131 ymin: -7.765417 xmax: 108.7504 ymax: -6.021667
## Geodetic CRS: WGS 84
## First 10 features:
## Date Station_Id Rainfall_station District Sub_District
## 1 2024-12-04 32010102a Ciriung Kab. Bogor Cibinong
## 2 2024-12-04 32010603a BP3K Jonggol Kab. Bogor Jonggol
## 3 2024-12-04 32010606a Dayeuh/Sukanagara Kab. Bogor Jonggol
## 4 2024-12-04 32010701a Klapanunggal Kab. Bogor Cileungsi
## 5 2024-12-04 32010704a Situ Tungilis Kab. Bogor Cileungsi
## 6 2024-12-04 32010802a Cariu Kab. Bogor Cariu
## 7 2024-12-04 32010803a Sinar Rasa_Tanjung rasa Kab. Bogor Tanjung sari
## 8 2024-12-04 32011002a Jampang Kab. Bogor Ciseeng
## 9 2024-12-04 32011102a Curug Serpong Kab. Bogor Gunung Sindur
## 10 2024-12-04 32011201a Situ Kemang Kab. Bogor Kemang
## Latitude Longitude Elev...m.. Missing_Data.... Rainfall_mm
## 1 -6.460417 106.8564 146 0.7755474 0.0
## 2 -6.467889 107.0659 137 0.0000000 2.0
## 3 -6.509972 107.0200 169 0.0000000 15.0
## 4 -6.414222 106.9621 90 0.0000000 11.0
## 5 -6.419806 107.0313 100 0.0000000 5.0
## 6 -6.504222 107.1329 128 0.0000000 2.0
## 7 -6.619778 107.1580 115 0.5018248 5.0
## 8 -6.481472 106.6947 122 0.0000000 9.2
## 9 -6.386444 106.7266 88 0.0000000 5.0
## 10 -6.495389 106.7416 144 1.3686131 19.0
## geometry
## 1 POINT (106.8564 -6.460417)
## 2 POINT (107.0659 -6.467889)
## 3 POINT (107.02 -6.509972)
## 4 POINT (106.9621 -6.414222)
## 5 POINT (107.0313 -6.419806)
## 6 POINT (107.1329 -6.504222)
## 7 POINT (107.158 -6.619778)
## 8 POINT (106.6947 -6.481472)
## 9 POINT (106.7266 -6.386444)
## 10 POINT (106.7416 -6.495389)# Mengubah CRS menjadi UTM zona 48S
data_utm <- st_transform(
data_sf,
crs = 32748
)
# Melihat hasil transformasi
data_utm## Simple feature collection with 458 features and 10 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 655978.5 ymin: 9140143 xmax: 913928.5 ymax: 9333927
## Projected CRS: WGS 84 / UTM zone 48S
## First 10 features:
## Date Station_Id Rainfall_station District Sub_District
## 1 2024-12-04 32010102a Ciriung Kab. Bogor Cibinong
## 2 2024-12-04 32010603a BP3K Jonggol Kab. Bogor Jonggol
## 3 2024-12-04 32010606a Dayeuh/Sukanagara Kab. Bogor Jonggol
## 4 2024-12-04 32010701a Klapanunggal Kab. Bogor Cileungsi
## 5 2024-12-04 32010704a Situ Tungilis Kab. Bogor Cileungsi
## 6 2024-12-04 32010802a Cariu Kab. Bogor Cariu
## 7 2024-12-04 32010803a Sinar Rasa_Tanjung rasa Kab. Bogor Tanjung sari
## 8 2024-12-04 32011002a Jampang Kab. Bogor Ciseeng
## 9 2024-12-04 32011102a Curug Serpong Kab. Bogor Gunung Sindur
## 10 2024-12-04 32011201a Situ Kemang Kab. Bogor Kemang
## Latitude Longitude Elev...m.. Missing_Data.... Rainfall_mm
## 1 -6.460417 106.8564 146 0.7755474 0.0
## 2 -6.467889 107.0659 137 0.0000000 2.0
## 3 -6.509972 107.0200 169 0.0000000 15.0
## 4 -6.414222 106.9621 90 0.0000000 11.0
## 5 -6.419806 107.0313 100 0.0000000 5.0
## 6 -6.504222 107.1329 128 0.0000000 2.0
## 7 -6.619778 107.1580 115 0.5018248 5.0
## 8 -6.481472 106.6947 122 0.0000000 9.2
## 9 -6.386444 106.7266 88 0.0000000 5.0
## 10 -6.495389 106.7416 144 1.3686131 19.0
## geometry
## 1 POINT (705304.9 9285525)
## 2 POINT (728476.8 9284609)
## 3 POINT (723378 9279975)
## 4 POINT (717014.7 9290591)
## 5 POINT (724668.6 9289943)
## 6 POINT (735879.2 9280560)
## 7 POINT (738593.2 9267765)
## 8 POINT (687403.5 9283259)
## 9 POINT (690975.7 9293757)
## 10 POINT (692588 9281702)# Mengambil koordinat X dan Y
koordinat <- st_coordinates(data_utm)
# Menambahkan koordinat ke dataset
data_utm$X <- koordinat[, 1]
data_utm$Y <- koordinat[, 2]
# Melihat hasil
head(data_utm)## Simple feature collection with 6 features and 12 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 705304.9 ymin: 9279975 xmax: 735879.2 ymax: 9290591
## Projected CRS: WGS 84 / UTM zone 48S
## Date Station_Id Rainfall_station District Sub_District Latitude
## 1 2024-12-04 32010102a Ciriung Kab. Bogor Cibinong -6.460417
## 2 2024-12-04 32010603a BP3K Jonggol Kab. Bogor Jonggol -6.467889
## 3 2024-12-04 32010606a Dayeuh/Sukanagara Kab. Bogor Jonggol -6.509972
## 4 2024-12-04 32010701a Klapanunggal Kab. Bogor Cileungsi -6.414222
## 5 2024-12-04 32010704a Situ Tungilis Kab. Bogor Cileungsi -6.419806
## 6 2024-12-04 32010802a Cariu Kab. Bogor Cariu -6.504222
## Longitude Elev...m.. Missing_Data.... Rainfall_mm geometry
## 1 106.8564 146 0.7755474 0 POINT (705304.9 9285525)
## 2 107.0659 137 0.0000000 2 POINT (728476.8 9284609)
## 3 107.0200 169 0.0000000 15 POINT (723378 9279975)
## 4 106.9621 90 0.0000000 11 POINT (717014.7 9290591)
## 5 107.0313 100 0.0000000 5 POINT (724668.6 9289943)
## 6 107.1329 128 0.0000000 2 POINT (735879.2 9280560)
## X Y
## 1 705304.9 9285525
## 2 728476.8 9284609
## 3 723378.0 9279975
## 4 717014.7 9290591
## 5 724668.6 9289943
## 6 735879.2 9280560# Memeriksa ringkasan koordinat
summary(data_utm$X)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 655978 729615 782802 787676 853372 913929summary(data_utm$Y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 9140143 9221179 9249824 9248276 9281652 9333927Berdasarkan hasil transformasi data, sebanyak 458 observasi berhasil dikonversi menjadi objek spasial (simple feature) dengan tipe geometri POINT, yang menunjukkan bahwa setiap observasi merepresentasikan satu lokasi stasiun pengamatan. Pada awalnya, data menggunakan sistem koordinat geografis WGS 84 (EPSG:4326) dengan rentang longitude 106,4131–108,7504 dan latitude -7,7654 hingga -6,0217. Selanjutnya, dilakukan transformasi sistem koordinat menjadi WGS 84 / UTM Zone 48S (EPSG:32748) sehingga koordinat lokasi dinyatakan dalam satuan meter. Hasil transformasi menghasilkan koordinat X dengan rentang 655.978,5–913.928,5 meter dan koordinat Y sebesar 9.140.143–9.333.927 meter. Penambahan variabel X dan Y memungkinkan perhitungan jarak antarstasiun dilakukan secara lebih tepat dalam satuan meter, sehingga data telah siap digunakan untuk proses interpolasi spasial.
Grid Prediksi
# Membuat convex hull dari seluruh titik stasiun
wilayah <- st_convex_hull(
st_union(data_utm)
)
# Mengubah menjadi sf
wilayah <- st_sf(
geometry = st_sfc(wilayah, crs = 32748)
)
# Melihat wilayah cakupan
plot(st_geometry(wilayah))
# Membuat grid dengan 100 baris dan 100 kolom
grid <- st_make_grid(
wilayah,
n = c(100, 100),
what = "centers"
)
# Mengubah grid menjadi objek sf
grid_sf <- st_sf(
grid_id = 1:length(grid),
geometry = grid
)
# Mempertahankan titik grid yang berada
# di dalam wilayah cakupan
grid_sf <- grid_sf[
lengths(st_within(grid_sf, wilayah)) > 0,
]
# Melihat jumlah titik prediksi
nrow(grid_sf)## [1] 7552# Mengambil koordinat grid
koordinat_grid <- st_coordinates(grid_sf)
# Menambahkan koordinat X dan Y
grid_sf$X <- koordinat_grid[, 1]
grid_sf$Y <- koordinat_grid[, 2]
# Melihat hasil grid
head(grid_sf)## Simple feature collection with 6 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 855889.8 ymin: 9141112 xmax: 876525.8 ymax: 9143050
## Projected CRS: WGS 84 / UTM zone 48S
## grid_id geometry X Y
## 83 83 POINT (868787.3 9141112) 868787.3 9141112
## 84 84 POINT (871366.8 9141112) 871366.8 9141112
## 85 85 POINT (873946.3 9141112) 873946.3 9141112
## 86 86 POINT (876525.8 9141112) 876525.8 9141112
## 178 178 POINT (855889.8 9143050) 855889.8 9143050
## 179 179 POINT (858469.3 9143050) 858469.3 9143050ggplot() +
geom_sf(
data = wilayah,
fill = NA,
color = "black"
) +
geom_sf(
data = grid_sf,
color = "red",
size = 0.3,
alpha = 0.5
) +
geom_sf(
data = data_utm,
color = "blue",
size = 1
) +
labs(
title = "Lokasi Pengamatan dan Grid Prediksi",
subtitle = "Biru: Stasiun pengamatan | Merah: Titik prediksi",
x = "Easting (m)",
y = "Northing (m)"
) +
theme_minimal()
Berdasarkan hasil pembentukan grid, diperoleh sebanyak 7.552 titik prediksi yang berada di dalam wilayah cakupan pengamatan. Grid dibentuk menggunakan pendekatan convex hull dari 458 stasiun pengamatan dengan pembagian awal sebanyak 100 baris dan 100 kolom. Titik-titik prediksi ditunjukkan dengan warna merah, sedangkan titik pengamatan ditunjukkan dengan warna biru. Visualisasi memperlihatkan bahwa grid prediksi tersebar di seluruh wilayah cakupan dan mencakup area di antara titik-titik stasiun pengamatan. Titik-titik tersebut akan digunakan sebagai lokasi untuk memperkirakan nilai curah hujan menggunakan berbagai metode interpolasi spasial.
Interpolasi Closest Observation
# Mengambil koordinat titik pengamatan
xy_obs <- st_coordinates(data_utm)
# Mengambil koordinat titik prediksi
xy_grid <- st_coordinates(grid_sf)
# Menghitung jarak setiap grid
# terhadap seluruh titik pengamatan
jarak <- as.matrix(
dist(rbind(xy_obs, xy_grid))
)
# Mengambil bagian jarak grid ke observasi
jarak_grid <- jarak[
(nrow(xy_obs) + 1):nrow(jarak),
1:nrow(xy_obs)
]
# Mencari indeks titik terdekat
indeks_terdekat <- apply(
jarak_grid,
1,
which.min
)
# Menentukan nilai prediksi
grid_sf$Closest <- data_utm$Rainfall_mm[
indeks_terdekat
]
# Melihat hasil prediksi
head(grid_sf)## Simple feature collection with 6 features and 4 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 855889.8 ymin: 9141112 xmax: 876525.8 ymax: 9143050
## Projected CRS: WGS 84 / UTM zone 48S
## grid_id geometry X Y Closest
## 83 83 POINT (868787.3 9141112) 868787.3 9141112 0
## 84 84 POINT (871366.8 9141112) 871366.8 9141112 0
## 85 85 POINT (873946.3 9141112) 873946.3 9141112 0
## 86 86 POINT (876525.8 9141112) 876525.8 9141112 8
## 178 178 POINT (855889.8 9143050) 855889.8 9143050 12
## 179 179 POINT (858469.3 9143050) 858469.3 9143050 1Berdasarkan enam observasi pertama pada hasil prediksi, diperoleh nilai curah hujan sebesar 0, 0, 0, 8, 12, dan 1 mm. Nilai tersebut merupakan curah hujan dari stasiun pengamatan yang memiliki jarak paling dekat dengan masing-masing titik prediksi. Hal ini menunjukkan bahwa metode Closest Observation tidak menghitung rata-rata nilai dari beberapa stasiun, melainkan secara langsung menggunakan nilai dari satu stasiun terdekat.
ggplot() +
geom_raster(
data = grid_sf,
aes(
x = X,
y = Y,
fill = Closest
)
) +
geom_sf(
data = data_utm,
color = "black",
size = 0.3
) +
scale_fill_viridis_c(
name = "Curah Hujan (mm)"
) +
labs(
title = "Interpolasi Closest Observation",
x = "Easting (m)",
y = "Northing (m)"
) +
theme_minimal()
Berdasarkan hasil interpolasi Closest Observation, terlihat bahwa persebaran curah hujan di wilayah Jawa Barat menunjukkan variasi spasial yang cukup jelas. Wilayah bagian barat daya cenderung memiliki nilai curah hujan lebih tinggi, ditandai dengan warna hijau hingga kuning, sedangkan sebagian besar wilayah lainnya menunjukkan curah hujan yang relatif rendah dengan warna ungu. Pola prediksi terlihat membentuk area-area dengan batas yang cukup tegas karena setiap lokasi prediksi menggunakan nilai dari satu stasiun pengamatan terdekat. Dengan demikian, metode ini menghasilkan pola interpolasi yang mempertahankan nilai pengamatan asli tanpa melakukan perataan atau pembobotan antarstasiun.
summary(grid_sf$Closest)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 3.00 16.00 48.28 61.00 357.00sd(grid_sf$Closest)## [1] 71.99071Berdasarkan hasil statistik deskriptif, nilai prediksi curah hujan memiliki nilai minimum sebesar 0 mm, kuartil pertama sebesar 3 mm, median 16 mm, rata-rata 48,28 mm, kuartil ketiga 61 mm, dan maksimum 357 mm. Standar deviasi sebesar 71,99 mm menunjukkan bahwa hasil prediksi memiliki variasi nilai yang cukup besar. Rata-rata yang lebih tinggi dibandingkan median juga menunjukkan bahwa distribusi hasil prediksi cenderung menceng ke kanan, yang dipengaruhi oleh keberadaan beberapa lokasi dengan curah hujan tinggi.
Nearest Neighbors (NN)
library(RANN)## Warning: package 'RANN' was built under R version 4.5.3# Mengambil koordinat pengamatan
xy_obs <- st_coordinates(data_utm)
# Mengambil koordinat lokasi prediksi
xy_grid <- st_coordinates(grid_sf)
# Mencari 5 tetangga terdekat
nn_result <- nn2(
data = xy_obs,
query = xy_grid,
k = 5
)
# Menghitung rata-rata curah hujan
# dari 5 tetangga terdekat
grid_sf$NN <- apply(
nn_result$nn.idx,
1,
function(i) {
mean(data_utm$Rainfall_mm[i])
}
)
# Melihat hasil prediksi
head(grid_sf)## Simple feature collection with 6 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 855889.8 ymin: 9141112 xmax: 876525.8 ymax: 9143050
## Projected CRS: WGS 84 / UTM zone 48S
## grid_id geometry X Y Closest NN
## 83 83 POINT (868787.3 9141112) 868787.3 9141112 0 4.0
## 84 84 POINT (871366.8 9141112) 871366.8 9141112 0 3.8
## 85 85 POINT (873946.3 9141112) 873946.3 9141112 0 3.8
## 86 86 POINT (876525.8 9141112) 876525.8 9141112 8 3.8
## 178 178 POINT (855889.8 9143050) 855889.8 9143050 12 4.8
## 179 179 POINT (858469.3 9143050) 858469.3 9143050 1 4.8Berdasarkan enam observasi pertama, diperoleh nilai prediksi curah hujan sebesar 4,0; 3,8; 3,8; 3,8; 4,8; dan 4,8 mm. Nilai tersebut merupakan rata-rata curah hujan dari lima stasiun pengamatan terdekat pada masing-masing lokasi prediksi. Hasil ini menunjukkan bahwa metode NN menggabungkan informasi dari beberapa stasiun sehingga nilai prediksi tidak harus sama dengan nilai salah satu stasiun pengamatan.
ggplot() +
geom_raster(
data = grid_sf,
aes(
x = X,
y = Y,
fill = NN
)
) +
geom_sf(
data = data_utm,
color = "black",
size = 0.3
) +
scale_fill_viridis_c(
name = "Curah Hujan (mm)"
) +
labs(
title = "Interpolasi Nearest Neighbors (k = 5)",
x = "Easting (m)",
y = "Northing (m)"
) +
theme_minimal()
Berdasarkan peta hasil interpolasi Nearest Neighbors dengan \(k=5\), terlihat bahwa pola persebaran curah hujan masih menunjukkan konsentrasi nilai tinggi di wilayah barat daya, sedangkan sebagian besar wilayah lainnya memiliki curah hujan relatif rendah. Dibandingkan metode Closest Observation, hasil interpolasi NN terlihat lebih halus dengan perubahan nilai yang lebih bertahap. Hal ini terjadi karena nilai prediksi diperoleh dari rata-rata lima stasiun pengamatan terdekat, sehingga pengaruh nilai ekstrem dari satu stasiun dapat dikurangi melalui proses perataan.
summary(grid_sf$NN)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 4.63 19.90 48.17 59.20 325.40sd(grid_sf$NN)## [1] 67.3742Hasil statistik deskriptif menunjukkan bahwa nilai prediksi memiliki minimum 0 mm, kuartil pertama 4,63 mm, median 19,90 mm, rata-rata 48,17 mm, kuartil ketiga 59,20 mm, dan maksimum 325,40 mm. Standar deviasi sebesar 67,37 mm menunjukkan bahwa hasil prediksi masih memiliki variasi yang cukup besar. Dibandingkan Closest Observation, nilai maksimum dan standar deviasi NN lebih rendah, yang mengindikasikan bahwa penggunaan rata-rata lima tetangga terdekat menghasilkan prediksi yang relatif lebih halus dan mengurangi pengaruh nilai curah hujan ekstrem.
Inverse Distance Weighting (IDW)
# Mengambil koordinat pengamatan
xy_obs <- as.matrix(
st_drop_geometry(data_utm)[, c("X", "Y")]
)
# Mengambil koordinat lokasi prediksi
xy_grid <- as.matrix(
st_drop_geometry(grid_sf)[, c("X", "Y")]
)
# Memastikan tipe data numerik
storage.mode(xy_obs) <- "double"
storage.mode(xy_grid) <- "double"
# Memeriksa struktur
str(xy_obs)## num [1:458, 1:2] 705305 728477 723378 717015 724669 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:458] "1" "2" "3" "4" ...
## ..$ : chr [1:2] "X" "Y"str(xy_grid)## num [1:7552, 1:2] 868787 871367 873946 876526 855890 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:7552] "83" "84" "85" "86" ...
## ..$ : chr [1:2] "X" "Y"Berdasarkan hasil matriks koordinat, terdapat 458 titik pengamatan dan 7.552 titik prediksi yang berhasil diproses. Seluruh lokasi prediksi memperoleh nilai estimasi curah hujan berdasarkan 12 tetangga terdekat dengan pembobotan invers jarak berpangkat dua. Hal ini menunjukkan bahwa setiap lokasi prediksi mempertimbangkan kontribusi beberapa stasiun pengamatan, dengan stasiun yang lebih dekat memiliki pengaruh lebih besar terhadap nilai estimasi.
# Mencari 12 tetangga terdekat
idw_result <- nn2(
data = xy_obs,
query = xy_grid,
k = 12
)
# Menghitung prediksi IDW
grid_sf$IDW <- sapply(
seq_len(nrow(xy_grid)),
function(j) {
# Jarak ke 12 tetangga terdekat
d <- idw_result$nn.dists[j, ]
# Indeks tetangga terdekat
idx <- idw_result$nn.idx[j, ]
# Mengambil nilai curah hujan
nilai <- data_utm$Rainfall_mm[idx]
# Jika terdapat jarak nol
if (any(d == 0)) {
return(nilai[which.min(d)])
}
# Menghitung bobot IDW
bobot <- 1 / (d^2)
# Menghitung prediksi
sum(bobot * nilai) / sum(bobot)
}
)ggplot() +
geom_raster(
data = grid_sf,
aes(x = X, y = Y, fill = IDW)
) +
geom_sf(
data = data_utm,
color = "black",
size = 0.3
) +
scale_fill_viridis_c(
name = "Curah Hujan (mm)"
) +
labs(
title = "Interpolasi Inverse Distance Weighting (IDW)",
subtitle = "Power = 2, Neighbors = 12",
x = "Easting (m)",
y = "Northing (m)"
) +
theme_minimal()
Berdasarkan peta hasil interpolasi IDW dengan parameter power sebesar 2 dan 12 tetangga terdekat, terlihat bahwa wilayah barat daya Jawa Barat cenderung memiliki curah hujan lebih tinggi dibandingkan wilayah lainnya. Area dengan curah hujan tinggi ditunjukkan oleh warna hijau hingga kuning, sedangkan wilayah dengan curah hujan rendah didominasi warna ungu. Perubahan nilai curah hujan pada peta terlihat lebih halus dibandingkan metode Closest Observation karena IDW memberikan bobot lebih besar kepada stasiun yang jaraknya lebih dekat dengan lokasi prediksi. Dengan demikian, nilai prediksi dipengaruhi oleh jarak dan nilai curah hujan dari 12 stasiun pengamatan terdekat.
summary(grid_sf$IDW)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00354 4.49275 18.73374 47.65929 60.71293 350.97800sd(grid_sf$IDW)## [1] 66.15567Hasil statistik deskriptif menunjukkan bahwa nilai prediksi IDW memiliki minimum sebesar 0,00354 mm, kuartil pertama 4,49275 mm, median 18,73374 mm, rata-rata 47,65929 mm, kuartil ketiga 60,71293 mm, dan maksimum 350,978 mm. Standar deviasi sebesar 66,15567 mm menunjukkan bahwa hasil prediksi masih memiliki variasi curah hujan yang cukup besar. Dibandingkan metode Closest Observation dan NN, IDW menghasilkan nilai maksimum dan standar deviasi yang lebih rendah. Hal ini menunjukkan bahwa pembobotan berdasarkan jarak menghasilkan permukaan prediksi yang relatif lebih halus dengan pengaruh nilai ekstrem yang lebih terkendali.
Ensemble Interpolation
# Menggabungkan tiga hasil interpolasi
grid_sf$Ensemble <- (
grid_sf$Closest +
grid_sf$NN +
grid_sf$IDW
) / 3
# Melihat hasil prediksi
head(grid_sf)## Simple feature collection with 6 features and 7 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 855889.8 ymin: 9141112 xmax: 876525.8 ymax: 9143050
## Projected CRS: WGS 84 / UTM zone 48S
## grid_id geometry X Y Closest NN IDW
## 83 83 POINT (868787.3 9141112) 868787.3 9141112 0 4.0 2.2446509
## 84 84 POINT (871366.8 9141112) 871366.8 9141112 0 3.8 1.0557909
## 85 85 POINT (873946.3 9141112) 873946.3 9141112 0 3.8 0.4017983
## 86 86 POINT (876525.8 9141112) 876525.8 9141112 8 3.8 4.2317771
## 178 178 POINT (855889.8 9143050) 855889.8 9143050 12 4.8 9.6180596
## 179 179 POINT (858469.3 9143050) 858469.3 9143050 1 4.8 4.3985455
## Ensemble
## 83 2.081550
## 84 1.618597
## 85 1.400599
## 86 5.343926
## 178 8.806020
## 179 3.399515Berdasarkan enam observasi pertama, diperoleh nilai prediksi Ensemble sebesar 2,081550; 1,618597; 1,400599; 5,343926; 8,806020; dan 3,399515 mm. Nilai tersebut diperoleh dari rata-rata hasil prediksi Closest Observation, NN, dan IDW pada masing-masing lokasi. Sebagai contoh, pada grid pertama, nilai prediksi sebesar 2,081550 mm diperoleh dari rata-rata nilai Closest sebesar 0 mm, NN sebesar 4 mm, dan IDW sebesar 2,244651 mm. Hal ini menunjukkan bahwa Ensemble menggabungkan informasi dari ketiga metode untuk menghasilkan satu estimasi curah hujan.
ggplot() +
geom_raster(
data = grid_sf,
aes(
x = X,
y = Y,
fill = Ensemble
)
) +
geom_sf(
data = data_utm,
color = "black",
size = 0.3
) +
scale_fill_viridis_c(
name = "Curah Hujan (mm)"
) +
labs(
title = "Interpolasi Ensemble",
subtitle = "Closest Observation + NN + IDW",
x = "Easting (m)",
y = "Northing (m)"
) +
theme_minimal()
Berdasarkan peta hasil interpolasi Ensemble yang menggabungkan metode Closest Observation, Nearest Neighbors, dan IDW dengan bobot yang sama, terlihat bahwa wilayah barat daya Jawa Barat masih menunjukkan konsentrasi curah hujan relatif tinggi, sedangkan sebagian besar wilayah lainnya memiliki curah hujan lebih rendah. Permukaan prediksi yang dihasilkan memperlihatkan pola perubahan nilai yang relatif halus karena menggabungkan karakteristik ketiga metode interpolasi. Penggabungan tersebut memungkinkan hasil prediksi memanfaatkan nilai dari titik terdekat, rata-rata beberapa tetangga, serta pembobotan berdasarkan jarak sehingga menghasilkan estimasi yang mengombinasikan ketiga pendekatan tersebut.
summary(grid_sf$Ensemble)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00692 4.06802 18.17391 48.03424 60.70285 344.11811sd(grid_sf$Ensemble)## [1] 67.48812Hasil statistik deskriptif menunjukkan bahwa nilai prediksi Ensemble memiliki minimum sebesar 0,00692 mm, kuartil pertama 4,06802 mm, median 18,17391 mm, rata-rata 48,03424 mm, kuartil ketiga 60,70285 mm, dan maksimum 344,11811 mm. Standar deviasi sebesar 67,48812 mm menunjukkan bahwa hasil prediksi masih memiliki variasi curah hujan yang cukup besar. Nilai rata-rata yang lebih tinggi dibandingkan median menunjukkan bahwa distribusi hasil prediksi masih cenderung menceng ke kanan. Hal ini mengindikasikan bahwa penggabungan ketiga metode tetap mempertahankan variasi spasial curah hujan, termasuk pengaruh beberapa lokasi dengan curah hujan tinggi.
Ordinary Kriging
library(gstat)
# Membuat variogram empiris
variogram_emp <- variogram(
Rainfall_mm ~ 1,
locations = ~ X + Y,
data = st_drop_geometry(data_utm)
)
# Melihat hasil variogram
head(variogram_emp)## np dist gamma dir.hor dir.ver id
## 1 589 4719.023 534.1467 0 0 var1
## 2 1486 11062.531 639.3882 0 0 var1
## 3 2276 18159.840 945.1056 0 0 var1
## 4 2992 25169.898 1266.1623 0 0 var1
## 5 3540 32332.649 1472.6924 0 0 var1
## 6 3898 39477.388 1974.0331 0 0 var1plot(
variogram_emp,
main = "Variogram Empiris Curah Hujan"
)
Berdasarkan grafik variogram empiris, terlihat bahwa nilai semivarians cenderung meningkat seiring bertambahnya jarak antarstasiun pengamatan. Pada jarak yang relatif dekat, semivarians masih rendah, kemudian meningkat hingga mencapai lebih dari 4.000 pada jarak sekitar 100 km. Pola ini menunjukkan adanya indikasi ketergantungan spasial pada data curah hujan, yaitu perbedaan nilai curah hujan cenderung semakin besar ketika jarak antarstasiun semakin jauh. Dengan demikian, variogram empiris memberikan gambaran mengenai hubungan antara jarak geografis dan tingkat keragaman nilai curah hujan yang menjadi dasar pembentukan model Kriging.
# Membuat model variogram awal
variogram_awal <- vgm(
psill = var(data_utm$Rainfall_mm) * 0.8,
model = "Exp",
range = 20000,
nugget = var(data_utm$Rainfall_mm) * 0.2
)
# Fitting variogram
variogram_fit <- fit.variogram(
variogram_emp,
variogram_awal
)
# Melihat parameter hasil fitting
variogram_fit## model psill range
## 1 Nug 304.9447 0
## 2 Exp 50311.5499 1294147plot(
variogram_emp,
variogram_fit,
main = "Fitting Variogram Exponential"
)
Berdasarkan hasil fitting, diperoleh model variogram Exponential dengan nilai nugget sebesar 304,9447, partial sill sebesar 50.311,5499, dan range sebesar 1.294.147 meter. Grafik menunjukkan bahwa model teoritis mengikuti kecenderungan peningkatan semivarians pada variogram empiris. Namun, nilai range hasil fitting jauh lebih besar dibandingkan rentang jarak yang ditampilkan pada grafik variogram empiris. Hal ini menunjukkan bahwa model belum memperlihatkan struktur sill yang tercapai dalam rentang pengamatan, sehingga kecocokan model variogram perlu dievaluasi lebih lanjut.
# Membuat model Ordinary Kriging
kriging_model <- gstat(
formula = Rainfall_mm ~ 1,
locations = ~ X + Y,
data = st_drop_geometry(data_utm),
model = variogram_fit
)
# Melakukan prediksi
kriging_result <- predict(
kriging_model,
newdata = st_drop_geometry(grid_sf)
)## [using ordinary kriging]# Menyimpan hasil prediksi
grid_sf$Kriging <- kriging_result$var1.pred
# Melihat hasil
head(grid_sf)## Simple feature collection with 6 features and 8 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 855889.8 ymin: 9141112 xmax: 876525.8 ymax: 9143050
## Projected CRS: WGS 84 / UTM zone 48S
## grid_id geometry X Y Closest NN IDW
## 83 83 POINT (868787.3 9141112) 868787.3 9141112 0 4.0 2.2446509
## 84 84 POINT (871366.8 9141112) 871366.8 9141112 0 3.8 1.0557909
## 85 85 POINT (873946.3 9141112) 873946.3 9141112 0 3.8 0.4017983
## 86 86 POINT (876525.8 9141112) 876525.8 9141112 8 3.8 4.2317771
## 178 178 POINT (855889.8 9143050) 855889.8 9143050 12 4.8 9.6180596
## 179 179 POINT (858469.3 9143050) 858469.3 9143050 1 4.8 4.3985455
## Ensemble Kriging
## 83 2.081550 11.021776
## 84 1.618597 7.793921
## 85 1.400599 5.567366
## 86 5.343926 6.066742
## 178 8.806020 27.947313
## 179 3.399515 23.523954Berdasarkan enam observasi pertama, diperoleh nilai prediksi Kriging sebesar 11,021776; 7,793921; 5,567366; 6,066742; 27,947313; dan 23,523954 mm. Nilai tersebut merupakan estimasi curah hujan pada lokasi grid yang tidak memiliki pengamatan langsung. Hasil ini menunjukkan bahwa Ordinary Kriging menghasilkan nilai prediksi berdasarkan kombinasi berbobot dari beberapa stasiun pengamatan dengan mempertimbangkan struktur variogram yang telah dibentuk.
ggplot() +
geom_raster(
data = grid_sf,
aes(
x = X,
y = Y,
fill = Kriging
)
) +
geom_sf(
data = data_utm,
color = "black",
size = 0.3
) +
scale_fill_viridis_c(
name = "Curah Hujan (mm)"
) +
labs(
title = "Interpolasi Ordinary Kriging",
subtitle = "Exponential Variogram",
x = "Easting (m)",
y = "Northing (m)"
) +
theme_minimal()
Berdasarkan peta hasil Ordinary Kriging menggunakan model variogram Exponential, terlihat bahwa wilayah barat daya Jawa Barat cenderung memiliki curah hujan lebih tinggi dibandingkan wilayah lainnya. Sementara itu, sebagian besar wilayah tengah, utara, dan timur menunjukkan nilai curah hujan yang relatif rendah. Permukaan prediksi terlihat lebih halus dibandingkan metode Closest Observation karena Kriging memperhitungkan struktur ketergantungan spasial melalui variogram. Dengan demikian, nilai prediksi tidak hanya dipengaruhi oleh jarak, tetapi juga oleh struktur variasi spasial yang dimodelkan melalui variogram.
summary(grid_sf$Kriging)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.261 3.815 20.046 48.371 61.891 332.409sd(grid_sf$Kriging)## [1] 66.13495Hasil statistik deskriptif menunjukkan bahwa nilai prediksi Ordinary Kriging memiliki minimum sebesar -2,261 mm, kuartil pertama 3,815 mm, median 20,046 mm, rata-rata 48,371 mm, kuartil ketiga 61,891 mm, dan maksimum 332,409 mm. Standar deviasi sebesar 66,13495 mm menunjukkan adanya variasi nilai prediksi yang cukup besar. Nilai minimum yang negatif perlu diperhatikan karena curah hujan secara fisik tidak mungkin bernilai negatif. Hal ini dapat terjadi akibat karakteristik model Kriging yang tidak membatasi nilai prediksi pada rentang nonnegatif. Nilai tersebut perlu dievaluasi dalam proses validasi model.
5 Fold Cross Validation
set.seed(123)
# Jumlah observasi
n <- nrow(data_utm)
# Membuat 5 fold
fold <- sample(
rep(1:5, length.out = n)
)
# Menambahkan fold ke data
data_utm$fold <- fold
# Melihat jumlah data setiap fold
table(data_utm$fold)##
## 1 2 3 4 5
## 92 92 92 91 91Berdasarkan hasil pembagian data menggunakan metode 5-Fold Cross Validation, sebanyak 458 observasi berhasil dibagi menjadi lima kelompok dengan ukuran yang relatif seimbang. Fold pertama, kedua, dan ketiga masing-masing terdiri atas 92 observasi, sedangkan fold keempat dan kelima masing-masing terdiri atas 91 observasi. Pada setiap iterasi, satu fold digunakan sebagai data pengujian (testing data), sementara empat fold lainnya digunakan sebagai data pelatihan (training data). Dengan demikian, seluruh observasi akan memperoleh kesempatan menjadi data pengujian sebanyak satu kali, sehingga evaluasi performa metode interpolasi dapat dilakukan secara lebih menyeluruh.
Cross Validation Kelima Metode
library(RANN)
library(gstat)
prediksi_deterministik <- function(train, test) {
xy_train <- as.matrix(
st_drop_geometry(train)[, c("X", "Y")]
)
xy_test <- as.matrix(
st_drop_geometry(test)[, c("X", "Y")]
)
z_train <- train$Rainfall_mm
# Closest Observation
co <- nn2(
data = xy_train,
query = xy_test,
k = 1
)
pred_co <- z_train[co$nn.idx[, 1]]
# Nearest Neighbors (k = 5)
nn <- nn2(
data = xy_train,
query = xy_test,
k = 5
)
pred_nn <- rowMeans(
matrix(
z_train[nn$nn.idx],
nrow = nrow(xy_test),
ncol = 5
)
)
# IDW (12 neighbors, power = 2)
idw <- nn2(
data = xy_train,
query = xy_test,
k = 12
)
pred_idw <- sapply(
seq_len(nrow(xy_test)),
function(j) {
d <- idw$nn.dists[j, ]
idx <- idw$nn.idx[j, ]
if (any(d == 0)) {
return(z_train[idx[which.min(d)]])
}
w <- 1 / d^2
sum(w * z_train[idx]) / sum(w)
}
)
# Ensemble
pred_ensemble <- (
pred_co + pred_nn + pred_idw
) / 3
data.frame(
CO = pred_co,
NN = pred_nn,
IDW = pred_idw,
Ensemble = pred_ensemble
)
}# Menyimpan hasil cross validation
hasil_cv <- data.frame()
for (k in 1:5) {
# Membagi data training dan testing
train <- data_utm[data_utm$fold != k, ]
test <- data_utm[data_utm$fold == k, ]
# Mengambil koordinat
xy_train <- as.matrix(
st_drop_geometry(train)[, c("X", "Y")]
)
xy_test <- as.matrix(
st_drop_geometry(test)[, c("X", "Y")]
)
z_train <- train$Rainfall_mm
# Closest Observation
co <- nn2(
data = xy_train,
query = xy_test,
k = 1
)
pred_co <- z_train[co$nn.idx[, 1]]
# Nearest Neighbors (k = 5)
nn <- nn2(
data = xy_train,
query = xy_test,
k = 5
)
pred_nn <- rowMeans(
matrix(
z_train[nn$nn.idx],
nrow = nrow(xy_test),
ncol = 5
)
)
# IDW (k = 12, power = 2)
idw <- nn2(
data = xy_train,
query = xy_test,
k = 12
)
pred_idw <- sapply(
seq_len(nrow(xy_test)),
function(j) {
d <- idw$nn.dists[j, ]
idx <- idw$nn.idx[j, ]
if (any(d == 0)) {
return(z_train[idx[which.min(d)]])
}
w <- 1 / d^2
sum(w * z_train[idx]) / sum(w)
}
)
# Ensemble
pred_ensemble <- (
pred_co + pred_nn + pred_idw
) / 3
# Ordinary Kriging
v_emp <- variogram(
Rainfall_mm ~ 1,
locations = ~ X + Y,
data = st_drop_geometry(train)
)
v_awal <- vgm(
psill = var(z_train) * 0.8,
model = "Exp",
range = 20000,
nugget = var(z_train) * 0.2
)
v_fit <- fit.variogram(
v_emp,
v_awal
)
krig_model <- gstat(
formula = Rainfall_mm ~ 1,
locations = ~ X + Y,
data = st_drop_geometry(train),
model = v_fit
)
krig_pred <- predict(
krig_model,
newdata = st_drop_geometry(test)
)
# Menggabungkan hasil prediksi
hasil_fold <- data.frame(
Fold = k,
Aktual = test$Rainfall_mm,
Closest = pred_co,
NN = pred_nn,
IDW = pred_idw,
Ensemble = pred_ensemble,
Kriging = krig_pred$var1.pred
)
hasil_cv <- rbind(
hasil_cv,
hasil_fold
)
}## Warning in fit.variogram(v_emp, v_awal): No convergence after 200 iterations:
## try different initial values?## [using ordinary kriging]
## [using ordinary kriging]## Warning in fit.variogram(v_emp, v_awal): No convergence after 200 iterations:
## try different initial values?## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]# Melihat hasil
head(hasil_cv)## Fold Aktual Closest NN IDW Ensemble Kriging
## 1 1 5 3.5 4.1 4.467416 4.022472 4.634771
## 2 1 2 3.0 3.1 4.036858 3.378953 3.769678
## 3 1 9 2.0 9.6 8.233318 6.611106 8.164779
## 4 1 8 7.0 6.8 8.060511 7.286837 9.684001
## 5 1 0 5.0 8.1 6.364251 6.488084 7.329685
## 6 1 6 12.5 6.3 8.064984 8.954995 8.808553# Memeriksa jumlah observasi
nrow(hasil_cv)## [1] 458# Memastikan jumlah observasi
nrow(hasil_cv)## [1] 458# Memeriksa missing value
colSums(is.na(hasil_cv))## Fold Aktual Closest NN IDW Ensemble Kriging
## 0 0 0 0 0 0 0# Melihat struktur hasil
str(hasil_cv)## 'data.frame': 458 obs. of 7 variables:
## $ Fold : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Aktual : num 5 2 9 8 0 6 7 14 252 295 ...
## $ Closest : num 3.5 3 2 7 5 12.5 6 33 61 305 ...
## $ NN : num 4.1 3.1 9.6 6.8 8.1 ...
## $ IDW : num 4.47 4.04 8.23 8.06 6.36 ...
## $ Ensemble: num 4.02 3.38 6.61 7.29 6.49 ...
## $ Kriging : num 4.63 3.77 8.16 9.68 7.33 ...# Menghitung RMSE setiap metode
rmse_cv <- data.frame(
Metode = c("Closest Observation", "Nearest Neighbor",
"IDW", "Ensemble", "Ordinary Kriging"),
RMSE = c(
sqrt(mean((hasil_cv$Aktual - hasil_cv$Closest)^2)),
sqrt(mean((hasil_cv$Aktual - hasil_cv$NN)^2)),
sqrt(mean((hasil_cv$Aktual - hasil_cv$IDW)^2)),
sqrt(mean((hasil_cv$Aktual - hasil_cv$Ensemble)^2)),
sqrt(mean((hasil_cv$Aktual - hasil_cv$Kriging)^2))
)
)
# Mengurutkan dari RMSE terkecil
rmse_cv <- rmse_cv[order(rmse_cv$RMSE), ]
rmse_cv## Metode RMSE
## 2 Nearest Neighbor 27.07995
## 5 Ordinary Kriging 27.88770
## 4 Ensemble 29.13232
## 3 IDW 30.78124
## 1 Closest Observation 34.32431Berdasarkan hasil 5-fold cross-validation, metode Nearest Neighbor menghasilkan nilai RMSE terkecil sebesar 27,07995 mm, sehingga menjadi metode dengan performa prediksi terbaik di antara kelima metode yang diuji. Ordinary Kriging menempati peringkat kedua dengan RMSE sebesar 27,88770 mm, yang menunjukkan kemampuan prediksi cukup baik dan mendekati Nearest Neighbor. Selanjutnya, Ensemble menghasilkan RMSE sebesar 29,13232 mm, diikuti IDW sebesar 30,78124 mm. Sementara itu, Closest Observation memiliki RMSE terbesar, yaitu 34,32431 mm, sehingga menghasilkan kesalahan prediksi paling tinggi. Hasil ini menunjukkan bahwa penggunaan rata-rata lima tetangga terdekat memberikan prediksi curah hujan yang lebih akurat dibandingkan metode lainnya pada data yang digunakan.
library(ggplot2)
ggplot(rmse_cv, aes(x = reorder(Metode, RMSE),
y = RMSE)) +
geom_col(fill = "steelblue") +
coord_flip() +
labs(
title = "Perbandingan RMSE Metode Interpolasi",
x = "Metode Interpolasi",
y = "RMSE (mm)"
) +
theme_minimal()
Berdasarkan diagram batang perbandingan nilai RMSE, metode Nearest Neighbor menghasilkan nilai RMSE terendah sebesar 27,07995 mm, sehingga menjadi metode dengan tingkat kesalahan prediksi paling kecil. Sebaliknya, metode Closest Observation menghasilkan nilai RMSE tertinggi sebesar 34,32431 mm. Ordinary Kriging dan Ensemble memiliki nilai RMSE yang relatif berdekatan, masing-masing sebesar 27,88770 mm dan 29,13232 mm, sedangkan IDW menghasilkan RMSE sebesar 30,78124 mm. Secara keseluruhan, grafik menunjukkan bahwa metode Nearest Neighbor memberikan performa interpolasi terbaik berdasarkan evaluasi 5-fold cross-validation pada data curah hujan yang digunakan.
Kesimpulan
Berdasarkan hasil interpolasi spasial curah hujan pada 458 stasiun pengamatan di Jawa Barat tanggal 4 Desember 2024, diperoleh lima metode interpolasi, yaitu Closest Observation, Nearest Neighbor, Inverse Distance Weighting (IDW), Ensemble, dan Ordinary Kriging. Kelima metode digunakan untuk memperkirakan nilai curah hujan pada lokasi yang tidak memiliki pengamatan langsung berdasarkan informasi dari lokasi-lokasi di sekitarnya.
Evaluasi performa metode dilakukan menggunakan 5-fold cross-validation dengan Root Mean Square Error (RMSE). Hasil evaluasi menunjukkan bahwa metode Nearest Neighbor menghasilkan RMSE terkecil sebesar 27,07995 mm, diikuti Ordinary Kriging sebesar 27,88770 mm, Ensemble sebesar 29,13232 mm, IDW sebesar 30,78124 mm, dan Closest Observation sebesar 34,32431 mm. Dengan demikian, Nearest Neighbor dengan lima tetangga terdekat merupakan metode yang memberikan performa prediksi terbaik di antara kelima metode yang diuji pada data ini.
Hasil tersebut menunjukkan bahwa pendekatan rata-rata nilai curah hujan dari lima titik pengamatan terdekat mampu menghasilkan prediksi yang relatif lebih akurat dibandingkan metode lainnya. Meskipun demikian, Ordinary Kriging juga menunjukkan performa yang kompetitif dengan nilai RMSE yang tidak jauh berbeda dari Nearest Neighbor. Oleh karena itu, pemilihan metode interpolasi tidak hanya perlu mempertimbangkan nilai RMSE, tetapi juga karakteristik data dan asumsi yang mendasari setiap metode.