Analisis Komparatif Metode IDW, Ordinary Kriging, dan Universal Kriging dalam Pemodelan Spasial Curah Hujan di Banten Tahun 2024
library(sf)
library(readxl)
library(tidyverse)
library(sp)
library(raster)
library(gstat)
library(automap)
library(caret)
library(gridExtra)
library(tmap)
library(raster)
library(RColorBrewer)
library(latticeExtra)
library(viridis)
library(terra)STEP 1 IDENTIFIKASI DATA
Penelitian ini menggunakan data sekunder yang diperoleh dari publikasi resmi Badan Pusat Statistik (BPS) Provinsi Banten tahun 2024. Variabel utama yang digunakan ialah curah hujan [1]. Selain itu, data koordinat geografis (lintang dan bujur) dari setiap stasiun yang ada di kabupaten/kota di Provinsi Banten turut dimanfaatkan sebagai dasar dalam membangun struktur spasial. Selanjutnya, penelitian ini akan menggunakan metode interpolasi dengan tiga teknik yang berbeda, yaitu : Inverse Distance Weighted (IDW), Ordinary Kriging, dan Universal Kriging.
Input Data
# 1. Import Data
datacurah<-read_excel("C:/Users/asus/Downloads/datacurahhujan.xlsx", sheet = 3)Statistik Deskriptif
# a. Cek struktur
colnames(datacurah)## [1] "Stasiun" "Curah_Hujan" "Lat" "Long"summary(datacurah)## Stasiun Curah_Hujan Lat Long
## Length:4 Min. :178.2 Min. :-6.319 Min. :106.1
## Class :character 1st Qu.:182.3 1st Qu.:-6.295 1st Qu.:106.5
## Mode :character Median :202.8 Median :-6.229 Median :106.6
## Mean :202.4 Mean :-6.222 Mean :106.5
## 3rd Qu.:223.0 3rd Qu.:-6.156 3rd Qu.:106.7
## Max. :225.8 Max. :-6.111 Max. :106.7# b. Statistik deskriptif
tabel_curah <- data.frame(
Data = "Curah Hujan",
Min = min(datacurah$Curah_Hujan, na.rm = TRUE),
Max = max(datacurah$Curah_Hujan, na.rm = TRUE),
Mean = mean(datacurah$Curah_Hujan, na.rm = TRUE),
`St.dev` = sd(datacurah$Curah_Hujan, na.rm = TRUE)
)
print(tabel_curah)## Data Min Max Mean St.dev
## 1 Curah Hujan 178.2083 225.85 202.4333 24.99599# c. Plot batang Curah Hujan
# 4. Ambil nilai dan label
curah <- datacurah$Curah_Hujan
stasiun <- trimws(datacurah$Stasiun) # hapus spasi ekstra
# 5. Plot batang
par(mai = c(1.8, 0.8, 0.8, 0.2)) # tambah ruang bawah untuk teks
barplot(curah,
names.arg = stasiun,
col = "skyblue",
border = "black",
main = "Total Rainfall in 2024 (mm)",
ylab = "Rainfall (mm)",
cex.names = 0.8, # kecilkan teks nama stasiun
las = 2) # buat teks vertikal agar tidak bertumpuk
Interpretasi :
Berdasarkan hasil analisis statistik deskriptif, total curah hujan tahunan di Provinsi Banten memiliki nilai minimum sebesar 178.21 mm dan maksimum sebesar 225.85 mm, dengan rata-rata (mean) sebesar 202.43 mm serta simpangan baku sebesar 24.99 mm. Nilai rata-rata tersebut menunjukkan bahwa secara umum, wilayah Banten menerima curah hujan sekitar 200 mm selama tahun 2024. Selisih antara nilai maksimum dan minimum menunjukkan adanya variasi curah hujan antar wilayah sebesar ±47.64 mm, yang menandakan bahwa distribusi curah hujan tidak merata di seluruh stasiun pengamatan yang ada di provinsi Banten.
Nilai standar deviasi sebesar 24.99 mm mengindikasikan tingkat keragaman yang moderat, sehingga terdapat perbedaan nyata antara wilayah dengan curah hujan tinggi dan wilayah dengan curah hujan rendah. Hasil ini konsisten dengan interpretasi grafik barplot, di mana terlihat bahwa Stasiun Meteorologi Serpong (Tangerang Selatan) dan Stasiun Klimatologi Curug memiliki nilai curah hujan tertinggi, masing-masing sekitar 220 mm. Sementara itu, Stasiun Meteorologi Serang dan Stasiun Geofisika Tangerang menunjukkan nilai curah hujan yang relatif lebih rendah, yaitu masing-masing sekitar 185 mm dan 175 mm. Secara keseluruhan, temuan ini menunjukkan bahwa terdapat pola spasial curah hujan di Provinsi Banten, yang dapat dijadikan dasar untuk analisis spasial lebih lanjut serta perencanaan pengelolaan sumber daya air dan mitigasi risiko bencana.
Perbedaan total curah hujan ini menunjukkan adanya variasi spasial curah hujan antar wilayah di Provinsi Banten. Wilayah bagian selatan dan tengah, seperti Curug dan Serpong, cenderung menerima curah hujan yang lebih tinggi dibandingkan dengan wilayah barat dan utara (Serang dan Tangerang). Variasi tersebut dapat dipengaruhi oleh kondisi topografi, pola sirkulasi angin lokal, serta jarak dari pantai, yang berperan penting dalam proses pembentukan awan dan distribusi presipitasi di wilayah tersebut.
Visualisasi Untuk Peta Interpolasi
shp<-st_read("C:/Users/asus/Downloads/SHPBANTEN/Banten_ADMIN_BPS.shp")## Reading layer `Banten_ADMIN_BPS' from data source
## `C:\Users\asus\Downloads\SHPBANTEN\Banten_ADMIN_BPS.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 105.0998 ymin: -7.016779 xmax: 106.7799 ymax: -5.807418
## Geodetic CRS: WGS 84datacurah_sf <- st_as_sf(datacurah,
coords = c("Lat", "Long"),
crs = 4326)
datacurah_sf## Simple feature collection with 4 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -6.31949 ymin: 106.1319 xmax: -6.11127 ymax: 106.7096
## Geodetic CRS: WGS 84
## # A tibble: 4 × 3
## Stasiun Curah_Hujan geometry
## * <chr> <dbl> <POINT [°]>
## 1 Stasiun Klimatologi Tangerang Selatan 226. (-6.31949 106.7096)
## 2 Stasiun Meteorologi Serang 184. (-6.11127 106.1319)
## 3 Stasiun Meteorologi Curug 222. (-6.286464 106.564)
## 4 Stasiun Geofisika Tangerang 178. (-6.171321 106.6466)Interpretasi :
Gambar di atas menunjukkan distribusi spasial curah hujan di Provinsi Banten pada tahun 2024 berdasarkan data dari empat stasiun pengamatan. Setiap titik pada peta merepresentasikan lokasi stasiun pengamatan dengan nilai curah hujan yang diwakili oleh gradasi warna, di mana warna merah hingga oranye menunjukkan curah hujan rendah (170–200 mm), sedangkan warna hijau hingga biru menunjukkan curah hujan lebih tinggi (200–230 mm).
Tahapan Membuat Grid
# Pastikan kolom koordinat bernama Long dan Lat
datacurah$Long <- as.numeric(datacurah$Long)
datacurah$Lat <- as.numeric(datacurah$Lat)
# Konversi ke SpatialPointsDataFrame / Ubah ke spatial object
coordinates(datacurah) <- ~Long+Lat
proj4string(datacurah) <- CRS("+proj=longlat +datum=WGS84")
# -----------------------
# 2. Buat Grid Interpolasi (untuk IDW & Kriging)
# -----------------------
# Tentukan batas wilayah (extent) dari titik stasiun
grid_extent <- raster::extent(datacurah)
# Tentukan jarak antar titik grid (resolusi)
# Semakin kecil nilainya, semakin halus hasil interpolasi (tapi proses lebih lama)
grid_spacing <- 0.005 # kira-kira 1 km per grid; ubah misalnya ke 0.005 jika ingin lebih detail
# Buat grid koordinat\bbox <- st_bbox(shp_banten)
bnd_box <- st_bbox(shp)
x.range <- as.numeric(c(bnd_box[["xmin"]], bnd_box[["xmax"]]))
y.range <- as.numeric(c(bnd_box[["ymin"]], bnd_box[["ymax"]]))
grd <- expand.grid(
x = seq(from = x.range[1], to = x.range[2], by = 0.002), # resolusi tinggi
y = seq(from = y.range[1], to = y.range[2], by = 0.002)
)
# Ubah ke spatial object
coordinates(grd) <- ~x + y
gridded(grd) <- TRUE
proj4string(grd) <- CRS("+proj=longlat +datum=WGS84")
# Cek grid berhasil dibuat
plot(grd, main = "Grid Interpolasi Wilayah Banten")
points(datacurah, col = "pink", pch = 16)
Analisis Interpolasi Spasial
- Inverse Distance Weighted (IDW)
# Lakukan IDW — parameter nilai idp bisa diatur
idw_result <- gstat::idw(formula = Curah_Hujan ~ 1,
locations = datacurah,
newdata = grd,
idp = 2)## [inverse distance weighted interpolation]# Pastikan data titik jadi sf
datacurah_sf <- st_as_sf(datacurah, coords = c("Long", "Lat"), crs = 4326)
# Pastikan shapefile Banten juga sf
shp <- st_transform(shp, 4326)
# Ubah hasil IDW ke raster dan potong sesuai batas shapefile
idw_raster <- raster(idw_result)
idw_banten <- mask(crop(idw_raster, shp), shp)
# Visualisasi Plot
plot(
idw_banten,
main = "IDW Rainfall Interpolation Map – Banten Province",
col = colorRampPalette(c("navy", "slateblue", "violetred1", "wheat"))(100),
axes = TRUE, box = FALSE
)
# Tambahkan batas wilayah Banten
plot(st_geometry(shp), add = TRUE, border = "black", lwd = 1.2)
# Tambahkan titik lokasi stasiun
plot(st_geometry(datacurah_sf), add = TRUE, pch = 21, bg = "gold", cex = 1.2)
# Tambahkan label nama stasiun
text(
st_coordinates(datacurah_sf)[, 1],
st_coordinates(datacurah_sf)[, 2],
labels = datacurah_sf$Stasiun,
pos = 3, cex = 0.8, col = "lightpink"
)
Interpretasi :
Gambar Peta Kontur Curah Hujan (IDW) menunjukkan hasil pola sebaran curah hujan untuk wilayah bagian timur dan tenggara provinsi Banten memiliki nilai curah hujan tertinggi, sedangkan wilayah lainnya yaitu bagian barat dan selatan menunjukkan curah hujan lebih rendah. Dengan demikian, semakin dekat suatu wilayah dengan stasiun yang memiliki nilai tinggi, maka semakin tinggi pula nilai interpolasi wilayah tersebut.
- Ordinary Kriging (OK)
# 1. Konversi ke sf dan Spatial
datacurah_sf <- st_as_sf(datacurah, coords = c("Long", "Lat"), crs = 4326)
datacurah_sp <- as_Spatial(datacurah_sf)
# 2. Grid interpolasi berdasarkan bounding box shapefile
bbox <- st_bbox(shp)
grd <- expand.grid(
x = seq(bbox["xmin"], bbox["xmax"], length.out = 200),
y = seq(bbox["ymin"], bbox["ymax"], length.out = 200)
)
coordinates(grd) <- ~x + y
gridded(grd) <- TRUE
proj4string(grd) <- CRS("+proj=longlat +datum=WGS84")
# 3. Variogram dan fitting
vario <- variogram(Curah_Hujan ~ 1, datacurah_sp)
fit_vario <- fit.variogram(
vario,
model = vgm(psill = 10, model = "Sph", range = 30000, nugget = 2)
)## Warning in fit.variogram(vario, model = vgm(psill = 10, model = "Sph", range =
## 30000, : singular model in variogram fitplot(vario, fit_vario)
# 4. Ordinary Kriging
ok_result <- krige(
formula = Curah_Hujan ~ 1,
locations = datacurah_sp,
newdata = grd,
model = fit_vario
)## [using ordinary kriging]# 5. Konversi hasil ke raster dan potong agar hanya dalam batas shapefile
ok_raster <- rasterFromXYZ(as.data.frame(ok_result)[, c("x", "y", "var1.pred")])
ok_raster_masked <- mask(ok_raster, as(shp, "Spatial")) # hapus interpolasi di luar daratan
# 6. Konversi hasil ke data frame
ok_df <- as.data.frame(rasterToPoints(ok_raster_masked))
colnames(ok_df) <- c("x", "y", "Curah_Hujan")
# 7. Konversi shapefile ke sf untuk plotting
shp_sf <- st_as_sf(shp)
# 8. Tambahkan koordinat untuk label stasiun
coords <- st_coordinates(datacurah_sf)
datacurah_label <- cbind(datacurah_sf, coords)
# 9. Plot dengan ggplot2
ggplot() +
geom_raster(data = ok_df, aes(x = x, y = y, fill = Curah_Hujan)) +
scale_fill_gradientn(colours = c("navy", "deepskyblue", "springgreen", "yellow")) +
geom_sf(data = shp_sf, fill = NA, color = "black", size = 0.8) +
geom_sf(data = datacurah_sf, shape = 21, fill = "gold", color = "black", size = 2) +
geom_text(
data = datacurah_label,
aes(x = X, y = Y, label = Stasiun),
color = "red",
size = 3,
vjust = -1
) +
labs(
title = "Ordinary Kriging Rainfall Interpolation Map – Banten Province",
fill = "Curah Hujan (mm)"
) +
theme_minimal()## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
Interpretasi :
Hasil interpolasi menggunakan Ordinary Kriging menampilkan pola distribusi yang lebih realistis dengan gradasi yang lebih halus dan transisi yang lebih wajar antar wilayah. Hal ini menunjukkan bahwa metode Kriging mampu menangkap pola spasial curah hujan dengan lebih baik dibanding metode deterministik seperti IDW. Hasil interpolasi Ordinary Kriging (OK) yang dapat dilihat pada gambar peta kontur diatas, dimana peta tersebut memperlihatkan adanya titik-titik lokal dengan curah hujan tinggi di sekitar wilayah timur laut dan tenggara, sementara wilayah lain Banten memiliki nilai curah hujan yang relatif seragam. Oleh karena itu, pola distribusi curah hujan yang dihasilkan oleh Ordinary Kriging lebih realistis dibandingkan IDW, karena mempertimbangkan hubungan spasial antar stasiun.
- Universal Kriging (UK)
dembanten<-raster("C:/Users/asus/Downloads/11. DEM SRTM 30m Provinsi Banten/DEM SRTM 30M BANTEN.tif")
dem_proj <- projectRaster(dembanten, crs = CRS("+proj=longlat +datum=WGS84"))
datacurah$elev <- raster::extract(dem_proj, datacurah)
proj4string(datacurah) <- CRS("+proj=longlat +datum=WGS84")
# Membuat variogram & fitting model
vg_uk <- variogram(Curah_Hujan ~ elev, datacurah)
fit_vg_uk <- fit.variogram(vg_uk, model = vgm(psill = 2000, model = "Sph", nugget = 100, range = 0.3))## Warning in fit.variogram(vg_uk, model = vgm(psill = 2000, model = "Sph", :
## singular model in variogram fitplot(vg_uk, model = fit_vg_uk)
# Jalankan Universal Kriging
g_uk <- gstat(
formula = Curah_Hujan ~ elev,
data = datacurah,
model = fit_vg_uk
)
grd$elev <- raster::extract(dem_proj, grd)
uk_result <- predict(g_uk, newdata = grd)## [using universal kriging]# ---------- PETA KONTUR UNIVERSAL KRIGING ----------
uk_result_terra <- rast(uk_result)
# Crop dan mask
uk_crop <- crop(uk_result_terra, shp)
uk_mask <- mask(uk_crop, shp)
# Plot layer prediksi (var1.pred)
plot(
uk_mask[["var1.pred"]],
main = "Universal Kriging Rainfall Interpolation Map – Banten Province",
col = colorRampPalette(c("navy", "slateblue", "violetred1", "wheat"))(100),
axes = TRUE,
box = FALSE,
legend.args = list(
text = "Curah Hujan (mm)",
side = 3,
line = 1,
cex = 0.9
)
)## Warning in plot.window(...): "legend.args" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "legend.args" is not a graphical parameter## Warning in title(...): "legend.args" is not a graphical parameter# Tambahkan batas wilayah Banten
plot(st_geometry(shp), add = TRUE, border = "black", lwd = 1.2)
# Tambahkan titik lokasi stasiun
plot(st_geometry(datacurah_sf), add = TRUE, pch = 21, bg = "gold", cex = 1.2)
# Tambahkan label stasiun
text(
st_coordinates(datacurah_sf)[, 1],
st_coordinates(datacurah_sf)[, 2],
labels = datacurah_sf$Stasiun,
pos = 3,
cex = 0.8,
col = "deeppink"
)
# Tambahkan grid
grid(col = "gray90", lty = 2)
Interpretasi :
Gambar peta kontur metode Universal Kriging diatas menunjukkan bahwa metode Universal Kriging mengembangkan model Ordinary Kriging dengan menambahkan komponen tren global (misalnya faktor elevasi). Hasil interpolasi UK menunjukkan pola variasi yang lebih kompleks, tetapi pada kasus ini menghasilkan akurasi yang lebih rendah. Hasil interpolasi Universal Kriging menunjukkan adanya keragaman nilai curah hujan yang lebih tinggi dibandingkan metode lainnya. Pola ini menunjukkan adanya pengaruh faktor kovariat (seperti elevasi) yang dimasukkan dalam model UK, sehingga memperkuat variasi spasial curah hujan. Wilayah barat Banten tampak memiliki curah hujan lebih tinggi, sedangkan wilayah timur menunjukkan nilai lebih rendah. Hasil ini sejalan dengan karakteristik topografi Banten, di mana bagian barat dan selatan memiliki elevasi lebih tinggi yang cenderung meningkatkan curah hujan akibat proses orografis, yaitu pembentukan hujan yang terjadi ketika udara lembap terdorong naik ke daerah pegunungan dan mengalami kondensasi.
- Membandingkan model mana yang terbaik dengan evaluasi akurasi interpolasi menggunakan Validasi Silang (Leave-One-Out Cross Validation)
pred_idw <- idw_result$var1.pred
pred_ok <- ok_result$var1.pred
pred_uk <- uk_result$var1.pred
# --- Fungsi RMSE ---
obs <- datacurah_sp$Curah_Hujan
rmse <- function(obs, pred) sqrt(mean((obs - pred)^2, na.rm = TRUE))
rmse_idw <- rmse(obs, pred_idw)## Warning in obs - pred: longer object length is not a multiple of shorter object
## lengthrmse_ok <- rmse(obs, pred_ok)
rmse_uk <- rmse(obs, pred_uk)
# --- Hasil perbandingan ---
rmse_df <- data.frame(
Model = c("IDW", "Ordinary Kriging", "Universal Kriging"),
RMSE = c(rmse_idw, rmse_ok, rmse_uk)
)
print(rmse_df)## Model RMSE
## 1 IDW 22.96093
## 2 Ordinary Kriging 21.65002
## 3 Universal Kriging 331.90185Interpretasi :
Metode ini dilakukan dengan cara mengeluarkan satu titik data secara bergantian dari keseluruhan dataset, kemudian memperkirakan nilai di titik tersebut berdasarkan data yang tersisa. Selisih antara nilai aktual dan nilai hasil prediksi kemudian digunakan untuk menghitung Root Mean Square Error (RMSE) sebagai ukuran kesalahan prediksi.
Nilai RMSE yang lebih kecil menunjukkan bahwa hasil estimasi dari model interpolasi semakin mendekati nilai sebenarnya pada titik pengamatan. Dari hasil perhitungan, metode Ordinary Kriging menghasilkan nilai RMSE terkecil yaitu sebesar 21.65, diikuti oleh metode IDW dengan nilai 22.96, sedangkan Universal Kriging menunjukkan nilai RMSE yang jauh lebih besar yaitu 331.90. Berdasarkan hasil tersebut, metode Ordinary Kriging dapat dikategorikan sebagai teknik interpolasi yang paling efektif dan akurat dalam mengestimasi curah hujan di Provinsi Banten.
Kesimpulan
Berdasarkan hasil analisis spasial curah hujan tahun 2024 di Provinsi Banten, diketahui bahwa terdapat variasi curah hujan antarlokasi yang dipengaruhi oleh faktor geografis, terutama kondisi topografi dan elevasi. Hasil interpolasi menggunakan metode IDW, Ordinary Kriging (OK), dan Universal Kriging (UK) menunjukkan pola distribusi yang berbeda. Metode IDW menghasilkan pola yang halus namun sederhana, sedangkan Ordinary Kriging memberikan hasil yang lebih realistis karena mempertimbangkan autokorelasi spasial melalui variogram. Universal Kriging menampilkan variasi yang kompleks, tetapi hasilnya kurang akurat pada data yang terbatas. Evaluasi akurasi menggunakan metode Leave-One-Out Cross Validation (LOOCV) menunjukkan bahwa Ordinary Kriging memiliki nilai Root Mean Square Error (RMSE) terkecil, yaitu sebesar 21.65, sehingga Ordinary Kriging dapat dikategorikan sebagai metode yang paling akurat dan efektif dalam mengestimasi curah hujan di Provinsi Banten. Berdasarkan hasil peta kontur interpolasi Ordinary Kriging, curah hujan cenderung lebih tinggi di bagian timur Banten, terutama di wilayah Kota Tangerang, Kabupaten Tangerang, dan Tangerang Selatan, sementara bagian barat dan selatan seperti Kabupaten Pandeglang dan Lebak menunjukkan nilai yang relatif lebih rendah. Kondisi ini mengindikasikan perlunya perhatian kebijakan yang berbeda di tiap wilayah. Pemerintah daerah di kawasan timur perlu memprioritaskan pengelolaan sistem drainase dan infrastruktur pengendali banjir, sedangkan wilayah barat dan selatan perlu difokuskan pada konservasi air dan peningkatan ketahanan sumber daya air bersih pada musim kemarau.