1 Pendahuluan

Curah hujan merupakan salah satu unsur lingkungan yang penting untuk dikaji karena berkaitan dengan ketersediaan air, pertanian, risiko banjir, dan kekeringan. Laporan ini menganalisis curah hujan bulanan pada stasiun pengamatan BMKG untuk periode September 2026. Tujuannya adalah mendeskripsikan data, memetakan nilai pada lokasi stasiun, membandingkan hasil interpolasi Inverse Distance Weighted (IDW), Ordinary Kriging (OK), dan Universal Kriging (UK), serta menentukan model dengan kesalahan prediksi paling kecil melalui validasi silang Leave-One-Out Cross Validation (LOOCV).

Data pada berkas yang diberikan mencakup stasiun di Pulau Jawa. Karena itu, hasilnya tidak mewakili seluruh Indonesia, melainkan cakupan stasiun yang tersedia dalam berkas.

2 Metode Penelitian

2.1 Data dan variabel

Variabel utama adalah Curah Hujan (mm) untuk periode September 2026. Koordinat lintang dan bujur digunakan sebagai lokasi titik. Nilai kosong yang pada data sumber ditulis “-” diperlakukan sebagai data hilang, bukan sebagai nol.

2.2 Metode analisis

  1. Statistik deskriptif: jumlah data valid, rata-rata, median, simpangan baku, minimum, maksimum, dan kuartil.
  2. Peta titik: menampilkan curah hujan pada lokasi stasiun yang memiliki data.
  3. Interpolasi IDW: nilai di lokasi tak teramati diperkirakan sebagai rata-rata berbobot jarak; titik yang lebih dekat mendapat bobot lebih besar.
  4. Ordinary Kriging (OK): interpolasi geostatistik dengan asumsi rerata lokal konstan tetapi tidak diketahui.
  5. Universal Kriging (UK): kriging yang memasukkan tren spasial koordinat.
  6. LOOCV: satu observasi dikeluarkan bergantian, lalu diprediksi dari observasi lainnya. Model dibandingkan dengan MAE, RMSE, dan bias. MAE/RMSE lebih kecil menunjukkan prediksi lebih akurat pada data validasi.

2.3 Paket

pkgs <- c("readxl", "dplyr", "ggplot2", "sf", "gstat", "knitr")
for (p in pkgs) {
  if (!requireNamespace(p, quietly = TRUE)) install.packages(p)
}
library(readxl)
library(dplyr)
library(ggplot2)
library(sf)
library(gstat)
library(knitr)

3 Hasil dan Pembahasan

3.1 Input dan pembersihan data

Simpan berkas Curah_Hujan_BMKG_Sep2026_Gabungan (1).xlsx di folder kerja yang sama dengan berkas Rmd, atau ubah nilai file_data sesuai lokasi berkas Anda.

file_data <- "~/Documents/intan/SL/Curah_Hujan_BMKG_Sep2026_Gabungan (1).xlsx"
data_raw <- read_excel(file_data, sheet = "Data Gabungan")

# Standarisasi nama kolom agar lebih mudah dipanggil
names(data_raw) <- make.names(names(data_raw))
# Kolom hasil make.names: Nama.Stasiun, Lintang, Bujur, Curah.Hujan..mm.
data <- data_raw |>
  transmute(
    id = `ID.WMO`,
    stasiun = Nama.Stasiun,
    lintang = Lintang,
    bujur = Bujur,
    elevasi = Elevasi..m.,
    periode = Periode,
    curah_hujan = Curah.Hujan..mm.
  )

# Pastikan variabel numerik benar-benar numerik
data$curah_hujan <- as.numeric(data$curah_hujan)
data$lintang <- as.numeric(data$lintang)
data$bujur <- as.numeric(data$bujur)

cat("Jumlah stasiun pada berkas:", nrow(data), "\n")
## Jumlah stasiun pada berkas: 35
cat("Data curah hujan valid:", sum(!is.na(data$curah_hujan)), "\n")
## Data curah hujan valid: 29
cat("Data curah hujan hilang:", sum(is.na(data$curah_hujan)), "\n")
## Data curah hujan hilang: 6
data_valid <- data |> filter(!is.na(curah_hujan),
                             !is.na(lintang), !is.na(bujur))

3.2 Statistik deskriptif

x <- data_valid$curah_hujan
statistik <- data.frame(
  Ukuran = c("Jumlah data valid", "Rata-rata (mm)", "Median (mm)",
             "Simpangan baku (mm)", "Minimum (mm)", "Kuartil 1 (mm)",
             "Kuartil 3 (mm)", "Maksimum (mm)"),
  Nilai = c(length(x), mean(x), median(x), sd(x), min(x),
            unname(quantile(x, 0.25)), unname(quantile(x, 0.75)), max(x))
)
knitr::kable(statistik, digits = 2, caption = "Statistik deskriptif curah hujan")
Statistik deskriptif curah hujan
Ukuran Nilai
Jumlah data valid 29.00
Rata-rata (mm) 19.06
Median (mm) 0.40
Simpangan baku (mm) 56.03
Minimum (mm) 0.00
Kuartil 1 (mm) 0.00
Kuartil 3 (mm) 7.10
Maksimum (mm) 290.20 
cat("Terdapat", length(x), "stasiun dengan data valid dari", nrow(data),
    "stasiun pada berkas. Rata-rata curah hujan adalah",
    round(mean(x), 2), "mm, sedangkan median", round(median(x), 2),
    "mm. Perbedaan rata-rata dan median menunjukkan distribusi menceng ke kanan,",
    "yang dipengaruhi oleh nilai tinggi pada sebagian stasiun. Nilai maksimum",
    round(max(x), 2), "mm tercatat di",
    data_valid$stasiun[which.max(data_valid$curah_hujan)], ".")

Terdapat 29 stasiun dengan data valid dari 35 stasiun pada berkas. Rata-rata curah hujan adalah 19.06 mm, sedangkan median 0.4 mm. Perbedaan rata-rata dan median menunjukkan distribusi menceng ke kanan, yang dipengaruhi oleh nilai tinggi pada sebagian stasiun. Nilai maksimum 290.2 mm tercatat di Stasiun Klimatologi Jawa Barat .

3.3 Peta titik stasiun

ggplot(data_valid, aes(x = bujur, y = lintang,
                       color = curah_hujan, size = curah_hujan)) +
  geom_point(alpha = 0.85) +
  scale_color_viridis_c(option = "C", name = "Curah hujan\n(mm)") +
  scale_size_continuous(range = c(2, 8), guide = "none") +
  coord_equal() +
  labs(title = "Curah Hujan pada Stasiun BMKG",
       subtitle = "September 2026; titik menunjukkan lokasi stasiun",
       x = "Bujur (derajat)", y = "Lintang (derajat)") +
  theme_minimal()

Peta ini memperlihatkan nilai pada titik pengamatan, bukan nilai untuk seluruh area. Area di antara stasiun akan diperkirakan melalui interpolasi. Tidak tersedia batas administrasi/peta dasar pada berkas yang diberikan, sehingga visualisasi awal menggunakan koordinat stasiun.

3.4 Persiapan koordinat untuk interpolasi

Interpolasi jarak sebaiknya dilakukan pada sistem koordinat terproyeksi dalam satuan meter, bukan langsung pada derajat lintang-bujur. Stasiun berada di sekitar Pulau Jawa; UTM zona 49S digunakan sebagai pendekatan proyeksi kerja. Periksa cakupan koordinat bila data diganti.

pts_ll <- st_as_sf(data_valid, coords = c("bujur", "lintang"),
                   crs = 4326, remove = FALSE)
pts <- st_transform(pts_ll, 32749) # WGS84 / UTM zone 49S
coords <- st_coordinates(pts)
pts$x <- coords[, 1]
pts$y <- coords[, 2]

# Grid prediksi di dalam kotak pembatas lokasi stasiun
grid_xy <- expand.grid(
  x = seq(min(pts$x), max(pts$x), length.out = 100),
  y = seq(min(pts$y), max(pts$y), length.out = 100)
)
grid_sf <- st_as_sf(
  grid_xy,
  coords = c("x", "y"),
  crs = st_crs(pts),
  remove = FALSE
)

3.5 Interpolasi IDW, Ordinary Kriging, dan Universal Kriging

# IDW, power = 2
idw_fit <- gstat::idw(curah_hujan ~ 1, locations = pts,
                      newdata = grid_sf, idp = 2)
## [inverse distance weighted interpolation]
# Variogram empiris dan model teoritis untuk kriging
vg_emp <- gstat::variogram(curah_hujan ~ 1, data = pts)
vg_fit <- gstat::fit.variogram(
  vg_emp, model = gstat::vgm(psill = var(pts$curah_hujan, na.rm = TRUE) * 0.8,
                            model = "Exp", range = 50000,
                            nugget = var(pts$curah_hujan, na.rm = TRUE) * 0.2)
)
plot(vg_emp, vg_fit, main = "Variogram empiris dan model eksponensial")

ok_fit <- gstat::krige(curah_hujan ~ 1, locations = pts,
                       newdata = grid_sf, model = vg_fit)
## [using ordinary kriging]
# Universal Kriging dengan tren koordinat linear
# Universal Kriging dengan tren koordinat linear
uk_fit <- gstat::krige(
  curah_hujan ~ x + y,
  locations = pts,
  newdata = grid_sf,
  model = vg_fit
)
## [using universal kriging]
# Fungsi untuk menggambar prediksi pada grid
plot_pred <- function(obj, kolom, judul) {
  xy <- st_coordinates(obj)
  z <- obj[[kolom]]
  ggplot(data.frame(x = xy[,1], y = xy[,2], z = z),
         aes(x, y, fill = z)) +
    geom_raster() +
    geom_point(data = pts, aes(x = x, y = y),
               inherit.aes = FALSE, size = 1.2, color = "black") +
    scale_fill_viridis_c(option = "C", name = "mm") +
    coord_equal() +
    labs(title = judul, x = "UTM Easting (m)", y = "UTM Northing (m)") +
    theme_minimal()
}
plot_pred(idw_fit, "var1.pred", "Interpolasi IDW (power = 2)")

plot_pred(ok_fit, "var1.pred", "Interpolasi Ordinary Kriging")

plot_pred(uk_fit, "var1.pred", "Interpolasi Universal Kriging")

Catatan interpretasi: permukaan interpolasi hanya merupakan estimasi di dalam kotak pembatas stasiun, bukan batas wilayah administratif atau peta risiko resmi. Hasil di tepi area dan area dengan stasiun jarang perlu ditafsirkan hati-hati. Kriging juga bergantung pada pemilihan dan kecocokan model variogram.

3.6 Validasi silang Leave-One-Out Cross Validation

# IDW LOOCV manual: tiap observasi diprediksi dari stasiun lainnya
idw_cv_pred <- rep(NA_real_, nrow(pts))
for (i in seq_len(nrow(pts))) {
  train_i <- pts[-i, ]
  test_i <- pts[i, ]
  pred_i <- gstat::idw(curah_hujan ~ 1, locations = train_i,
                       newdata = test_i, idp = 2)
  idw_cv_pred[i] <- pred_i$var1.pred
}
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
# Kriging LOOCV dari gstat
ok_cv <- gstat::krige.cv(curah_hujan ~ 1, locations = pts,
                         model = vg_fit, nfold = nrow(pts))
uk_cv <- gstat::krige.cv(curah_hujan ~ x + y, locations = pts,
                         model = vg_fit, nfold = nrow(pts))

metrik <- function(obs, pred, nama) {
  e <- obs - pred
  data.frame(Model = nama,
             MAE = mean(abs(e), na.rm = TRUE),
             RMSE = sqrt(mean(e^2, na.rm = TRUE)),
             Bias = mean(e, na.rm = TRUE))
}
hasil_cv <- dplyr::bind_rows(
  metrik(pts$curah_hujan, idw_cv_pred, "IDW (p = 2)"),
  metrik(ok_cv$observed, ok_cv$var1.pred, "Ordinary Kriging"),
  metrik(uk_cv$observed, uk_cv$var1.pred, "Universal Kriging")
) |> arrange(RMSE)

knitr::kable(hasil_cv, digits = 3,
             caption = "Perbandingan akurasi berdasarkan LOOCV")
Perbandingan akurasi berdasarkan LOOCV
Model MAE RMSE Bias
Ordinary Kriging 19.980 48.279 1.344
Universal Kriging 20.926 48.528 1.541
IDW (p = 2) 18.667 50.340 2.616

Model dengan RMSE paling kecil dipilih sebagai model terbaik berdasarkan kriteria RMSE; MAE digunakan sebagai ukuran pendukung dan bias untuk melihat kecenderungan prediksi. Keputusan akhir harus mengikuti tabel hasil yang dihasilkan ketika dokumen dijalankan.

3.7 Kesimpulan

Kesimpulan akhir diisi berdasarkan keluaran RMarkdown: (1) ringkasan pola curah hujan, (2) pola spasial pada peta titik dan permukaan interpolasi, dan (3) model terbaik berdasarkan RMSE LOOCV. Jumlah data hilang dilaporkan secara eksplisit dan tidak diubah menjadi nol. Hasil interpolasi bersifat estimasi serta dibatasi oleh sebaran stasiun yang tersedia.

4 Referensi data

Badan Meteorologi, Klimatologi, dan Geofisika (BMKG). Data curah hujan stasiun periode September 2026, dihimpun pada berkas Curah_Hujan_BMKG_Sep2026_Gabungan (1).xlsx yang digunakan dalam analisis ini.