Perbandingan Metode Interpolasi Spasial IDW, Ordinary Kriging, dan Universal Kriging dalam Pemetaan Curah Hujan Tahunan di Indonesia Tahun 2024

1 Pendahuluan

1.1 Latar Belakang

Curah hujan merupakan salah satu unsur iklim yang berperan penting dalam pengelolaan sumber daya air, sektor pertanian, serta mitigasi bencana hidrometeorologi seperti banjir dan kekeringan. Informasi mengenai pola dan sebaran curah hujan diperlukan untuk memahami variasi kondisi lingkungan antarwilayah dan mendukung pengambilan keputusan yang berkaitan dengan pengelolaan air serta perencanaan wilayah. Di Indonesia, kondisi geografis dan iklim yang beragam menyebabkan curah hujan berbeda antarwilayah. Variasi tersebut antara lain dipengaruhi oleh angin monsun dan kondisi topografi. Sebagai contoh, curah hujan maksimum di sekitar Gunung Merapi secara konsisten terpusat di sisi barat, yang dikaitkan dengan angin monsun barat dan pengangkatan massa udara di lereng barat (Sofia et al., 2018).

Pengukuran curah hujan umumnya dilakukan melalui stasiun pengamatan yang tersebar pada lokasi tertentu. Namun, nilai yang diperoleh dari stasiun tersebut belum dapat memberikan informasi secara langsung mengenai kondisi curah hujan pada seluruh wilayah, khususnya pada lokasi yang tidak memiliki stasiun pengamatan. Oleh karena itu, diperlukan metode interpolasi spasial untuk memperkirakan nilai curah hujan pada lokasi yang tidak terukur berdasarkan informasi dari stasiun-stasiun yang tersedia.

Interpolasi spasial dapat dilakukan menggunakan beberapa metode dengan pendekatan yang berbeda. Inverse Distance Weighted (IDW) memperkirakan nilai pada lokasi tertentu berdasarkan nilai pengamatan di sekitarnya dengan memberikan bobot yang dipengaruhi oleh jarak. Sementara itu, Ordinary Kriging (OK) memanfaatkan struktur ketergantungan spasial yang dimodelkan melalui variogram dengan asumsi rata-rata lokal konstan tetapi tidak diketahui. Adapun Universal Kriging (UK) mempertimbangkan adanya tren spasial dalam proses pendugaan. Perbedaan pendekatan tersebut memungkinkan setiap metode menghasilkan estimasi dengan tingkat ketelitian yang berbeda. Pada data curah hujan bulanan Selandia Baru dari 294 pos hujan, metode geostatistik menghasilkan galat yang lebih kecil dibandingkan IDW pada seluruh bulan, sedangkan penambahan elevasi dan koordinat sebagai informasi tambahan hanya memberikan peningkatan kecil (Caloiero et al., 2021). Hal ini menunjukkan bahwa kinerja metode interpolasi perlu dibandingkan pada data yang dianalisis, bukan diasumsikan.

Penelitian ini menggunakan data curah hujan tahunan stasiun pengamatan di Indonesia dari Southeast Asia Climate Assessment & Dataset (SACA&D) yang dikelola BMKG, bagian dari International Climate Assessment & Dataset (ICA&D) (van den Besselaar et al., 2015). Analisis mencakup pemeriksaan dan pembersihan data, eksplorasi karakteristik curah hujan, visualisasi sebaran stasiun, serta pemetaan hasil interpolasi menggunakan metode IDW, Ordinary Kriging, dan Universal Kriging. Kinerja ketiga metode dibandingkan menggunakan Leave-One-Out Cross-Validation (LOOCV), yaitu prosedur validasi silang yang dilakukan dengan mengeluarkan satu pengamatan secara bergantian dari proses pemodelan, kemudian membandingkan nilai prediksi dengan nilai pengamatan sebenarnya. Melalui analisis tersebut, diharapkan diperoleh gambaran mengenai pola sebaran curah hujan tahunan di Indonesia serta informasi mengenai perbedaan kinerja ketiga metode interpolasi.

1.2 Penelitian Terdahulu dan Celah Penelitian

Interpolasi curah hujan di Indonesia telah banyak dikaji pada skala lokal dan regional. Sanusi et al. (2024) membandingkan IDW (pangkat 2) dengan Ordinary Kriging berbagai model semivariogram (spherical, gaussian, exponential) pada curah hujan harian maksimum dari 56 pos hujan di Sulawesi Selatan periode 1986–2021. Ordinary Kriging dengan semivariogram gaussian memberikan galat terkecil (RMSE 57,45; MAE 46,49), sedikit lebih baik daripada IDW (RMSE 59,96; MAE 48,10). Universal Kriging dalam penelitian tersebut hanya dibahas pada tinjauan pustaka dan tidak diuji. Tahitu et al. (2026) menerapkan IDW untuk curah hujan tahun 2023 di Kabupaten Langkat berdasarkan 25 titik pengamatan, dengan variasi pangkat jarak 1 sampai 6,8. Pangkat 6,8 memberikan RMSE terkecil (49,99 mm/tahun), dengan evaluasi menggunakan lima titik uji dan tanpa pembanding metode kriging.

Penelitian lain berfokus pada Ordinary Kriging. Bahtiyar et al. (2014) mengestimasi curah hujan bulan Februari di Kota Semarang dengan membandingkan model variogram exponential, gaussian, dan spherical. Rohma (2022) menerapkan Ordinary Kriging dengan teknik jackknife pada curah hujan bulanan Malang Raya tahun 2016 dan menyimpulkan bahwa metode tersebut dapat digunakan baik pada data berdistribusi normal maupun tidak normal. Najwa et al. (2025) mengembangkan SAR-X Kriging, yaitu gabungan model spatial autoregressive dengan variabel eksogen dan Ordinary Kriging, untuk prediksi curah hujan di Pulau Jawa. Galat MAPE SAR-X Kriging (3,82%) hanya sedikit lebih kecil daripada Ordinary Kriging (3,86%).

Dari penelitian-penelitian tersebut terdapat beberapa celah yang akan dibahas

1. Cakupan wilayah dan jumlah stasiun. Penelitian terdahulu pada umumnya mencakup satu kabupaten, kota, atau provinsi dengan puluhan pos hujan, sedangkan penelitian ini mencakup seluruh Indonesia dengan ribuan stasiun.

2. Perbandingan tiga metode. Penelitian terdahulu yang ditelaah membandingkan IDW dengan Ordinary Kriging (Sanusi et al., 2024) atau hanya memakai salah satu metode (Tahitu et al., 2026; Bahtiyar et al., 2014; Rohma, 2022). Penelitian ini membandingkan IDW, Ordinary Kriging, dan Universal Kriging secara langsung pada data yang sama.

3. Skema validasi. Evaluasi pada penelitian terdahulu menggunakan jumlah titik uji yang kecil (Tahitu et al., 2026) atau tidak dijelaskan sebagai validasi silang (Sanusi et al., 2024). Penelitian ini mengevaluasi seluruh stasiun dengan LOOCV sehingga setiap stasiun menjadi titik uji.

1.3 Rumusan Masalah

Berdasarkan latar belakang tersebut, rumusan masalah dalam penelitian ini adalah sebagai berikut.

1. Bagaimana karakteristik dan sebaran spasial curah hujan tahunan di Indonesia berdasarkan data stasiun pengamatan?

2. Bagaimana hasil pemetaan curah hujan tahunan menggunakan metode IDW, Ordinary Kriging, dan Universal Kriging?

3. Bagaimana perbandingan kinerja ketiga metode interpolasi berdasarkan hasil LOOCV, dan metode manakah yang memberikan kinerja prediksi terbaik?

1.4 Tujuan

Penelitian ini bertujuan untuk:

1. Mendeskripsikan karakteristik data dan sebaran spasial curah hujan tahunan di Indonesia.

2. Menghasilkan peta estimasi curah hujan tahunan menggunakan metode IDW, Ordinary Kriging, dan Universal Kriging.

3. Mengevaluasi dan membandingkan kinerja ketiga metode interpolasi menggunakan LOOCV, serta menentukan metode dengan kinerja prediksi terbaik berdasarkan ukuran kesalahan prediksi yang digunakan.

1.5 Manfaat Penelitian

Penelitian ini diharapkan dapat memberikan manfaat berupa:

Manfaat metodologis, yaitu memberikan pemahaman mengenai penerapan dan perbandingan metode interpolasi spasial untuk data curah hujan.

1.6 Batasan Penelitian

Penelitian ini dibatasi pada data curah hujan tahunan tahun 2024 dari stasiun pengamatan di Indonesia yang memenuhi kriteria kelengkapan data dan validasi koordinat. Metode interpolasi yang dibandingkan adalah IDW, Ordinary Kriging, dan Universal Kriging. Evaluasi kinerja prediksi dilakukan menggunakan LOOCV berdasarkan ukuran kesalahan prediksi yang diterapkan dalam analisis. Hasil penelitian merepresentasikan estimasi berdasarkan data dan model yang digunakan, sehingga interpretasinya perlu mempertimbangkan persebaran stasiun pengamatan serta keterbatasan metode interpolasi.

2 Metode Penelitian

2.1 Jenis dan Pendekatan Penelitian

Penelitian ini merupakan penelitian kuantitatif dengan pendekatan analisis deskriptif dan spasial. Analisis deskriptif digunakan untuk menggambarkan karakteristik data curah hujan tahunan, sedangkan analisis spasial dilakukan untuk memetakan estimasi curah hujan pada lokasi yang tidak memiliki pengamatan langsung. Metode interpolasi yang digunakan meliputi Inverse Distance Weighted (IDW), Ordinary Krigin* (OK), dan Universal Kriging (UK). Hasil ketiga metode dibandingkan menggunakan Leave-One-Out Cross-Validation (LOOCV). Alur analisis secara berurutan adalah persiapan data, kontrol kualitas data, analisis deskriptif, visualisasi spasial, interpolasi, dan evaluasi model (Gambar 2.1).

Gambar 2.1. Alur penelitian

2.2 Data dan Sumber Data

Data yang digunakan adalah data curah hujan tahunan (annual precipitation) tahun 2024 dari stasiun pengamatan di Indonesia. Data bersumber dari Southeast Asia Climate Assessment & Dataset (SACA&D) yang dikelola oleh BMKG dan merupakan bagian dari International Climate Assessment & Dataset (ICA&D) (van den Besselaar et al., 2015). Dua jenis berkas dari sumber tersebut digunakan, yaitu stations.txt sebagai sumber identitas, koordinat, dan ketinggian stasiun, serta berkas indeks curah hujan indexRR sebagai sumber nilai curah hujan.

Variabel utama yang dianalisis adalah ANNUAL, yaitu jumlah curah hujan tahunan. Pada berkas sumber, nilai indeks curah hujan (RR) dinyatakan dalam satuan 0,01 mm sehingga seluruh nilai ANNUAL dibagi 100 agar dinyatakan dalam milimeter (mm). Koordinat lintang (LAT) dan bujur (LON) digunakan untuk menentukan posisi stasiun. Data difokuskan pada tahun 2024 dan selanjutnya diseleksi berdasarkan ketersediaan nilai curah hujan serta validitas koordinat.

Pada tahun 2024 tercatat 7.472 stasiun Indonesia dengan data indeks curah hujan, tetapi 2.864 di antaranya tidak memiliki nilai ANNUAL. Jumlah stasiun akhir yang digunakan dalam analisis, setelah pembersihan data dan pemeriksaan koordinat, dilaporkan pada subbab Kontrol Kualitas Data (Bab 3). Jumlah tersebut merupakan stasiun yang memenuhi kriteria pengolahan data, bukan jumlah seluruh stasiun yang tersedia pada sumber data.

2.3 Tahap Pengolahan Data

2.3.1 Pengumpulan dan Penggabungan Data

Data curah hujan dibaca dari berkas indexRR, kemudian dipilih berdasarkan identitas stasiun yang berada di wilayah Indonesia. Data tersebut digabungkan dengan informasi pada berkas stations.txt menggunakan identitas stasiun (STAID) sebagai kunci penggabungan. Selanjutnya, data diseleksi untuk memperoleh pengamatan tahun 2024.

2.3.2 Pembersihan dan Validasi Data

Pembersihan data dilakukan dalam urutan berikut.

1. Kelengkapan nilai. Stasiun tanpa nilai ANNUAL dikeluarkan dari data analisis.

2. Duplikasi identitas. Pengamatan dengan STAID ganda diperiksa; hanya satu baris per stasiun yang dipertahankan.

3. Konversi koordinat. Koordinat lintang dan bujur yang semula ditulis dalam format derajat, menit, dan detik dikonversi ke derajat desimal.

4. Pemeriksaan koordinat. Koordinat diperiksa dengan tiga kriteria: (a) berada di luar kotak batas Indonesia (lintang −11° sampai 6° dan bujur 95° sampai 141°); (b) nilai lintang sama dengan nilai bujur, yang mengindikasikan kesalahan input; dan (c) berjarak lebih dari 50 km dari daratan Indonesia berdasarkan poligon batas negara.

5. Koreksi tanda lintang. Pemeriksaan menunjukkan adanya stasiun yang tanda lintangnya tertulis utara (+) pada metadata sumber, padahal nama dan bujurnya menunjukkan lokasi di selatan khatulistiwa. Koreksi dilakukan dengan membalik tanda lintang hanya jika memenuhi aturan eksplisit: lintang melebihi batas atas kotak batas dengan kebalikannya masuk rentang Indonesia dan bujur wajar, atau stasiun jatuh di luar daratan Indonesia sedangkan titik dengan tanda lintang dibalik berada dalam 25 km dari daratan. Seluruh stasiun yang dikoreksi didaftarkan agar dapat diverifikasi.

6. Stasiun yang dikeluarkan. Stasiun yang tetap gagal pemeriksaan koordinat setelah koreksi dikeluarkan dari analisis dan dicatat sebagai temuan kualitas data.

7. Koordinat kembar. Stasiun dengan koordinat yang sama dalam proyeksi analisis dikurangi menjadi satu titik, karena koordinat kembar membuat matriks kovariansi pada kriging menjadi singular.

8. Pemeriksaan nilai ekstrem. Nilai ANNUAL yang berada di luar batas pagar dalam metode rentang antarkuartil (interquartile range, IQR), yaitu di bawah \(Q_1 - 1{,}5\times\mathrm{IQR}\) atau di atas \(Q_3 + 1{,}5\times\mathrm{IQR}\), diidentifikasi. Nilai tersebut tidak langsung dibuang karena curah hujan ekstrem di Indonesia dapat bersifat nyata. Stasiun hanya dikeluarkan apabila terbukti salah ukur atau salah input berdasarkan informasi sumber.

Data yang tampak mencurigakan tidak dikoreksi hanya berdasarkan nama stasiun, melainkan berdasarkan aturan di atas dan dilaporkan agar dapat diperiksa terhadap informasi sumber. Data bersih disimpan dan digunakan pada seluruh analisis berikutnya. Rangkaian pemeriksaan koordinat diringkas pada Gambar 2.2.

Gambar 2.2. Alur kontrol kualitas koordinat dan nilai

2.3.3 Analisis Deskriptif dan Eksplorasi Spasial

Analisis deskriptif dilakukan untuk mengetahui karakteristik curah hujan tahunan melalui ukuran statistik berupa rata-rata, median, simpangan baku, koefisien variasi, nilai minimum dan maksimum, kuartil, serta kemencengan (skewness). Normalitas distribusi diuji dengan uji Shapiro-Wilk, dan distribusi data diamati melalui histogram dan boxplot.

Eksplorasi spasial dilakukan dengan memetakan lokasi stasiun berdasarkan nilai curah hujan, baik dengan skala warna kontinu maupun kelas kuantil. Untuk melihat adanya tren spasial berskala luas, dilakukan regresi linear dengan curah hujan tahunan (ANNUAL, dalam mm) sebagai variabel respons dan koordinat \(X\) (arah timur-barat, easting) serta \(Y\) (arah utara-selatan, northing) sebagai variabel penjelas, yaitu \(\mathrm{ANNUAL}=\beta_0+\beta_1X+\beta_2Y+\varepsilon\). \(X\) dan \(Y\) adalah koordinat stasiun hasil proyeksi Lambert Azimuthal Equal Area dalam satuan kilometer, bukan variabel lingkungan tersendiri, dan keeratan hubungan dengan ketinggian stasiun dilihat melalui koefisien korelasi pada stasiun dengan ketinggian lebih dari nol (nilai ketinggian 0 diperlakukan sebagai data hilang karena banyak stasiun pedalaman yang tercatat 0 m). Informasi tersebut digunakan untuk memahami cakupan data dan mempertimbangkan keterbatasan hasil interpolasi.

2.4 Metode Interpolasi Spasial

Seluruh metode interpolasi menggunakan semua stasiun sebagai data masukan (tanpa pembatasan radius pencarian) dan menduga nilai curah hujan pada titik-titik grid prediksi yang dijelaskan pada Subbab 2.5.

2.4.1 Inverse Distance Weighted (IDW)

Inverse Distance Weighted merupakan metode interpolasi deterministik yang memperkirakan nilai pada suatu lokasi sebagai rata-rata tertimbang nilai pengamatan di sekitarnya, dengan bobot yang menurun terhadap jarak (Caloiero et al., 2021). Stasiun yang lebih dekat memperoleh bobot lebih besar dibandingkan stasiun yang lebih jauh. Estimasi IDW dirumuskan sebagai:

\[ \hat{Z}(s_0)= \frac{\displaystyle\sum_{i=1}^{n}d_{i0}^{-p}Z(s_i)} {\displaystyle\sum_{i=1}^{n}d_{i0}^{-p}} \]

dengan:

\(\hat{Z}(s_0)\) adalah nilai curah hujan yang diestimasi pada lokasi \(s_0\)

\(Z(s_i)\) adalah nilai curah hujan pada stasiun ke-\(i\)

\(d_{i0}\) adalah jarak antara stasiun ke-\(i\) dan lokasi estimasi

\(p\) adalah parameter pangkat jarak yang mengatur besarnya pengaruh jarak terhadap bobot

\(n\) adalah jumlah stasiun pengamatan yang digunakan dalam estimasi

Pangkat jarak \(p\) sering ditetapkan sama dengan 2 (Caloiero et al., 2021; Sanusi et al., 2024), tetapi dapat pula dipilih dengan membandingkan galat prediksi pada beberapa nilai \(p\) (Tahitu et al., 2026). Pada penelitian ini, nilai \(p\) dipilih dari 1 sampai 6 berdasarkan RMSE LOOCV terkecil. IDW digunakan sebagai metode pembanding bagi metode kriging.

2.4.2 Ordinary Kriging (OK)

Ordinary Kriging merupakan metode interpolasi geostatistik yang memperhitungkan ketergantungan spasial antarobservasi melalui variogram. Metode ini mengasumsikan bahwa nilai rata-rata tidak diketahui tetapi konstan pada wilayah lokal yang dianalisis. Estimasi diperoleh sebagai kombinasi linear berbobot dari nilai pengamatan:

\[\hat{Z}(s_0)=\sum_{i=1}^{n}\lambda_i Z(s_i),\qquad \sum_{i=1}^{n}\lambda_i=1\]

dengan bobot \(\lambda_i\) ditentukan dengan meminimumkan ragam galat prediksi berdasarkan struktur variasi spasial yang dinyatakan oleh semivariogram:

\[\gamma(h)=\tfrac{1}{2}\,E\!\left[\{Z(s)-Z(s+h)\}^2\right]\]

Tahapan penerapan Ordinary Kriging meliputi

1. penyusunan semivariogram empiris pada selang jarak (lag) selebar 75 km dengan jarak maksimum 1.500 km

2. pencocokan model semivariogram teoretis dengan metode kuadrat terkecil terboboti

3. penggunaan model terpilih untuk menghasilkan estimasi dan ragam galat pada lokasi yang tidak teramati.

Tiga model kandidat dibandingkan, yaitu spherical, exponential, dan gaussian, sebagaimana lazim dilakukan pada kajian curah hujan (Bahtiyar et al., 2014; Sanusi et al., 2024). Model dengan jumlah kuadrat galat pencocokan terkecil dipilih. Ketiga model memiliki parameter nugget (\(c_0\)), partial sill (\(c\)), dan range (\(a\)). Sebagai contoh, model exponential dinyatakan sebagai:

\[\gamma(h)=c_0+c\left(1-e^{-h/a}\right),\quad h>0\]

2.4.3 Universal Kriging (UK)

Universal Kriging merupakan metode interpolasi geostatistik yang mempertimbangkan adanya tren spasial pada nilai pengamatan. Berbeda dengan Ordinary Kriging yang mengasumsikan rata-rata lokal konstan, Universal Kriging memungkinkan nilai rata-rata berubah mengikuti fungsi koordinat atau komponen tren lainnya:

\[Z(s)=\mu(s)+\varepsilon(s),\qquad \mu(s)=\beta_0+\beta_1 X+\beta_2 Y\]

dengan \(\varepsilon(s)\) adalah komponen galat acak yang memiliki struktur ketergantungan spasial. Pada penelitian ini, tren dimodelkan sebagai fungsi linear dari koordinat \(X\) dan \(Y\) dalam satuan kilometer. Ketinggian stasiun tidak digunakan sebagai komponen tren karena nilainya tidak tersedia pada seluruh titik grid prediksi dan banyak stasiun memiliki nilai ketinggian 0 yang tidak dapat dipastikan keakuratannya. Struktur variasi spasial pada residu tren dimodelkan melalui semivariogram dengan prosedur pemilihan model yang sama seperti pada Ordinary Kriging. Pemilihan model semivariogram pada Universal Kriging juga telah diterapkan pada kajian di Indonesia (Sofieyanti & Sirodj, 2024).

2.5 Sistem Koordinat dan Pemetaan

Koordinat stasiun diproyeksikan ke sistem koordinat planar Lambert Azimuthal Equal Area yang berpusat pada 2° LS dan 118° BT dengan satuan kilometer (Snyder, 1987), sehingga jarak antartitik dapat dihitung secara konsisten di seluruh wilayah Indonesia. Pengolahan data spasial dilakukan menggunakan konsep simple features (Pebesma, 2018).

Prediksi dibuat pada grid beraturan berjarak 25 km. Titik grid hanya dipertahankan apabila berada di dalam poligon batas daratan Indonesia, yaitu batas negara dari data Natural Earth skala menengah yang diakses melalui paket rnaturalearth dan disederhanakan dengan toleransi 5 km untuk mempercepat perhitungan. Hasil interpolasi disajikan sebagai peta estimasi curah hujan untuk masing-masing metode. Untuk metode kriging, disajikan pula peta simpangan baku galat prediksi sebagai ukuran ketidakpastian; IDW tidak menghasilkan ukuran tersebut.

Hasil pemetaan perlu diinterpretasikan dengan memperhatikan persebaran stasiun pengamatan. Daerah dengan stasiun yang relatif jarang dapat memiliki tingkat ketidakpastian estimasi yang berbeda dari daerah dengan pengamatan yang lebih rapat.

2.6 Evaluasi Model dengan Leave-One-Out Cross-Validation (LOOCV)

Evaluasi model dilakukan menggunakan Leave-One-Out Cross-Validation (LOOCV). Pada prosedur ini, satu pengamatan dikeluarkan secara bergantian, kemudian nilainya diprediksi berdasarkan pengamatan lainnya. Nilai prediksi selanjutnya dibandingkan dengan nilai aktual pada stasiun yang dikeluarkan tersebut. Prosedur ini serupa dengan teknik jackknife yang pernah diterapkan pada kriging curah hujan (Rohma, 2022). Seluruh stasiun menjadi titik uji sehingga evaluasi tidak bergantung pada pemilihan sejumlah kecil titik uji (Gambar 2.3).

Gambar 2.3. Prosedur evaluasi LOOCV

Pada IDW, prediksi LOOCV dihitung langsung dengan mengeluarkan stasiun yang dievaluasi dari penjumlahan bobot. Pada OK dan UK dengan seluruh stasiun sebagai tetangga, prediksi LOOCV dihitung secara eksak dengan satu kali inversi matriks kriging berukuran \(n\times n\), tanpa mengulang proses kriging sebanyak \(n\) kali (Dubrule, 1983). Jika \(A\) adalah matriks sistem kriging yang diperluas dengan basis tren (\(\mathbf{K}\) matriks kovariansi antarstasiun dan \(\mathbf{F}\) matriks basis tren) dan \(\mathbf{y}=(Z_1,\dots,Z_n,0,\dots,0)^{\top}\), maka galat LOOCV pada stasiun ke-\(i\) adalah:

\[ e_i=\frac{\left[A^{-1}\mathbf{y}\right]_i}{\left[A^{-1}\right]_{ii}},\qquad A=\begin{pmatrix}\mathbf{K}&\mathbf{F}\\ \mathbf{F}^{\top}&\mathbf{0}\end{pmatrix} \]

dan prediksi LOOCV adalah \(\hat{Z}_{-i}=Z_i-e_i\). Model semivariogram diestimasi satu kali dari seluruh data sehingga LOOCV pada metode kriging dapat sedikit lebih optimistis dibanding situasi sebenarnya, yaitu ketika model semivariogram harus diestimasi tanpa stasiun yang dievaluasi.

Kinerja model dievaluasi menggunakan ukuran kesalahan prediksi, yaitu Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), dan Mean Error (ME), sebagaimana digunakan pada validasi silang interpolasi curah hujan sebelumnya (Caloiero et al., 2021; Gräler et al., 2016). RMSE relatif, yaitu RMSE dibagi rata-rata curah hujan observasi dan dinyatakan dalam persen, digunakan untuk menggambarkan besarnya galat relatif terhadap tingkat curah hujan. Koefisien determinasi (\(R^2\)) dihitung sebagai kuadrat koefisien korelasi Pearson antara nilai aktual dan nilai prediksi LOOCV.

RMSE dirumuskan sebagai:

\[ \mathrm{RMSE}= \sqrt{\frac{1}{n}\sum_{i=1}^{n}(Z_i-\hat{Z}_i)^2} \]

MAE dirumuskan sebagai:

\[ \mathrm{MAE}= \frac{1}{n}\sum_{i=1}^{n}|Z_i-\hat{Z}_i| \]

ME dirumuskan sebagai:

\[ \mathrm{ME}= \frac{1}{n}\sum_{i=1}^{n}(\hat{Z}_i-Z_i) \]

dengan \(Z_i\) sebagai nilai aktual, \(\hat{Z}_i\) sebagai nilai prediksi LOOCV, dan \(n\) sebagai jumlah pengamatan yang dievaluasi. RMSE dan MAE yang lebih kecil menunjukkan kesalahan prediksi yang lebih rendah, sedangkan ME yang mendekati nol menunjukkan bahwa rata-rata kesalahan prediksi tidak condong ke satu arah. Nilai \(R^2\) digunakan sebagai ukuran tambahan untuk melihat kesesuaian nilai prediksi dengan nilai aktual.

Metode dengan nilai RMSE LOOCV terkecil dipilih sebagai metode dengan kinerja prediksi terbaik di antara metode yang dibandingkan. Namun, perbedaan nilai evaluasi tetap perlu dipertimbangkan secara proporsional dan tidak secara otomatis menunjukkan bahwa suatu metode unggul secara signifikan.

2.7 Perangkat Lunak

Pengolahan data, analisis statistik, pemodelan interpolasi, evaluasi prediksi, dan visualisasi dilakukan menggunakan perangkat lunak R versi 4.5.2. Paket yang digunakan mencakup dplyr untuk pengolahan data, ggplot2 untuk visualisasi, sf untuk pengolahan data spasial, gstat untuk pemodelan semivariogram dan interpolasi geostatistik (Gräler et al., 2016), rnaturalearth dan rnaturalearthdata untuk batas wilayah, serta knitr dan rmarkdown untuk penyusunan laporan.

3 Hasil dan Pembahasan

3.1 Data dan Persiapan

3.1.1 Ekstrak File ZIP

if (!file.exists(file_zip)) {
  stop("File zip tidak ditemukan. Ubah nilai 'file_zip' pada chunk 'setup'.")
}

folder_data <- tempfile()
dir.create(folder_data)
unzip(file_zip, exdir = folder_data)

3.1.2 Fungsi Pembantu

# Konversi koordinat DMS ("-06:11:00") ke derajat desimal
dms_to_decimal <- function(x) {
  x <- trimws(as.character(x))

  vapply(x, function(z) {
    if (is.na(z) || z == "") return(NA_real_)

    tanda  <- if (substr(z, 1, 1) == "-") -1 else 1
    bagian <- strsplit(sub("^[+-]", "", z), ":", fixed = TRUE)[[1]]
    if (length(bagian) != 3) return(NA_real_)

    tanda * (as.numeric(bagian[1]) +
               as.numeric(bagian[2]) / 60 +
               as.numeric(bagian[3]) / 3600)
  }, numeric(1), USE.NAMES = FALSE)
}

# Ambil STAID dari header file indexRR
cari_staid <- function(header) {
  baris_id <- grep("STAID:", header, value = TRUE)
  if (length(baris_id) == 0) return(NA_integer_)
  as.integer(sub(".*STAID:\\s*([0-9]+).*", "\\1", baris_id[1]))
}

3.1.3 Metadata Stasiun

file_stasiun <- list.files(
  folder_data, pattern = "^stations\\.txt$",
  recursive = TRUE, full.names = TRUE, ignore.case = TRUE
)
if (length(file_stasiun) == 0) stop("File stations.txt tidak ditemukan.")

stations <- read.csv(
  file_stasiun[1], skip = 15, header = TRUE,
  strip.white = TRUE, fill = TRUE
)
names(stations) <- trimws(names(stations))
stations[] <- lapply(stations, function(x) if (is.character(x)) trimws(x) else x)

stations_id <- stations %>%
  filter(toupper(COUNTRYNAME) == "INDONESIA") %>%
  mutate(STAID = as.integer(STAID))

cat("Jumlah stasiun Indonesia dalam metadata:",
    n_distinct(stations_id$STAID), "\n")

## Jumlah stasiun Indonesia dalam metadata: 8890

3.1.4 Data Curah Hujan

files_rr <- list.files(
  folder_data, pattern = "^indexRR.*\\.txt$",
  recursive = TRUE, full.names = TRUE, ignore.case = TRUE
)
cat("Jumlah file curah hujan ditemukan:", length(files_rr), "\n")

## Jumlah file curah hujan ditemukan: 10061

if (length(files_rr) == 0) stop("File indexRR tidak ditemukan.")

baca_hujan <- function(file) {
  header <- readLines(file, n = 30, warn = FALSE)
  staid  <- cari_staid(header)

  # Hanya stasiun Indonesia
  if (is.na(staid) || !staid %in% stations_id$STAID) return(NULL)

  data <- read.table(file, skip = 30, header = FALSE, na.strings = "-999999")

  if (ncol(data) != 21) {
    warning("Jumlah kolom tidak sesuai: ", basename(file))
    return(NULL)
  }

  names(data) <- c(
    "SOUID", "YEAR", "ANNUAL", "ONDJFM", "AMJJAS", "DJF", "MAM", "JJA",
    "SON", "JAN", "FEB", "MAR", "APR", "MAY", "JUN", "JUL", "AUG",
    "SEP", "OCT", "NOV", "DEC"
  )
  data$STAID <- staid
  data
}

data_indonesia <- bind_rows(lapply(files_rr, baca_hujan))
if (nrow(data_indonesia) == 0) stop("Tidak ada data curah hujan stasiun Indonesia.")

cat("Jumlah baris data curah hujan:", nrow(data_indonesia), "\n")

## Jumlah baris data curah hujan: 596917

cat("Jumlah stasiun dengan data curah hujan:",
    n_distinct(data_indonesia$STAID), "\n")

## Jumlah stasiun dengan data curah hujan: 8890

3.1.5 Sumber dan Definisi Variabel

cat("=== Header stations.txt ===\n")

## === Header stations.txt ===

cat(readLines(file_stasiun[1], n = 15, warn = FALSE), sep = "\n")

## SOUTHEAST ASIA CLIMATE ASSESSMENT & DATASET (SACA&D), file created on 02-12-2025
## THESE DATA CAN BE USED FREELY PROVIDED THAT THE FOLLOWING SOURCE IS ACKNOWLEDGED: 
## 
## van den Besselaar and Coauthors, 2015. 
## International Climate Assessment & Dataset (ICA&D): climate services across borders.
## Bulletin of the American Meteorological Society, January, 16-21,doi:10.1175/BAMS-D-13-00249.1 
## Data and metadata available at  http://sacad.database.bmkg.go.id/ 
## 
## 01-05 STAID       : Station identifier
## 07-46 STATIONNAME : Station Name
## 48-87 COUNTRYNAME : Country Name
## 89-97 LAT         : Latitude in degrees:minutes:seconds (+:North, -:South)
## 99-109 LON        : Longitude in degrees:minutes:seconds (+:East, -:West)
## 111-115 HGT       : Station elevation in meters

file_contoh <- NULL
for (f in files_rr) {
  id <- cari_staid(readLines(f, n = 30, warn = FALSE))
  if (!is.na(id) && id %in% stations_id$STAID) { file_contoh <- f; break }
}
if (!is.null(file_contoh)) {
  cat("\n=== Header salah satu file indexRR (stasiun Indonesia) ===\n")
  cat(readLines(file_contoh, n = 30, warn = FALSE), sep = "\n")
}

## 
## === Header salah satu file indexRR (stasiun Indonesia) ===
## SOUTHEAST ASIA CLIMATE ASSESSMENT & DATASET (SACA&D), file created on: 02-12-2025
## THESE DATA CAN BE USED FREELY PROVIDED THAT THE FOLLOWING SOURCE IS ACKNOWLEDGED: 
## 
## van den Besselaar and Coauthors, 2015. 
## International Climate Assessment & Dataset (ICA&D): climate services across borders. 
## Bulletin of the American Meteorological Society, January, 16-21,doi:10.1175/BAMS-D-13-00249.1 
## Data and metadata available at http://sacad.database.bmkg.go.id/ 
## 
## FILE FORMAT (MISSING VALUE CODE = -999999):
## 
## THIS FILE HOLDS DATA FOR STATION INDONESIA                                STASIUN METEOROLOGI CUT BAU MAIMUN SALEH (STAID:      1)
## LATITUDE: dd:mm:ss +05:52:36
## LONGITUDE: ddd:mm:ss +095:20:16
## ELEVATION: m  126
## 
## DATA IS FOR INDEX Precipitation sum (RR) with unit 0.01 Precipitation (mm)
## 
## 01-06 SOUID: Source identifier
## 08-11 YEAR : YYYY
## 13-20 ANNUAL DATA VALUES
## 21-28 ONDJFM DATA VALUES
## 29-36 AMJJAS DATA VALUES
## 37-44 DJF DATA VALUES
## 45-52 MAM DATA VALUES
## 53-60 JJA DATA VALUES
## 61-68 SON DATA VALUES
## 69-76 JANUARY DATA VALUES
## etc. 
## 157-164 DECEMBER DATA VALUES

Ringkasan variabel dan sumber data

Keterangan	Isi
Variabel	Curah hujan tahunan (kolom ANNUAL); satuan asli 0,01 mm, dikonversi ke mm (dibagi 100)
Tipe data	Titik (stasiun pengamatan)
Tahun	2024
Cakupan	Indonesia
Sumber data	Southeast Asia Climate Assessment & Dataset (SACA&D), BMKG, https://sacad.bmkg.go.id/. Diakses: Rabu, 14 Oktober, 2026, 22:06

3.1.6 Gabungkan Data dan Pilih Tahun 2024

data_indonesia <- data_indonesia %>% left_join(stations_id, by = "STAID")

data_indonesia_2024 <- data_indonesia %>% filter(YEAR == 2024)

cat("Jumlah baris tahun 2024:", nrow(data_indonesia_2024), "\n")

## Jumlah baris tahun 2024: 7472

cat("Jumlah stasiun tahun 2024:", n_distinct(data_indonesia_2024$STAID), "\n")

## Jumlah stasiun tahun 2024: 7472

data_spasial <- data_indonesia_2024 %>%
  select(STAID, YEAR, ANNUAL) %>%
  left_join(
    stations_id %>% select(STAID, STATIONNAME, COUNTRYNAME, LAT, LON, HGT),
    by = "STAID"
  )

cat("Jumlah baris data spasial:", nrow(data_spasial), "\n")

## Jumlah baris data spasial: 7472

3.2 Kontrol Kualitas Data

Urutan pemeriksaan:

1. data hilang dan duplikasi

2. konversi koordinat

3. validasi koordinat terhadap wilayah Indonesia

4. pemeriksaan outlier nilai curah hujan. Hanya titik yang gagal validasi koordinat yang dikeluarkan, karena posisinya salah dan akan merusak interpolasi.

Outlier nilai tidak otomatis dibuang karena curah hujan tinggi di Indonesia bisa nyata.

3.2.1 Data Hilang dan Duplikasi

cat("Jumlah NA setiap kolom:\n")

## Jumlah NA setiap kolom:

print(colSums(is.na(data_spasial)))

##       STAID        YEAR      ANNUAL STATIONNAME COUNTRYNAME         LAT 
##           0           0        2864           0           0           0 
##         LON         HGT 
##           0           0

# Satuan asli indeks RR di berkas SACA&D adalah 0,01 mm -> bagi 100 agar menjadi mm
data_analisis <- data_spasial %>%
  filter(!is.na(ANNUAL)) %>%
  mutate(ANNUAL = ANNUAL / 100)
cat("\nBaris dengan ANNUAL kosong yang dibuang:",
    nrow(data_spasial) - nrow(data_analisis), "\n")

## 
## Baris dengan ANNUAL kosong yang dibuang: 2864

duplikasi_stasiun <- data_analisis %>% count(STAID) %>% filter(n > 1)

if (nrow(duplikasi_stasiun) == 0) {
  cat("Tidak ada duplikasi STAID.\n")
} else {
  cat("Ada STAID ganda; hanya baris pertama tiap stasiun dipertahankan:\n")
  print(duplikasi_stasiun)
  data_analisis <- data_analisis %>% group_by(STAID) %>% slice(1) %>% ungroup()
}

## Tidak ada duplikasi STAID.

3.2.2 Konversi Koordinat

String koordinat asli disimpan di LAT_raw dan LON_raw agar kesalahan input di metadata bisa ditelusuri.

data_analisis <- data_analisis %>%
  mutate(
    LAT_raw = LAT, LON_raw = LON,
    LAT = dms_to_decimal(LAT),
    LON = dms_to_decimal(LON)
  )

cat("Koordinat gagal dikonversi (NA):",
    sum(is.na(data_analisis$LAT) | is.na(data_analisis$LON)), "\n")

## Koordinat gagal dikonversi (NA): 0

data_analisis <- data_analisis %>% filter(!is.na(LAT), !is.na(LON))

summary(data_analisis[, c("LAT", "LON")])

##       LAT               LON          
##  Min.   :-10.862   Min.   : -0.6203  
##  1st Qu.: -7.626   1st Qu.:106.4649  
##  Median : -6.523   Median :110.7644  
##  Mean   : -4.570   Mean   :110.9432  
##  3rd Qu.: -1.976   3rd Qu.:114.2837  
##  Max.   :  8.043   Max.   :140.8508

3.2.3 Validasi Koordinat

Titik yang jatuh di luar Indonesia dicek dengan tiga cara, dari yang paling sederhana:

1. Kotak batas (bounding box): lintang di luar -11 sampai 6 atau bujur di luar 95 sampai 141.
2. Lintang sama dengan bujur: tanda khas salah input (mis. kedua kolom terisi nilai yang sama atau tertukar).
3. Jarak ke poligon Indonesia: titik yang berjarak lebih dari 50 km dari daratan Indonesia. Titik yang gagal akan dicoba dikoreksi (tanda lintang) dan, jika tetap gagal, dikeluarkan dari analisis. Karena poligon resolusi menengah tidak memuat pulau-pulau yang sangat kecil, daftar stasiun yang dikeluarkan perlu diperiksa secara manual.

Jika koordinat yang benar sudah diverifikasi dari sumber yang bisa dipertanggungjawabkan (mis. BMKG atau halaman stasiun di situs penyedia data), isi tabel koreksi_manual di bawah. Titik tersebut akan diperbaiki, bukan dibuang.

# Isi hanya jika koordinat benar sudah diverifikasi. Contoh format:
# koreksi_manual <- tibble(STAID = 00000L, LAT = <lintang desimal>, LON = <bujur desimal>)
koreksi_manual <- tibble(STAID = integer(), LAT = numeric(), LON = numeric())

if (nrow(koreksi_manual) > 0) {
  data_analisis <- data_analisis %>%
    left_join(koreksi_manual %>% rename(LAT_k = LAT, LON_k = LON), by = "STAID") %>%
    mutate(LAT = coalesce(LAT_k, LAT), LON = coalesce(LON_k, LON)) %>%
    select(-LAT_k, -LON_k)
}

Sebagian stasiun di Jawa tercatat bertanda lintang utara (+) di metadata sumber, padahal nama stasiun dan bujurnya menunjukkan lokasi di selatan khatulistiwa (mis. stasiun Yogyakarta tercatat +07:54, seharusnya -07:54). Aturan koreksi: jika lintang melebihi batas atas, kebalikan tandanya masuk rentang Indonesia, dan bujurnya wajar, maka tanda lintang dibalik. Daftar hasil koreksi ditampilkan agar bisa diverifikasi dengan sumber. Aturan ini hanya menangkap kasus yang lintangnya melewati batas atas; salah tanda di lintang yang kecil (mendekati khatulistiwa) hanya bisa dideteksi lewat pemeriksaan jarak ke poligon di bawah.

koreksi_tanda <- data_analisis %>%
  filter(
    LAT > bbox_lat[2], -LAT >= bbox_lat[1],
    LON >= bbox_lon[1], LON <= bbox_lon[2],
    !STAID %in% koreksi_manual$STAID
  ) %>%
  transmute(STAID, STATIONNAME, LAT_raw, LON_raw,
            LAT_sebelum = LAT, LAT_sesudah = -LAT, LON)

cat("Stasiun dengan tanda lintang dibalik:", nrow(koreksi_tanda), "\n")

## Stasiun dengan tanda lintang dibalik: 11

kable(as.data.frame(koreksi_tanda), digits = 3,
      caption = "Stasiun yang tanda lintangnya dikoreksi")

Stasiun yang tanda lintangnya dikoreksi

STAID	STATIONNAME	LAT_raw	LON_raw	LAT_sebelum	LAT_sesudah	LON
9522	MANGGIS	+06:46:22	+106:48:48	6.773	-6.773	106.813
9523	CIKARANGNGLUWUNG_JPG KULON	+07:15:15	+106:37:27	7.254	-7.254	106.624
9538	BONGAS_GEGERBITUNG	+06:59:24	+107:01:41	6.990	-6.990	107.028
11839	PACIRAN	+06:52:36	+112:23:26	6.877	-6.877	112.391
14403	DISTAN KOTA YOGYAKARTA	+07:54:00	+110:27:41	7.900	-7.900	110.461
14404	KEPATIHAN YOGYAKARTA	+07:59:11	+110:22:37	7.986	-7.986	110.377
14405	POLBANGTAN YOGYAKARTA	+07:53:30	+110:36:50	7.892	-7.892	110.614
14406	STPN YOGYAKARTA	+08:02:34	+110:30:30	8.043	-8.043	110.508
14407	KEBUN BENIH TEGALREJO	+07:57:24	+110:27:06	7.957	-7.957	110.452
14451	CIUTARA	+06:47:24	+106:46:48	6.790	-6.790	106.780
14452	CILETUH	+07:15:00	+106:31:48	7.250	-7.250	106.530

data_analisis <- data_analisis %>%
  mutate(LAT = ifelse(STAID %in% koreksi_tanda$STAID, -LAT, LAT))

# Poligon batas Indonesia (opsional; perlu paket rnaturalearth + rnaturalearthdata)
indonesia_ll <- tryCatch(
  rnaturalearth::ne_countries(scale = "medium", country = "Indonesia",
                              returnclass = "sf"),
  error = function(e) NULL
)

data_analisis <- data_analisis %>%
  mutate(
    flag_bbox = LAT < bbox_lat[1] | LAT > bbox_lat[2] |
      LON < bbox_lon[1] | LON > bbox_lon[2],
    flag_lat_sama_lon = abs(LAT - LON) < 1e-8
  )

# Poligon diproyeksikan ke km dan DISEDERHANAKAN (toleransi 5 km) satu kali di sini.
# Poligon batas resolusi tinggi (mis. GADM/shp) berisi jutaan titik dan membuat
# proses spasial berjalan sangat lama; 5 km cukup untuk skala nasional.
indonesia_km <- NULL
if (!is.null(indonesia_ll)) {
  indonesia_km <- indonesia_ll %>%
    st_make_valid() %>%
    st_union() %>%
    st_transform(crs_km) %>%
    st_simplify(preserveTopology = TRUE, dTolerance = 5)
}

if (!is.null(indonesia_km)) {
  pts_cek <- st_as_sf(data_analisis, coords = c("LON", "LAT"),
                      crs = 4326, remove = FALSE) %>%
    st_transform(crs_km)
  dekat <- st_is_within_distance(
    pts_cek, indonesia_km, dist = 50, sparse = FALSE   # 50 km (satuan CRS = km)
  )[, 1]
  data_analisis$flag_poligon <- !dekat
} else {
  data_analisis$flag_poligon <- NA
  message("Poligon Indonesia tidak tersedia; pemeriksaan jarak ke poligon dilewati.")
}

kable(
  data_analisis %>%
    summarise(
      `Di luar bounding box` = sum(flag_bbox),
      `Lintang = bujur` = sum(flag_lat_sama_lon),
      `>50 km dari poligon` = sum(flag_poligon, na.rm = TRUE)
    ),
  caption = "Jumlah stasiun per jenis pemeriksaan"
)

Jumlah stasiun per jenis pemeriksaan

Di luar bounding box	Lintang = bujur	>50 km dari poligon
1	1	8

3.2.3.1 Koreksi Tanda Lintang Berdasarkan Poligon

Pemeriksaan bounding box hanya menangkap salah tanda lintang yang jatuh di luar -11 sampai 6. Salah tanda pada lintang kecil (mis. stasiun di Bengkulu tercatat +03:18 padahal seharusnya -03:18) tetap berada di dalam bounding box, tetapi terdeteksi oleh pemeriksaan jarak ke poligon karena jatuh di luar daratan Indonesia (mis. di Semenanjung Malaysia). Stasiun yang gagal pemeriksaan poligon dikoreksi jika lintang dengan tanda dibalik jatuh dalam 25 km dari daratan Indonesia. Hasil koreksi ditampilkan agar bisa diverifikasi; stasiun yang tetap gagal setelah koreksi dikeluarkan dari analisis.

koreksi_tanda2 <- data_analisis[0, ]

if (!is.null(indonesia_km)) {
  kandidat <- data_analisis %>%
    filter(flag_poligon %in% TRUE, !flag_bbox, !flag_lat_sama_lon)

  if (nrow(kandidat) > 0) {
    pts_flip <- kandidat %>%
      mutate(LAT_flip = -LAT) %>%
      st_as_sf(coords = c("LON", "LAT_flip"), crs = 4326, remove = FALSE) %>%
      st_transform(crs_km)

    flip_ok <- st_is_within_distance(
      pts_flip, indonesia_km, dist = 25, sparse = FALSE
    )[, 1]

    koreksi_tanda2 <- kandidat[flip_ok, ]
  }

}

cat("Stasiun dikoreksi (tanda lintang dibalik, hasil di dekat daratan):",
    nrow(koreksi_tanda2), "\n")

## Stasiun dikoreksi (tanda lintang dibalik, hasil di dekat daratan): 7

data_analisis <- data_analisis %>%
  mutate(
    LAT = ifelse(STAID %in% koreksi_tanda2$STAID, -LAT, LAT),
    flag_poligon = ifelse(STAID %in% koreksi_tanda2$STAID, FALSE, flag_poligon)
  )

kable(
  koreksi_tanda2 %>%
    transmute(STAID, STATIONNAME, LAT_raw, LON_raw,
              LAT_sebelum = LAT, LAT_sesudah = -LAT, LON),
  digits = 3,
  caption = "Stasiun yang tanda lintangnya dikoreksi berdasarkan pemeriksaan poligon"
)

Stasiun yang tanda lintangnya dikoreksi berdasarkan pemeriksaan poligon

STAID	STATIONNAME	LAT_raw	LON_raw	LAT_sebelum	LAT_sesudah	LON
14015	KARANG CAYO	+04:16:36	+102:53:12	4.277	-4.277	102.887
14020	KERKAP	+03:18:51	+102:09:59	3.314	-3.314	102.166
14022	NAPAL JUNGUR	+03:57:36	+102:30:18	3.960	-3.960	102.505
14023	SUKARAMI	+04:07:09	+102:31:24	4.119	-4.119	102.523
14025	RENA JAYA	+03:07:18	+101:35:43	3.122	-3.122	101.595
14253	AIR PADANG	+03:16:12	+102:30:00	3.270	-3.270	102.500
14254	PINANG RAYA	+03:11:24	+101:31:48	3.190	-3.190	101.530

Daftar stasiun yang masih bermasalah setelah koreksi. Kolom LAT_raw dan LON_raw menunjukkan teks asli di metadata, sehingga terlihat apakah masalahnya ada di file sumber atau di proses konversi.

bermasalah <- data_analisis %>%
  filter(flag_bbox | flag_lat_sama_lon | flag_poligon %in% TRUE) %>%
  select(STAID, STATIONNAME, LAT_raw, LON_raw, LAT, LON, HGT, ANNUAL,
         flag_bbox, flag_lat_sama_lon, flag_poligon)

kable(as.data.frame(bermasalah), digits = 3,
      caption = "Stasiun yang gagal salah satu pemeriksaan koordinat")

Stasiun yang gagal salah satu pemeriksaan koordinat

STAID	STATIONNAME	LAT_raw	LON_raw	LAT	LON	HGT	ANNUAL	flag_bbox	flag_lat_sama_lon	flag_poligon
12744	BERIWIT	-00:37:13	-000:37:13	-0.62	-0.62	0	4106	TRUE	TRUE	TRUE

ggplot(data_analisis, aes(LON, LAT)) +
  geom_point(color = "steelblue", size = 2) +
  geom_point(
    data = filter(data_analisis, flag_bbox | flag_lat_sama_lon),
    color = "red", size = 3
  ) +
  labs(title = "Pemeriksaan Koordinat Stasiun", x = "Bujur", y = "Lintang") +
  theme_minimal()

Titik merah: stasiun yang gagal pemeriksaan bounding box atau lintang = bujur

Metadata mentah untuk stasiun yang gagal bounding box, untuk dicocokkan dengan sumber aslinya:

id_bermasalah <- data_analisis %>%
  filter(flag_bbox | flag_lat_sama_lon) %>%
  pull(STAID)

kable(
  stations %>%
    filter(as.integer(STAID) %in% id_bermasalah) %>%
    select(STAID, STATIONNAME, COUNTRYNAME, LAT, LON, HGT),
  caption = "Metadata asli stasiun yang gagal validasi koordinat"
)

Metadata asli stasiun yang gagal validasi koordinat

STAID	STATIONNAME	COUNTRYNAME	LAT	LON	HGT
12744	BERIWIT	INDONESIA	-00:37:13	-000:37:13	0

Titik yang gagal bounding box, berciri lintang = bujur, atau masih berjarak lebih dari 50 km dari daratan Indonesia setelah koreksi tanda dikeluarkan dari analisis spasial dan dicatat sebagai temuan kualitas data. Setelah itu, stasiun dengan koordinat kembar dikurangi menjadi satu titik.

data_final <- data_analisis %>%
  filter(!(flag_bbox | flag_lat_sama_lon | flag_poligon %in% TRUE)) %>%
  mutate(HGT = suppressWarnings(as.numeric(HGT))) %>%
  select(STAID, STATIONNAME, LAT, LON, HGT, ANNUAL)

n_awal_ann   <- nrow(data_analisis)
n_gagal_koord <- n_awal_ann - nrow(data_final)

# Koordinat kembar (dalam proyeksi km) membuat matriks kriging singular:
# pertahankan satu titik untuk setiap koordinat
xy_cek <- st_coordinates(
  st_transform(st_as_sf(data_final, coords = c("LON", "LAT"), crs = 4326), crs_km)
)
dup_xy <- duplicated(round(xy_cek, 3))
n_dup  <- sum(dup_xy)
data_final <- data_final[!dup_xy, ]

cat("Stasiun sebelum validasi koordinat   :", n_awal_ann, "\n")

## Stasiun sebelum validasi koordinat   : 4608

cat("Dikeluarkan (gagal validasi koordinat):", n_gagal_koord, "\n")

## Dikeluarkan (gagal validasi koordinat): 1

cat("Dikeluarkan (koordinat kembar)        :", n_dup, "\n")

## Dikeluarkan (koordinat kembar)        : 44

cat("Stasiun dipakai untuk analisis        :", nrow(data_final), "\n")

## Stasiun dipakai untuk analisis        : 4563

Ringkasan hasil kontrol kualitas: dari 4.608 stasiun yang memiliki nilai ANNUAL, sebanyak 18 stasiun dikoreksi tanda lintangnya, 1 stasiun dikeluarkan karena koordinat tidak valid, dan 44 stasiun dikeluarkan karena berkoordinat kembar, sehingga 4.563 stasiun dipakai dalam seluruh analisis berikutnya.

3.2.4 Outlier Nilai Curah Hujan (Metode IQR)

Q1 <- quantile(data_final$ANNUAL, 0.25)
Q3 <- quantile(data_final$ANNUAL, 0.75)
IQR_value <- IQR(data_final$ANNUAL)

batas_bawah <- Q1 - 1.5 * IQR_value
batas_atas  <- Q3 + 1.5 * IQR_value

cat("Batas bawah:", round(batas_bawah, 1), "\n")

## Batas bawah: 27

cat("Batas atas:", round(batas_atas, 1), "\n")

## Batas atas: 4917.6

outlier_hujan <- data_final %>%
  filter(ANNUAL < batas_bawah | ANNUAL > batas_atas) %>%
  arrange(desc(ANNUAL))

cat("Jumlah outlier nilai:", nrow(outlier_hujan), "\n")

## Jumlah outlier nilai: 127

kable(as.data.frame(outlier_hujan), digits = 2,
      caption = "Stasiun dengan nilai ANNUAL di luar batas IQR")

Stasiun dengan nilai ANNUAL di luar batas IQR

STAID	STATIONNAME	LAT	LON	HGT	ANNUAL
9502	BATU KARUT	-6.99	106.68	0	10300.5
9504	CIKELAT	-7.38	106.46	61	9803.5
10023	REMBANG/LOSARI	-7.30	109.52	17	9685.0
8302	BEUTONG	4.26	96.42	2	9389.0
8201	SUNGAI MAS	4.61	96.20	0	7900.0
13234	GIMPU	-1.62	120.04	44	7789.0
14319	SASAK	0.01	99.70	0	7784.0
14417	TAMPINGAN	-7.10	110.29	0	7193.4
10688	LEBAKBARANG	-7.13	109.66	65	7180.0
14229	IWAKA	-4.45	136.63	2	7106.0
12567	TOHO	0.43	109.27	2	7096.0
13052	BP3K SINONSAYANG	1.12	124.34	6	6998.2
14298	PENGKADAN	0.49	112.35	4	6850.0
92	STASIUN METEOROLOGI MOZEZ KILANGIN	-4.53	136.89	30	6846.7
14505	CIAUL	-6.92	106.94	66	6782.0
10724	BELIK	-7.19	109.33	0	6737.0
10729	BONGAS	-7.19	109.49	49	6598.0
8723	PATAMUAN TANDIKAT	0.50	100.25	21	6594.4
10728	SIPEDANG	-7.14	109.49	47	6439.0
14421	KETENGER	-7.32	109.23	0	6418.0
12726	HANAU	-2.54	112.10	0	6306.0
10727	WATUKUMPUL 1	-7.16	109.42	55	6284.0
9497	GUNUNG BUNDER	-6.65	106.70	0	6170.0
14318	SUKAMENANTI	0.10	99.81	19	6127.5
8370	HAPESONG	1.47	99.07	0	6097.0
8715	SEI GERINGGING	-0.42	100.09	23	6069.5
9613	RANCABALI	-7.16	107.36	164	6065.0
14321	LUNANG	-2.00	100.98	1	5973.6
10135	TANJUNGSARI	-7.45	109.98	71	5938.0
10154	TUK MUDAL	-7.33	109.93	86	5914.0
9608	RANCASUNI	-6.94	107.74	72	5901.5
9257	BATU ROTO	-3.47	102.32	48	5876.5
8763	PECONINA	-1.59	101.14	94	5857.5
14056	BP3K TANA LILI	-2.51	120.60	0	5852.0
12736	TEWAH	-1.08	113.70	7	5848.4
10730	WATUKUMPUL BPP	-7.16	109.42	0	5836.0
14301	RUMAH KAY	-3.45	128.53	0	5799.0
14	STASIUN METEOROLOGI MINANGKABAU	-0.79	100.29	6	5776.1
9459	PAMIJAHAN	-6.62	106.64	31	5760.0
8721	KANDANG AMPEK	0.50	100.37	39	5742.1
8337	SIMPANG KIRI	2.64	97.98	0	5714.5
11199	JERUK/GUNUNG RINGGIT	-7.99	113.17	69	5683.8
8766	UPTB TALAMAU	0.23	100.97	0	5673.0
13456	SINGKALA / LAMASI TIMUR	-3.35	120.16	0	5653.0
8724	TIKU	0.37	100.91	1	5644.6
13478	STAMET MASAMBA	-2.55	120.32	5	5621.0
12033	BANJAR IRIGASI	-6.57	106.41	21	5600.5
57	STASIUN METEOROLOGI ANDI JEMMA	-2.55	120.32	43	5596.3
8340	SULTAN DAULAT	2.85	97.88	0	5595.0
14233	SILAT HULU	0.25	111.99	4	5585.0
10807	PENJALIN	-7.33	109.06	34	5578.0
9291	PINANG PELAPIS	-3.00	102.16	67	5569.5
9133	KIKIM SELATAN	-3.73	103.24	0	5561.5
13038	BP3K TAMAKO	3.44	125.56	2	5546.1
14306	LEKSULA	-3.76	126.52	0	5539.2
9625	SINGAJAYA	-7.48	107.92	59	5483.5
14144	SEUNAGAN TIMUR	4.27	96.35	0	5482.5
14256	PAGAR JATI	-3.60	102.36	8	5416.0
10079	SOMAGEDE	-7.53	109.56	11	5412.0
8722	SICAUNG	-0.55	100.29	11	5410.4
9228	ULU MANNA	-4.27	102.98	17	5388.5
9244	KEMUMU	-3.43	102.27	38	5377.0
12356	EBAN	-9.55	124.31	113	5359.0
12618	NANGA TEPUAI	0.38	112.29	4	5341.5
12616	NANGA SURUK	0.50	112.64	6	5340.0
8732	BPP BAWAN	-0.19	99.96	2	5334.5
14424	SINGASARI UPBAT / BBI	-7.39	109.18	0	5334.0
10719	SEMUGIH AFD/KT	-7.14	109.25	0	5324.0
14411	WATERPLAN SEMEN PADANG	-0.96	100.48	0	5322.6
13042	BP3K KENDAHE	3.69	125.40	11	5294.6
8717	PADANG LARING	0.33	100.09	32	5292.0
9538	BONGAS_GEGERBITUNG	-6.99	107.03	72	5291.0
10146	BEDAKAH	-7.35	109.97	122	5288.0
13314	PASIR PUTIH	-5.32	120.07	74	5285.0
10767	BUMIJAWA	-7.17	109.13	0	5277.0
13370	TOMPO LEMO - LEMO ( SDN.MENRONG )	-4.51	119.77	55	5271.0
8716	BPP PILUBANG	-0.54	100.10	1	5235.3
11991	PULOSARI	-6.33	105.92	0	5231.0
10768	BOJONG	-7.14	109.15	0	5226.0
8963	GUNUNG KERINCI	-1.89	101.58	101	5225.7
9529	CIGANGSA_SURADE	-7.33	106.55	14	5225.5
8718	ULAKAN TAPAKIS	-0.69	100.20	1	5224.0
13173	LALUNDU	-1.29	119.51	0	5219.5
8350	HUTABALANG	1.59	98.87	0	5215.7
9277	MEDAN JAYA	-2.99	101.49	2	5209.0
8933	KUANTAN MUDIK BPP	-0.63	101.42	7	5203.5
8767	KINALI	0.07	100.91	4	5202.5
9631	PAPANDAYAN	-7.37	107.71	0	5199.5
14025	RENA JAYA	-3.12	101.60	1000	5197.5
9458	PERKEB. CIANTEN	-6.72	106.60	95	5196.5
8346	MANDUAMAS	2.22	98.25	0	5185.0
8720	LIMAU PURUT	0.58	100.13	0	5176.4
10522	PANGKALAN	-6.96	110.68	0	5156.0
9285	RIMBO PENGADANG	-3.34	102.44	94	5140.0
13284	LAKATONG / CIKOANG	-5.54	119.45	1	5133.0
8162	LABUHAN HAJI BARAT	3.54	97.01	0	5123.0
9853	WANAYASA	-6.68	107.55	66	5122.0
12613	JONGKONG	0.66	112.28	3	5119.0
10674	BLADO SP	-7.07	109.83	53	5081.5
10769	KEMARON	-7.18	109.16	0	5078.0
8750	SIMPATI	-0.02	100.17	30	5052.9
14459	CIJAMPANG_PSDA	-7.13	107.15	64	5047.2
10156	WANGAN AJI	-7.30	109.92	97	5044.0
10741	SEMANGKIR AFD	-7.10	109.35	0	5040.0
14448	PERKEBUNAN CIGARU	-7.10	106.60	73	5024.5
12607	NANGA DEDAI	-0.05	111.57	3	5022.0
8758	LUBUK GADANG SELATAN	0.55	100.25	57	5001.5
13261	BORONG RAPOA	-5.39	120.04	71	4990.0
9333	TANJUNG BINGA	-2.60	107.64	1	4985.4
11988	BD.BUNUT	-6.27	106.11	0	4983.3
12030	BPP CIBEBER	-6.83	106.33	63	4979.5
8742	MUARA PAITI	0.32	100.39	13	4977.0
10153	MOJOTENGAH	-7.33	109.89	0	4969.0
8550	MUARA SOMA	0.69	99.36	0	4957.6
12608	NANGA MAU	0.06	112.49	4	4954.0
12690	BAAMANG TENGAH	-2.53	112.89	0	4952.6
9536	PERBAWATI	-6.51	106.57	112	4950.0
14142	KEB.SEUMANYAM DARUL MAKMUR	3.97	96.56	3	4949.8
12620	SEMITAU	0.50	112.00	3	4948.0
10132	KEPIL	-7.52	110.00	54	4939.0
14255	MANNA	-4.44	102.90	4	4936.0
9854	CIRACAS	-6.65	107.59	55	4935.5
8746	LUBUK SIKAPING	0.14	100.17	46	4934.4
10044	PENUSUPAN / PEJAWARAN	-7.25	109.80	0	4926.0
12622	LANJAK	0.04	112.23	4	4926.0
13297	EX SMPK BB. MALINO/BPP. TINGGI MONCONG	-5.25	119.85	102	4925.0
14018	WARA BARAT	-2.96	120.10	2	4918.0

Outlier nilai dipertahankan. Curah hujan Indonesia memang sangat bervariasi menurut topografi dan letak, sehingga nilai ekstrem belum tentu salah. Stasiun hanya dibuang jika terbukti salah ukur atau salah input.

3.2.5 Simpan Data Bersih

write.csv(data_final, "data_curah_hujan_2024_bersih.csv", row.names = FALSE)
cat("Tersimpan: data_curah_hujan_2024_bersih.csv (", nrow(data_final), "stasiun )\n")

## Tersimpan: data_curah_hujan_2024_bersih.csv ( 4563 stasiun )

3.3 Statistik Deskriptif

x <- data_final$ANNUAL
n_st <- length(x)
rata <- mean(x)
med  <- median(x)
sdv  <- sd(x)
cv   <- 100 * sdv / rata
skew <- mean((x - rata)^3) / (mean((x - rata)^2))^(3 / 2)

statistik_deskriptif <- data.frame(
  Statistik = c("Jumlah stasiun", "Rata-rata", "Median", "Simpangan baku",
                "Koefisien variasi (%)", "Minimum", "Kuartil 1", "Kuartil 3",
                "Maksimum", "Skewness"),
  Nilai = c(n_st, rata, med, sdv, cv, min(x),
            unname(quantile(x, 0.25)), unname(quantile(x, 0.75)),
            max(x), skew)
)
kable(statistik_deskriptif, digits = 2,
      caption = "Statistik deskriptif curah hujan tahunan 2024")

Statistik deskriptif curah hujan tahunan 2024

Statistik	Nilai
Jumlah stasiun	4563.00
Rata-rata	2556.18
Median	2380.50
Simpangan baku	1010.77
Koefisien variasi (%)	39.54
Minimum	413.00
Kuartil 1	1861.00
Kuartil 3	3083.65
Maksimum	10300.50
Skewness	1.29

stasiun_max <- data_final$STATIONNAME[which.max(x)]
stasiun_min <- data_final$STATIONNAME[which.min(x)]
hgt_max <- data_final$HGT[which.max(x)]
hgt_min <- data_final$HGT[which.min(x)]
q1 <- unname(quantile(x, 0.25)); q3 <- unname(quantile(x, 0.75))

if (n_st >= 3 && n_st <= 5000) {
  uji_normal <- shapiro.test(x)
  print(uji_normal)
}

## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.93068, p-value < 2.2e-16

ggplot(data_final, aes(ANNUAL)) +
  geom_histogram(bins = 20, fill = "lightblue", color = "white") +
  geom_vline(xintercept = rata, color = "red", linetype = "dashed") +
  geom_vline(xintercept = med, color = "darkgreen", linetype = "dashed") +
  labs(title = "Distribusi Curah Hujan Tahunan 2024",
       x = "Nilai ANNUAL", y = "Frekuensi") +
  theme_minimal()

Gambar 3.1. Histogram curah hujan tahunan 2024 (garis putus-putus: rata-rata dan median)

ggplot(data_final, aes(x = "", y = ANNUAL)) +
  geom_boxplot(fill = "lightgreen") +
  labs(title = "Boxplot Curah Hujan Tahunan 2024", x = NULL, y = "Nilai ANNUAL") +
  theme_minimal()

Gambar 3.2. Boxplot curah hujan tahunan 2024

Data yang dianalisis adalah curah hujan tahunan (annual precipitation) tahun 2024 dari 4.563 stasiun pengamatan di Indonesia setelah pembersihan data, koreksi dan validasi koordinat, serta penghapusan 44 titik berkoordinat kembar.

Rata-rata curah hujan tahunan sebesar 2.556,2 mm dengan simpangan baku 1.010,8 mm, atau koefisien variasi 39,5%, yang menunjukkan keragaman antarlokasi yang cukup besar. Rata-rata lebih tinggi daripada median (2.380,5 mm) dan skewness bernilai 1,29, sehingga distribusi menjulur ke kanan: sebagian besar stasiun berada pada rentang curah hujan menengah, sedangkan sejumlah kecil stasiun mencatat nilai yang jauh lebih tinggi. Uji Shapiro-Wilk menolak normalitas (W = 0,931; p < 2,2 × 10^-16). Pada ukuran sampel sebesar ini uji tersebut sangat sensitif sehingga hasilnya dibaca bersama histogram dan boxplot.

Berdasarkan batas pagar IQR (27,0 sampai 4.917,6 mm) terdapat 127 stasiun bernilai ekstrem. Stasiun tersebut tidak dibuang karena curah hujan tinggi di Indonesia dapat bersifat nyata, terutama di daerah pegunungan. Nilai tertinggi, 10.300,5 mm, tercatat di stasiun BATU KARUT (ketinggian tercatat 0 m), sedangkan nilai terendah, 413,0 mm, tercatat di stasiun SUMBERWARU/KALOR KORAN (ketinggian tercatat 3 m). Nilai-nilai ekstrem ini perlu diverifikasi terhadap sumber atau stasiun di sekitarnya karena dapat memengaruhi simpangan baku dan skewness; stasiun hanya dikeluarkan jika terbukti salah ukur atau salah input.

Keragaman yang besar antarstasiun menjadi alasan dilakukannya analisis spasial: statistik deskriptif belum menjelaskan bagaimana curah hujan tersebar secara geografis, sehingga analisis dilanjutkan dengan pemetaan dan interpolasi.

3.4 Visualisasi Spasial

pts_ll <- st_as_sf(data_final, coords = c("LON", "LAT"), crs = 4326, remove = FALSE)
pts_km <- st_transform(pts_ll, crs_km)

xy <- st_coordinates(pts_km)
pts_km$X <- xy[, 1]
pts_km$Y <- xy[, 2]

# Koordinat kembar sudah dikeluarkan pada chunk 'keluarkan'; pastikan tidak tersisa
stopifnot(!any(duplicated(round(st_coordinates(pts_km), 3))))

# Batas wilayah: poligon Indonesia jika tersedia, jika tidak convex hull stasiun
if (!is.null(indonesia_km)) {
  poly_km <- indonesia_km   # sudah disederhanakan pada chunk validasi-koordinat
} else {
  # Cadangan: gabungan lingkaran 75 km di sekitar stasiun. Lebih baik daripada
  # convex hull yang menutupi seluruh laut di antara pulau. Disarankan memasang
  # rnaturalearth + rnaturalearthdata (atau memakai shp) agar batasnya akurat.
  message("Poligon Indonesia tidak tersedia; memakai buffer 75 km di sekitar stasiun.")
  poly_km <- pts_km %>% st_geometry() %>% st_buffer(75) %>% st_union()
}

n_stasiun <- nrow(pts_km)
cat("Jumlah stasiun untuk visualisasi dan interpolasi:", n_stasiun, "\n")

## Jumlah stasiun untuk visualisasi dan interpolasi: 4563

ggplot() +
  geom_sf(data = st_sf(geometry = poly_km), fill = "grey95", color = "grey60",
          linewidth = 0.2) +
  geom_sf(data = pts_km, aes(color = ANNUAL), size = 2, alpha = 0.9) +
  scale_color_viridis_c(name = "ANNUAL") +
  coord_sf(datum = st_crs(4326)) +
  labs(title = "Sebaran Curah Hujan Tahunan Indonesia 2024") +
  theme_minimal()

Gambar 3.3. Peta sebaran curah hujan tahunan stasiun Indonesia 2024

pts_km$kelas <- cut(
  pts_km$ANNUAL,
  breaks = unique(quantile(pts_km$ANNUAL, probs = seq(0, 1, 0.2))),
  include.lowest = TRUE, dig.lab = 5
)

ggplot() +
  geom_sf(data = st_sf(geometry = poly_km), fill = "grey95", color = "grey60",
          linewidth = 0.2) +
  geom_sf(data = pts_km, aes(color = kelas), size = 2) +
  scale_color_viridis_d(name = "Kelas ANNUAL\n(kuantil)") +
  coord_sf(datum = st_crs(4326)) +
  labs(title = "Kelas Curah Hujan Tahunan Stasiun") +
  theme_minimal()

Gambar 3.4. Peta sebaran stasiun menurut kelas kuantil curah hujan

model_tren <- lm(ANNUAL ~ X + Y, data = st_drop_geometry(pts_km))
kable(summary(model_tren)$coefficients, digits = 4,
      caption = "Regresi tren ANNUAL terhadap koordinat (X, Y dalam km)")

Regresi tren ANNUAL terhadap koordinat (X, Y dalam km)

	Estimate	Std. Error	t value	Pr(>|t|)
(Intercept)	2661.0190	24.7957	107.3177	0.000
X	-0.0388	0.0183	-2.1212	0.034
Y	0.4595	0.0363	12.6519	0.000

r2_tren <- summary(model_tren)$r.squared
cat("R-kuadrat model tren:", round(r2_tren, 3), "\n")

## R-kuadrat model tren: 0.041

cor_hgt <- suppressWarnings(
  cor(pts_km$ANNUAL, ifelse(pts_km$HGT > 0, pts_km$HGT, NA_real_),
      use = "complete.obs")
)
if (length(cor_hgt) == 0 || is.na(cor_hgt)) cor_hgt <- NA_real_
prop_hgt0 <- mean(pts_km$HGT == 0, na.rm = TRUE)
cat("Proporsi stasiun dengan HGT = 0 (dianggap data hilang):", round(prop_hgt0, 3), "\n")

## Proporsi stasiun dengan HGT = 0 (dianggap data hilang): 0.285

t_tren <- summary(model_tren)$coefficients[, "t value"]
cat("Korelasi ANNUAL dengan ketinggian (hanya HGT > 0):", round(cor_hgt, 3), "\n")

## Korelasi ANNUAL dengan ketinggian (hanya HGT > 0): 0.17

Gambar di atas menampilkan lokasi stasiun yang diwarnai berdasarkan curah hujan tahunan. Stasiun sangat rapat di Pulau Jawa dan cukup rapat di Sumatra, sedangkan di Kalimantan bagian utara, Sulawesi bagian timur, Maluku, dan Papua stasiun jauh lebih jarang. Pola persebaran yang tidak merata ini penting karena ketelitian interpolasi pada suatu lokasi bergantung pada ketersediaan pengamatan di sekitarnya: estimasi di Jawa disusun dari banyak informasi lokal, sedangkan estimasi di Papua bertumpu pada sedikit stasiun yang berjauhan.

Secara visual, nilai curah hujan tinggi (hijau–kuning) tampak di sepanjang sisi barat Sumatra serta di beberapa titik di Kalimantan, Sulawesi, dan Papua bagian selatan-tengah, sedangkan nilai rendah (ungu) lebih banyak terlihat di sebagian Jawa bagian timur dan Nusa Tenggara. Pola yang beragam secara regional ini sejalan dengan temuan bahwa Indonesia terbagi atas beberapa wilayah curah hujan dengan pola musiman yang berbeda, yang dipengaruhi monsun dan suhu permukaan laut (Aldrian & Susanto, 2003). Penelitian ini tidak menguji penyebab klimatologis pola tersebut sehingga interpretasi bersifat deskriptif.

Regresi tren memakai ANNUAL (mm) sebagai variabel respons dan koordinat stasiun X (timur-barat) dan Y (utara-selatan), hasil proyeksi dalam kilometer, sebagai variabel penjelas. Model ini menghasilkan R² = 0,041, artinya tren berskala luas yang searah lemah dan hanya menjelaskan sekitar 4% keragaman data. Statistik t koefisien X dan Y masing-masing -2,12 dan 12,65; koefisien dapat signifikan secara statistik pada ribuan titik, tetapi kontribusinya terhadap keragaman tetap kecil. Korelasi curah hujan dengan ketinggian stasiun juga rendah (r = 0,17), dihitung hanya pada stasiun berketinggian lebih dari nol karena 28,5% stasiun tercatat berketinggian 0 m dan diperlakukan sebagai data hilang. Temuan ini penting untuk membaca hasil Universal Kriging pada bagian evaluasi.

3.5 Interpolasi Spasial

Seluruh perhitungan memakai proyeksi Lambert Azimuthal Equal Area dalam satuan kilometer agar jarak antartitik valid di seluruh Indonesia. Prediksi dibuat pada grid 25 km yang dipotong mengikuti batas wilayah.

3.5.1 Grid Prediksi

grid_pts <- st_make_grid(poly_km, cellsize = ukuran_sel_km, what = "centers")
grid_km  <- st_sf(geometry = grid_pts)
grid_km  <- grid_km[lengths(st_intersects(grid_km, poly_km)) > 0, ]

g_xy <- st_coordinates(grid_km)
grid_km$X <- g_xy[, 1]
grid_km$Y <- g_xy[, 2]

cat("Jumlah sel grid prediksi:", nrow(grid_km), "\n")

## Jumlah sel grid prediksi: 2963

3.5.2 Fungsi Pembantu Interpolasi

hitung_metrik <- function(obs, pred) {
  c(ME = mean(pred - obs),
    MAE = mean(abs(pred - obs)),
    RMSE = sqrt(mean((pred - obs)^2)),
    R2 = cor(obs, pred)^2)
}

# Pilih model variogram (Sph/Exp/Gau) dengan SSErr terkecil
pilih_model <- function(v) {
  kandidat <- c("Sph", "Exp", "Gau")

  hasil <- lapply(kandidat, function(m) {
    awal <- vgm(psill = max(v$gamma) * 0.8, model = m,
                range = max(v$dist) / 3, nugget = min(v$gamma))
    suppressWarnings(fit.variogram(v, awal))
  })

  sse <- sapply(hasil, function(f) attr(f, "SSErr"))
  singular <- sapply(hasil, function(f) isTRUE(attr(f, "singular")))
  sse[singular] <- Inf

  list(model = hasil[[which.min(sse)]],
       tabel = data.frame(Model = kandidat, SSErr = sse))
}

ambil_df <- function(sfobj, nama) {
  xy <- st_coordinates(sfobj)
  data.frame(X = xy[, 1], Y = xy[, 2],
             pred = sfobj$var1.pred, var = sfobj$var1.var, metode = nama)
}

3.5.3 Inverse Distance Weighted (IDW)

IDW menduga nilai di suatu lokasi sebagai rata-rata tertimbang nilai stasiun, dengan bobot \(1/d^p\). Pangkat \(p\) dipilih dari nilai yang memberi RMSE LOOCV terkecil.

D <- as.matrix(dist(st_coordinates(pts_km)))
z <- pts_km$ANNUAL

# Prediksi LOOCV untuk IDW: tiap stasiun diduga dari stasiun lain
idw_loocv_pred <- function(p) {
  W <- 1 / D^p
  diag(W) <- 0
  as.numeric(W %*% z) / rowSums(W)
}

tabel_idw <- do.call(rbind, lapply(1:6, function(p) {
  data.frame(idp = p, t(hitung_metrik(z, idw_loocv_pred(p))))
}))
idp_opt <- tabel_idw$idp[which.min(tabel_idw$RMSE)]

kable(tabel_idw, digits = 2, caption = "LOOCV IDW untuk berbagai pangkat jarak (p)")

LOOCV IDW untuk berbagai pangkat jarak (p)

idp	ME	MAE	RMSE	R2
1	-48.81	639.87	879.06	0.37
2	-28.93	455.52	678.14	0.56
3	-15.35	441.87	669.42	0.57
4	-9.89	456.75	693.47	0.55
5	-7.30	469.55	712.85	0.53
6	-6.10	478.43	726.97	0.52

cat("Pangkat optimal (RMSE terkecil): p =", idp_opt, "\n")

## Pangkat optimal (RMSE terkecil): p = 3

idw_grid <- idw(ANNUAL ~ 1, locations = pts_km, newdata = grid_km,
                idp = idp_opt, debug.level = 0)
pred_idw <- idw_loocv_pred(idp_opt)

Kinerja IDW pada berbagai pangkat jarak disajikan pada tabel di atas. RMSE terkecil diperoleh pada p = 3 (669,4 mm), sehingga nilai tersebut dipakai. Pada p = 1, RMSE mencapai 879,1 mm, dan setelah titik optimum RMSE kembali naik hingga 727,0 mm pada p = 6. Pada p rendah, estimasi terlalu dipengaruhi stasiun jauh sehingga terlalu halus; pada p tinggi, estimasi terlalu dikuasai stasiun terdekat sehingga permukaan menjadi bergerigi dan meniru nilai lokal. Karena bobot IDW hanya bergantung pada jarak, IDW cenderung membentuk gugusan nilai di sekitar stasiun yang bernilai ekstrem dan tidak menyediakan ukuran ketidakpastian.

3.5.4 Ordinary Kriging (OK)

v_emp <- variogram(ANNUAL ~ 1, pts_km, cutoff = cutoff_km, width = lebar_km)
hasil_ok <- pilih_model(v_emp)
vgm_ok <- hasil_ok$model
nama_ok <- as.character(vgm_ok$model[2])

kable(transform(hasil_ok$tabel, SSErr = formatC(SSErr, format = "e", digits = 3)),
      caption = "Perbandingan model variogram (OK); dipilih SSErr terkecil")

Perbandingan model variogram (OK); dipilih SSErr terkecil

Model	SSErr
Sph	1.296e+11
Exp	7.393e+10
Gau	2.011e+11

kable(as.data.frame(vgm_ok) %>% select(model, psill, range), digits = 2,
      caption = paste("Model variogram terpilih untuk OK:", nama_ok))

Model variogram terpilih untuk OK: Exp

model	psill	range
Nug	476921.7	0.00
Exp	574449.5	196.86

ggplot() +
  geom_point(data = v_emp, aes(dist, gamma), size = 2) +
  geom_line(data = variogramLine(vgm_ok, maxdist = cutoff_km), aes(dist, gamma),
            color = "red") +
  labs(x = "Jarak (km)", y = "Semivariansi") +
  theme_minimal()

Gambar 3.5. Variogram empiris dan model terpilih (Ordinary Kriging)

ok_grid <- krige(ANNUAL ~ 1, pts_km, newdata = grid_km, model = vgm_ok,
                 debug.level = 0)

Dari tiga model kandidat (spherical, exponential, gaussian), model Exp memiliki SSErr terkecil sehingga dipilih. Nugget sebesar 476.922 mm² atau 45,4% dari total sill (1.051.371 mm²) menunjukkan bahwa hampir separuh variasi tidak terstruktur pada jarak antarstasiun yang diamati. Hal ini dapat berasal dari variasi lokal curah hujan pada jarak dekat, galat pengukuran, atau kesalahan data. Parameter range sebesar 196,9 km pada model eksponensial berarti semivariansi mencapai sekitar 95% sill pada jarak kira-kira tiga kali nilai tersebut (± 591 km), sehingga ketergantungan spasial masih ada pada skala regional tetapi menurun cepat terhadap jarak. Karena nugget yang besar membatasi ketepatan prediksi di lokasi yang menyimpang dari tetangganya, kriging cenderung menghaluskan variasi lokal.

3.5.5 Universal Kriging (UK)

UK memakai tren orde pertama terhadap koordinat (ANNUAL ~ X + Y) sebagai komponen deterministik, lalu memodelkan galat (residu) tren dengan variogram. Ketinggian tidak dipakai sebagai kovariat karena nilainya tidak tersedia di seluruh grid prediksi.

v_res <- variogram(ANNUAL ~ X + Y, pts_km, cutoff = cutoff_km, width = lebar_km)
hasil_uk <- pilih_model(v_res)
vgm_uk <- hasil_uk$model
nama_uk <- as.character(vgm_uk$model[2])

kable(transform(hasil_uk$tabel, SSErr = formatC(SSErr, format = "e", digits = 3)),
      caption = "Perbandingan model variogram (UK); dipilih SSErr terkecil")

Perbandingan model variogram (UK); dipilih SSErr terkecil

Model	SSErr
Sph	9.700e+10
Exp	6.020e+10
Gau	1.321e+11

kable(as.data.frame(vgm_uk) %>% select(model, psill, range), digits = 2,
      caption = paste("Model variogram terpilih untuk UK:", nama_uk))

Model variogram terpilih untuk UK: Exp

model	psill	range
Nug	447389.5	0.00
Exp	524890.2	136.86

uk_grid <- krige(ANNUAL ~ X + Y, pts_km, newdata = grid_km, model = vgm_uk,
                 debug.level = 0)
n_negatif <- sum(uk_grid$var1.pred < 0, na.rm = TRUE)
cat("Sel grid UK dengan prediksi negatif:", n_negatif, "\n")

## Sel grid UK dengan prediksi negatif: 0

ggplot() +
  geom_point(data = v_res, aes(dist, gamma), size = 2) +
  geom_line(data = variogramLine(vgm_uk, maxdist = cutoff_km), aes(dist, gamma),
            color = "red") +
  labs(x = "Jarak (km)", y = "Semivariansi") +
  theme_minimal()

Gambar 3.6. Variogram residu dan model terpilih (Universal Kriging)

UK memodelkan tren orde pertama terhadap koordinat dan menyusun variogram pada residunya; model terpilih adalah Exp dengan nugget 46,0% dari total sill dan range 136,9 km (pada OK: 45,4% dan 196,9 km). Karena tren koordinat hanya menjelaskan 4,1% keragaman, perubahan struktur residu dibandingkan OK relatif kecil. Jumlah sel grid UK dengan prediksi negatif: 0.

3.5.6 Peta Hasil Interpolasi

nama_idw <- paste0("IDW (p = ", idp_opt, ")")
nama_okl <- paste0("Ordinary Kriging (", nama_ok, ")")
nama_ukl <- paste0("Universal Kriging (", nama_uk, ")")

df_pred <- bind_rows(
  ambil_df(idw_grid, nama_idw),
  ambil_df(ok_grid, nama_okl),
  ambil_df(uk_grid, nama_ukl)
) %>% mutate(metode = factor(metode, levels = c(nama_idw, nama_okl, nama_ukl)))

ggplot() +
  geom_tile(data = df_pred, aes(X, Y, fill = pred),
            width = ukuran_sel_km, height = ukuran_sel_km) +
  geom_sf(data = st_sf(geometry = poly_km), fill = NA, color = "grey30",
          linewidth = 0.2) +
  geom_sf(data = pts_km, size = 0.12, color = "black", alpha = 0.3) +
  scale_fill_viridis_c(name = "ANNUAL (mm)") +
  labs(x = "Bujur", y = "Lintang") +
  coord_sf(datum = st_crs(4326)) +
  facet_wrap(~metode, ncol = 1) +
  theme_minimal()

Gambar 3.7. Peta hasil interpolasi curah hujan tahunan 2024 (titik: stasiun)

Ketiga peta menampilkan pola regional yang sangat mirip: curah hujan relatif tinggi di sisi barat Sumatra dan beberapa wilayah Kalimantan serta Papua bagian selatan-tengah, dan lebih rendah di Jawa bagian timur serta Nusa Tenggara. Perbedaan utama ada pada kehalusan permukaan. IDW memunculkan titik-titik setempat bernilai tinggi di sekitar stasiun tertentu, sedangkan OK dan UK menghasilkan permukaan yang lebih halus karena memakai struktur korelasi spasial. Di daerah dengan stasiun yang sangat jarang, seperti bagian tengah Papua dan Kalimantan bagian utara, estimasi pada ketiga metode lebih banyak ditentukan oleh beberapa stasiun yang jauh sehingga perlu ditafsirkan dengan hati-hati. Perbedaan visual antarmetode tidak dengan sendirinya menunjukkan akurasi; pemilihan metode terbaik didasarkan pada LOOCV.

3.5.7 Ketidakpastian Hasil Interpolasi

Peta ketidakpastian kriging (simpangan baku galat prediksi). IDW tidak menghasilkan ukuran ketidakpastian.

df_sd <- df_pred %>%
  filter(metode != nama_idw) %>%
  mutate(sd = sqrt(pmax(var, 0)))

ggplot() +
  geom_tile(data = df_sd, aes(X, Y, fill = sd),
            width = ukuran_sel_km, height = ukuran_sel_km) +
  geom_sf(data = st_sf(geometry = poly_km), fill = NA, color = "grey30",
          linewidth = 0.2) +
  geom_sf(data = pts_km, size = 0.12, color = "black", alpha = 0.3) +
  scale_fill_viridis_c(name = "SD galat (mm)", option = "magma", direction = -1) +
  labs(x = "Bujur", y = "Lintang") +
  coord_sf(datum = st_crs(4326)) +
  facet_wrap(~metode, ncol = 1) +
  theme_minimal()

Gambar 3.8. Simpangan baku galat prediksi kriging

Simpangan baku galat prediksi kriging berkisar antara 692 dan 1.026 mm. Nilai terendah terlihat di wilayah padat stasiun seperti Jawa dan sebagian besar Sumatra, sedangkan nilai tertinggi di wilayah yang stasiunnya paling jarang, yaitu Papua bagian selatan-tengah dan Kalimantan bagian utara. Pola OK dan UK hampir identik. Dua hal perlu dicermati. Pertama, SD galat di wilayah padat pun masih relatif besar karena nugget yang besar, sehingga menambah stasiun di wilayah yang sudah rapat hanya sedikit menurunkan ketidakpastian. Kedua, selisih SD galat antara wilayah padat dan jarang tidak terlalu besar, sehingga peta ini lebih menegaskan daerah mana yang paling sedikit didukung data daripada memberikan batas kesalahan yang presisi. Simpangan baku galat dihitung berdasarkan model variogram yang dianggap benar sehingga tidak menjamin bahwa kesalahan aktual selalu berada dalam rentang tersebut; IDW tidak menghasilkan ukuran ketidakpastian yang sebanding.

3.6 Evaluasi dan Perbandingan Akurasi Metode (LOOCV)

Setiap stasiun secara bergantian dikeluarkan, diduga dari stasiun lain, lalu dibandingkan dengan nilai sebenarnya. Variogram dipilih memakai seluruh data sehingga hasil kriging sedikit lebih optimistis dibanding situasi sebenarnya. Catatan ini perlu disebutkan saat membahas hasil.

LOOCV untuk OK dan UK dihitung dengan rumus eksak dari satu kali inversi matriks kovariansi (galat LOOCV stasiun ke-\(i\) sama dengan elemen ke-\(i\) dari \(A^{-1}y\) dibagi elemen diagonal ke-\(i\) dari \(A^{-1}\), dengan \(A\) sistem kriging yang diperluas dengan basis tren; Dubrule, 1983), bukan dengan menjalankan kriging ulang untuk setiap stasiun. Hasilnya sama dengan krige.cv() berlingkungan global, tetapi selesai dalam hitungan menit untuk ribuan stasiun, sedangkan krige.cv() dengan 4.000 lebih stasiun sangat lambat.

# Kovariansi dari model variogram gstat (Nug, Exp, Sph, Gau): C(h) = sill - gamma(h)
kov_dari_vgm <- function(D, vgm_model) {
  K <- matrix(0, nrow(D), ncol(D))
  for (i in seq_len(nrow(vgm_model))) {
    m  <- as.character(vgm_model$model[i])
    ps <- vgm_model$psill[i]
    r  <- vgm_model$range[i]
    K <- K + switch(m,
      Nug = ps * (D == 0),
      Exp = ps * exp(-D / r),
      Sph = ps * ifelse(D < r, 1 - 1.5 * D / r + 0.5 * (D / r)^3, 0),
      Gau = ps * exp(-(D / r)^2),
      stop("Model variogram belum didukung: ", m)
    )
  }
  K
}

# LOOCV eksak untuk kriging dengan basis tren F (OK: kolom 1; UK: 1, X, Y)
loocv_kriging <- function(K, F, z) {
  diag(K) <- diag(K) * (1 + 1e-8)          # jaga agar K definit positif
  Kinv <- chol2inv(chol(K))
  KF   <- Kinv %*% F
  Sinv <- solve(crossprod(F, KF))
  d    <- diag(Kinv) - rowSums((KF %*% Sinv) * KF)
  Kz   <- as.numeric(Kinv %*% z)
  Az   <- Kz - as.numeric(KF %*% (Sinv %*% crossprod(F, Kz)))
  data.frame(observed = z, var1.pred = z - Az / d)
}

K_ok  <- kov_dari_vgm(D, vgm_ok)
cv_ok <- loocv_kriging(K_ok, matrix(1, nrow(D), 1), z)
rm(K_ok)

K_uk  <- kov_dari_vgm(D, vgm_uk)
F_uk  <- cbind(1, as.numeric(scale(pts_km$X)), as.numeric(scale(pts_km$Y)))
cv_uk <- loocv_kriging(K_uk, F_uk, z)
rm(K_uk)

tabel_cv <- bind_rows(
  data.frame(Metode = nama_idw, t(hitung_metrik(z, pred_idw))),
  data.frame(Metode = nama_okl,
             t(hitung_metrik(cv_ok$observed, cv_ok$var1.pred))),
  data.frame(Metode = nama_ukl,
             t(hitung_metrik(cv_uk$observed, cv_uk$var1.pred)))
) %>%
  mutate(`RMSE relatif (%)` = 100 * RMSE / mean(z)) %>%
  arrange(RMSE)

kable(tabel_cv, digits = 3,
      caption = "Perbandingan akurasi interpolasi berdasarkan LOOCV (urut RMSE)")

Perbandingan akurasi interpolasi berdasarkan LOOCV (urut RMSE)

Metode	ME	MAE	RMSE	R2	RMSE relatif (%)
Universal Kriging (Exp)	-1.047	441.933	654.067	0.584	25.588
Ordinary Kriging (Exp)	-0.844	448.887	661.519	0.575	25.879
IDW (p = 3)	-15.349	441.870	669.416	0.566	26.188

terbaik <- tabel_cv$Metode[1]

df_cv <- bind_rows(
  data.frame(obs = z, pred = pred_idw, metode = nama_idw),
  data.frame(obs = cv_ok$observed, pred = cv_ok$var1.pred, metode = nama_okl),
  data.frame(obs = cv_uk$observed, pred = cv_uk$var1.pred, metode = nama_ukl)
) %>% mutate(metode = factor(metode, levels = c(nama_idw, nama_okl, nama_ukl)))

ggplot(df_cv, aes(obs, pred)) +
  geom_point(alpha = 0.6, color = "steelblue") +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
  facet_wrap(~metode) +
  labs(x = "Pengamatan", y = "Prediksi LOOCV") +
  theme_minimal()

Gambar 3.9. Nilai pengamatan vs prediksi LOOCV

Berdasarkan RMSE LOOCV, urutan metode dari yang terbaik adalah: Universal Kriging (Exp) (654,1 mm), Ordinary Kriging (Exp) (661,5 mm), IDW (p = 3) (669,4 mm). Metode terbaik, Universal Kriging (Exp), memiliki RMSE 654,1 mm (25,6% dari rata-rata ANNUAL), MAE 441,9 mm, dan R² 0,584.

Selisih RMSE antara metode terbaik dan terburuk hanya 2,3%, sehingga ketiga metode berkinerja hampir setara dan keunggulan metode terbaik tidak boleh dibesar-besarkan; penelitian ini juga tidak menguji signifikansi statistik selisih tersebut. UK sedikit lebih baik daripada OK (selisih RMSE 1,1%), selaras dengan temuan bahwa tren koordinat lemah (R² = 0,041): pemodelan tren hanya menambah sedikit informasi sehingga peningkatan akurasinya juga kecil. OK memiliki RMSE lebih kecil daripada IDW tetapi MAE yang sedikit lebih besar (448,9 mm vs 441,9 mm); ini menunjukkan bahwa kriging lebih baik menekan galat besar (RMSE memberi bobot lebih pada galat besar), sedangkan galat tipikal kedua metode setara.

Mean error ketiga metode berkisar -15,35 sampai -0,84 mm, atau tidak lebih dari 0,6% dari rata-rata curah hujan, sehingga tidak ada bias sistematis yang berarti (ME negatif menandakan kecenderungan menaksir sedikit lebih rendah). RMSE relatif 25,6%–26,2% dan R² 0,57–0,58 menunjukkan bahwa sekitar 42% keragaman nilai pengamatan tidak dijelaskan oleh prediksi LOOCV metode terbaik. Besarnya nugget (45% dari sill) sejalan dengan batas kemampuan prediksi ini: variasi lokal yang tidak terstruktur secara spasial tidak dapat ditangkap oleh metode interpolasi mana pun yang hanya memakai informasi dari stasiun sekitar.

3.6.1 Hubungan Nilai Pengamatan dan Prediksi

Pada diagram pencar di atas, titik yang dekat garis diagonal 1:1 menunjukkan prediksi yang mendekati nilai aktual. Sebaran titik di sekitar garis relatif serupa pada ketiga metode, sesuai dengan selisih RMSE yang kecil. Nilai pengamatan yang tinggi cenderung berada di bawah garis 1:1, artinya nilai tinggi cenderung diprediksi lebih rendah; hal ini lazim pada interpolasi karena prediksi merupakan rata-rata terboboti nilai di sekitarnya sehingga menghaluskan nilai ekstrem. Stasiun bernilai sangat tinggi, termasuk yang perlu diverifikasi (stasiun BATU KARUT), termasuk yang paling sulit diprediksi dan dapat menaikkan RMSE. Dengan demikian, hasil interpolasi lebih tepat dipakai sebagai gambaran pola spasial regional daripada sebagai pengganti pengamatan langsung pada setiap lokasi.

3.7 Keterbatasan Analisis

Pertama, persebaran stasiun tidak merata sehingga ketelitian estimasi berbeda antarwilayah, terutama di Indonesia bagian timur. Kedua, tren pada UK hanya memakai koordinat X dan Y; ketinggian tidak dipakai karena tidak tersedia di seluruh titik grid dan 28,5% stasiun bernilai 0 m, serta faktor lain seperti monsun atau topografi tidak dimasukkan. Ketiga, variogram diestimasi dari seluruh data sebelum LOOCV sehingga hasil LOOCV kriging dapat sedikit optimistis dibanding situasi sebenarnya (Dubrule, 1983). Keempat, sebagian nilai ekstrem dan koordinat pada sumber data tidak dapat diverifikasi sepenuhnya; keputusan mempertahankan atau mengeluarkannya mengikuti aturan pada Bab 2. Kelima, penelitian ini hanya memakai satu tahun pengamatan (2024) sehingga hasilnya tidak mewakili kondisi tahun lain.

Pengembangan selanjutnya dapat mencakup kovariat lingkungan (misalnya ketinggian dari data DEM), prosedur validasi yang membentuk ulang model variogram pada setiap iterasi, serta penggunaan data beberapa tahun.

4 Kesimpulan

1. Data yang dianalisis adalah curah hujan tahunan 2024 dari 4.563 stasiun di Indonesia setelah validasi koordinat dan pembersihan data. Rata-rata curah hujan tahunan adalah 2.556,2 mm dengan simpangan baku 1.010,8 mm; distribusinya menjulur ke kanan.
2. Stasiun terkonsentrasi di Jawa dan Sumatra, sedangkan wilayah timur Indonesia jarang sehingga estimasi dan ketidakpastian di wilayah tersebut perlu ditafsirkan dengan hati-hati. Tren spasial linear terhadap koordinat lemah (R² = 0,041).
3. IDW (p = 3), Ordinary Kriging (Exp), dan Universal Kriging (Exp) menghasilkan pola spasial yang mirip. Berdasarkan LOOCV, metode dengan RMSE terkecil adalah Universal Kriging (Exp) (RMSE 654,1 mm; R² 0,584), tetapi selisih terhadap metode lain hanya 2,3%.
4. Keterbatasan utama: sebaran stasiun tidak merata, UK hanya memakai tren koordinat, variogram diestimasi dari seluruh data, dan hanya satu tahun pengamatan yang dianalisis.

5 Daftar Pustaka

Aldrian, E., & Susanto, R. D. (2003). Identification of three dominant rainfall regions within Indonesia and their relationship to sea surface temperature. International Journal of Climatology, 23(12), 1435–1452. https://doi.org/10.1002/joc.950

Bahtiyar, A. D. R., Hoyyi, A., & Yasin, H. (2014). Ordinary kriging dalam estimasi curah hujan di Kota Semarang. Jurnal Gaussian, 3(2), 151–159. https://ejournal3.undip.ac.id/index.php/gaussian/article/view/5900

Caloiero, T., Pellicone, G., Modica, G., & Guagliardi, I. (2021). Comparative analysis of different spatial interpolation methods applied to monthly rainfall as support for landscape management. Applied Sciences, 11(20), 9566. https://doi.org/10.3390/app11209566

Dubrule, O. (1983). Cross validation of kriging in a unique neighborhood. Mathematical Geology, 15, 687–699. https://doi.org/10.1007/BF01033232

Gräler, B., Pebesma, E., & Heuvelink, G. (2016). Spatio-temporal interpolation using gstat. The R Journal, 8(1), 204–218. https://doi.org/10.32614/RJ-2016-014

Najwa, S., Pratiwi, D. A., Ahdian, M. R., Indriani, A., Mindra, I. G. N., Falah, A. N., & Ruchjana, B. N. (2025). Integrating spatial autoregressive exogenous with ordinary kriging for improved rainfall prediction in Java: Enhancing accuracy with climate variables and spatial autocorrelation. InPrime: Indonesian Journal of Pure and Applied Mathematics, 7(1), 25–43. https://doi.org/10.15408/inprime.v7i1.42070

Pebesma, E. (2018). Simple features for R: Standardized support for spatial vector data. The R Journal, 10(1), 439–446. https://doi.org/10.32614/RJ-2018-009

Rohma, N. N. (2022). Estimation of ordinary kriging method with jackknife technique on rainfall data in Malang Raya. International Journal on Information and Communication Technology (IJoICT), 8(2), 22–39. https://doi.org/10.21108/ijoict.v8i2.678

Sanusi, W., Sidjara, S., Patahuddin, S., & Danial, M. (2024). A comparison of spatial interpolation methods for regionalizing maximum daily rainfall data in South Sulawesi, Indonesia. ITM Web of Conferences, 58, 04003. https://doi.org/10.1051/itmconf/20245804003

Snyder, J. P. (1987). Map projections: A working manual (U.S. Geological Survey Professional Paper 1395). U.S. Geological Survey. https://doi.org/10.3133/pp1395

Sofia, D. A., Sujono, J., & Legono, D. (2018). Analisis variabilitas spasial dan temporal curah hujan di wilayah Gunung Merapi. Jurnal Teknisia, XXIII(1). https://journal.uii.ac.id/teknisia/article/download/10679/8757/25443

Sofieyanti, N., & Sirodj, D. A. N. (2024). Penerapan metode universal kriging untuk mengestimasi laju penurunan muka tanah DKI Jakarta. Jurnal Gaussian, 13(2), 431–442. https://doi.org/10.14710/j.gauss.13.2.431-442

Tahitu, J. C., Titaley, J., & Manurung, T. (2026). Penerapan metode Inverse Distance Weighted (IDW) dalam interpolasi curah hujan (Studi kasus: Kabupaten Langkat). AITI: Jurnal Teknologi Informasi, 23(1), 61–72. https://doi.org/10.24246/aiti.v23i1.61-72

van den Besselaar, E. J. M., et al. (2015). International Climate Assessment & Dataset: Climate services across borders. Bulletin of the American Meteorological Society, 96(1), 16–21. https://doi.org/10.1175/BAMS-D-13-00249.1. Sumber data: SACA&D, http://sacad.database.bmkg.go.id/

6 Lampiran: Informasi Sesi

sessionInfo()

## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: Asia/Jakarta
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] knitr_1.51    gstat_2.1-6   sf_1.1-0      ggplot2_4.0.2 dplyr_1.2.0  
## 
## loaded via a namespace (and not attached):
##  [1] s2_1.1.9                sass_0.4.10             generics_0.1.4         
##  [4] class_7.3-23            KernSmooth_2.23-26      lattice_0.22-7         
##  [7] digest_0.6.39           magrittr_2.0.4          evaluate_1.0.5         
## [10] grid_4.5.2              RColorBrewer_1.1-3      rnaturalearth_1.2.0    
## [13] fastmap_1.2.0           jsonlite_2.0.0          e1071_1.7-17           
## [16] DBI_1.3.0               viridisLite_0.4.3       scales_1.4.0           
## [19] jquerylib_0.1.4         abind_1.4-8             cli_3.6.5              
## [22] rlang_1.1.7             units_1.0-0             intervals_0.15.5       
## [25] withr_3.0.2             cachem_1.1.0            yaml_2.3.12            
## [28] FNN_1.1.4.1             otel_0.2.0              parallel_4.5.2         
## [31] tools_4.5.2             spacetime_1.3-4         vctrs_0.7.1            
## [34] R6_2.6.1                zoo_1.8-15              proxy_0.4-29           
## [37] lifecycle_1.0.5         classInt_0.4-11         pkgconfig_2.0.3        
## [40] pillar_1.11.1           bslib_0.10.0            gtable_0.3.6           
## [43] glue_1.8.0              Rcpp_1.1.1              rnaturalearthdata_1.0.0
## [46] xfun_0.61               tibble_3.3.1            tidyselect_1.2.1       
## [49] rstudioapi_0.18.0       farver_2.1.2            htmltools_0.5.9        
## [52] labeling_0.4.3          rmarkdown_2.30          xts_0.14.1             
## [55] wk_0.9.5                compiler_4.5.2          S7_0.2.1               
## [58] stars_0.7-3             sp_2.2-3