Interpolasi Spasial Oksigen Terlarut (DO) Sungai dan Danau di Provinsi Aceh Tahun 2020
Perbandingan IDW, Ordinary Kriging, dan Universal Kriging dengan Validasi Silang Leave-One-Out
Rahayu Puspita Sari (3338240049)
Pendahuluan
Latar Belakang
Oksigen terlarut (dissolved oxygen, DO) adalah salah satu parameter utama kualitas air. Kadar DO yang rendah menandakan adanya beban bahan organik atau pencemaran yang menguras oksigen, dan dapat mengganggu kehidupan biota air. Kadar DO dipengaruhi suhu, arus dan turbulensi aliran, serta aktivitas di sekitar badan air, sehingga nilainya berbeda dari satu lokasi ke lokasi lain.
Pemantauan kualitas air sungai dan danau hanya dilakukan pada sejumlah lokasi sampling, sedangkan informasi kualitas air dibutuhkan juga di antara lokasi tersebut. Interpolasi spasial menjawab kebutuhan ini dengan menduga nilai pada lokasi yang tidak diamati dari nilai lokasi sekitarnya, berdasarkan Hukum Pertama Geografi Tobler: lokasi yang berdekatan cenderung bernilai mirip.
Penelitian ini menerapkan interpolasi spasial pada data DO sungai dan danau di Provinsi Aceh tahun 2020. Tiga metode dibandingkan, yaitu Inverse Distance Weighted (IDW), Ordinary Kriging (OK), dan Universal Kriging (UK), dan akurasinya dievaluasi dengan validasi silang Leave-One-Out (LOOCV).
Tujuan
- Menyajikan data numerik lingkungan bertipe titik (DO) di Provinsi Aceh beserta sumbernya.
- Meringkas karakteristik data melalui statistik deskriptif.
- Memvisualisasikan sebaran spasial DO dalam bentuk peta.
- Menerapkan interpolasi IDW, OK, dan UK, kemudian menentukan model terbaik berdasarkan LOOCV.
Metode Penelitian
Data dan Sumber
| Keterangan | Isi |
|---|---|
| Variabel | Oksigen terlarut / DO (mg/L) |
| Tipe data | Titik (lokasi sampling sungai dan danau, memiliki koordinat) |
| Wilayah | Provinsi Aceh |
| Tahun | 2020 |
| Sumber | Satu Data Indonesia (data.go.id), dataset Parameter Kualitas Air Sungai dan Danau, Pemerintah Provinsi Aceh, dihasilkan oleh DLHK |
| Waktu publikasi | Dipublikasikan 7 November 2024, dimodifikasi 31 Desember 2024 |
| Tautan | https://data.go.id/dataset/dataset/parameter-kualitas-air-sungai-dan-danau |
| Batas wilayah | Shapefile batas Provinsi Aceh |
Halaman sumber menyebut bahwa data masih dalam proses pemenuhan
Prinsip Satu Data Indonesia (SDI), sehingga kelengkapan metadata perlu
dipertimbangkan sebagai keterbatasan. Pada berkas, kolom koordinat
bernama Lintang_Selatan, tetapi nilainya adalah lintang
utara (Aceh berada di utara ekuator).
Pra-pemrosesan Data
Tahapan pra-pemrosesan:
- Koordinat berformat derajat-menit-detik (DMS) dikonversi menjadi
derajat desimal. Format yang tidak lazim, seperti
0950, diperbaiki menjadi095. - Nilai DO yang bertanda
<dibaca sebagai angka batas yang tertera (penyederhanaan). - Data dibatasi tahun 2020, lalu dirata-ratakan per lokasi unik (koordinat sama) agar setiap lokasi hanya memiliki satu nilai dan tidak ada koordinat kembar yang mengganggu perhitungan kriging.
- Hanya titik yang berada di dalam batas Provinsi Aceh yang dipakai.
# Sesuaikan dengan lokasi berkas di komputer Anda
file_csv <- "C:/Users/rahay/Downloads/data-kualitas-air-sungai.csv"
file_shp <- "C:/Users/rahay/Downloads/Provinsi Aceh-PROV/Provinsi Aceh-PROV/Provinsi Aceh-PROV_PL.shp"
data_air <- read.csv(
file_csv,
sep = ";",
fileEncoding = "latin1",
stringsAsFactors = FALSE,
check.names = FALSE,
na.strings = c("", "NA", "NULL")
)
batas_aceh <- st_read(file_shp, quiet = TRUE)
batas_aceh <- st_transform(batas_aceh, 4326)
cat("Dimensi data:", dim(data_air), "\n")## Dimensi data: 150 34dms_to_decimal <- function(x) {
vapply(as.character(x), function(z) {
if (is.na(z) || trimws(z) == "") return(NA_real_)
z <- gsub(",", ".", trimws(z), fixed = TRUE)
angka <- regmatches(z, gregexpr("[0-9]+(?:\\.[0-9]+)?", z, perl = TRUE))[[1]]
if (length(angka) < 2) return(NA_real_)
derajat <- as.numeric(angka[1])
# Perbaiki format seperti 0950 menjadi 095
if (derajat > 180 && derajat %% 10 == 0) derajat <- derajat / 10
menit <- as.numeric(angka[2])
detik <- if (length(angka) >= 3) as.numeric(angka[3]) else 0
if (menit >= 60 || detik >= 60) return(NA_real_)
derajat + menit / 60 + detik / 3600
}, numeric(1))
}
data_air$lat <- dms_to_decimal(data_air$Lintang_Selatan)
data_air$lon <- dms_to_decimal(data_air$Bujur_Timur)
# Konversi DO menjadi numerik
data_air$DO_num <- as.character(data_air$`Do_(mg/L)`)
data_air$DO_num <- gsub("<", "", data_air$DO_num, fixed = TRUE)
data_air$DO_num <- gsub(",", ".", data_air$DO_num, fixed = TRUE)
data_air$DO_num <- suppressWarnings(as.numeric(data_air$DO_num))
summary(data_air$lat)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.597 4.394 4.758 4.772 5.108 5.881summary(data_air$lon)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 95.32 95.72 96.41 96.43 96.98 97.90data_2020 <- data_air %>% filter(tahun == 2020)
tidak_valid <- !is.finite(data_2020$lat) | !is.finite(data_2020$lon) |
!is.finite(data_2020$DO_num)
data_analisis <- data_2020 %>%
filter(is.finite(lat), is.finite(lon), is.finite(DO_num),
lat >= 2, lat <= 7, lon >= 94, lon <= 99)
# Satu nilai per lokasi unik: rata-rata DO dari seluruh pengukuran di lokasi tersebut
data_titik <- data_analisis %>%
group_by(lon, lat) %>%
summarise(DO = mean(DO_num, na.rm = TRUE),
jumlah_pengukuran = n(), .groups = "drop") %>%
filter(is.finite(DO))
titik_do <- st_as_sf(data_titik, coords = c("lon", "lat"), crs = 4326, remove = FALSE)
# Hanya titik yang berada di dalam batas Aceh
wilayah_aceh <- st_union(batas_aceh)
di_dalam <- st_intersects(titik_do, wilayah_aceh, sparse = FALSE)[, 1]
titik_aceh <- titik_do[di_dalam, ]
# Titik di luar batas (untuk diperiksa)
st_drop_geometry(titik_do[!di_dalam, ])## # A tibble: 1 × 4
## lon lat DO jumlah_pengukuran
## * <dbl> <dbl> <dbl> <int>
## 1 97.9 5.85 5.20 2nilai_DO <- titik_aceh$DO
n <- length(nilai_DO)Berdasarkan hasil pengolahan data, terdapat 149 pengukuran kualitas air sungai di Provinsi Aceh pada tahun 2020. Seluruh pengukuran memiliki koordinat dan nilai Dissolved Oxygen (DO) yang valid berdasarkan pemeriksaan awal. Setelah pengukuran dengan koordinat yang sama dikelompokkan dan dirata-ratakan, diperoleh 79 lokasi unik. Dari jumlah tersebut, 78 lokasi berada di dalam batas administrasi Provinsi Aceh dan digunakan dalam analisis spasial, sedangkan satu lokasi berada di luar batas wilayah, yaitu pada koordinat 97,90° BT dan 5,85° LU, sehingga dikeluarkan dari analisis. Dengan demikian, analisis interpolasi DO dilakukan berdasarkan 78 lokasi pengamatan yang berada di dalam wilayah Aceh.
Statistik Deskriptif
Ringkasan data meliputi ukuran pemusatan (mean, median), penyebaran (simpangan baku, rentang, kuartil, koefisien variasi), bentuk distribusi (skewness, kurtosis), uji normalitas Shapiro-Wilk, serta identifikasi pencilan dengan aturan 1,5 × IQR. Hasilnya dilengkapi histogram dan boxplot.
Interpolasi IDW
IDW menduga nilai di lokasi \(s_0\) sebagai rata-rata terbobot dari nilai titik pengamatan, dengan bobot berbanding terbalik terhadap jarak:
\[\hat{Z}(s_0)=\frac{\sum_{i=1}^{n} w_i\,Z(s_i)}{\sum_{i=1}^{n} w_i},\qquad w_i=\frac{1}{d(s_0,s_i)^{p}}\]
dengan \(p\) adalah parameter power. Nilai \(p\) yang kecil menghasilkan permukaan halus, sedangkan \(p\) yang besar membuat dugaan didominasi titik terdekat. Nilai \(p\) dipilih dari kandidat \(\{0{,}25; 0{,}5; 0{,}75; 1; 1{,}5; 2; 2{,}5; 3; 4\}\) berdasarkan RMSE LOOCV terkecil.
Ordinary Kriging (OK)
Kriging menduga nilai sebagai kombinasi linear nilai terukur dengan bobot \(\lambda_i\) yang diturunkan dari struktur autokorelasi spasial:
\[\hat{Z}(s_0)=\sum_{i=1}^{n}\lambda_i Z(s_i),\qquad \sum_{i=1}^{n}\lambda_i=1\]
Struktur autokorelasi dinyatakan oleh semivariogram:
\[\gamma(h)=\frac{1}{2N(h)}\sum_{i=1}^{N(h)}\big[Z(s_i)-Z(s_i+h)\big]^2\]
Tiga model variogram dicoba (spherical, exponential, gaussian) dan dipilih yang RMSE LOOCV-nya terkecil. Parameter model adalah nugget (variasi pada jarak sangat dekat atau galat pengukuran), partial sill (variasi spasial yang terstruktur), dan range (jarak korelasi). Karena kriging bergantung pada jarak, koordinat diproyeksikan ke sistem koordinat lokal berbasis meter (azimuthal equidistant berpusat di 4,5°LU, 96,5°BT).
Universal Kriging (UK)
UK adalah kriging dengan rata-rata yang mengikuti tren:
\[Z(s)=\mu(s)+\varepsilon(s),\qquad \mu(s)=\beta_0+\beta_1 x+\beta_2 y\]
Tren linear terhadap koordinat (\(x\) dan \(y\)) dipakai karena belum ada kovariat lain pada data.
Validasi Silang Leave-One-Out (LOOCV)
Setiap titik disembunyikan bergiliran, nilainya diduga dari titik lain, lalu dibandingkan dengan nilai sebenarnya:
\[RMSE=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\big(Z(s_i)-\hat{Z}_{-i}(s_i)\big)^2},\quad MAE=\frac{1}{n}\sum_{i=1}^{n}\big|Z(s_i)-\hat{Z}_{-i}(s_i)\big|,\quad R^2=1-\frac{\sum_i (Z_i-\hat{Z}_{-i})^2}{\sum_i (Z_i-\bar{Z})^2}\]
Sebagai pembanding disertakan baseline, yaitu menduga setiap titik dengan rata-rata titik lainnya (tanpa informasi spasial). Metode interpolasi hanya bermanfaat bila lebih akurat daripada baseline. Untuk OK dan UK, model variogram difit dari seluruh data lalu dipakai pada setiap lipatan, sehingga RMSE sedikit lebih optimistis.
Pembatasan Area Prediksi
Prediksi hanya dibuat pada sel grid 3 km yang berjarak paling jauh 50 km dari titik pengamatan terdekat. Di luar radius itu tidak ada data pendukung, sehingga hasil hanya berupa ekstrapolasi dan ditampilkan sebagai area abu-abu.
Hasil dan Pembahasan
Statistik Deskriptif
statistik_DO <- data.frame(
Statistik = c("Jumlah titik", "Mean", "Median", "Standar deviasi", "Minimum",
"Kuartil 1", "Kuartil 3", "Maksimum", "Rentang", "CV (%)",
"Skewness", "Kurtosis (excess)", "Shapiro-Wilk p"),
Nilai = round(c(length(nilai_DO), mean(nilai_DO), median(nilai_DO), sd(nilai_DO),
min(nilai_DO), quantile(nilai_DO, 0.25), quantile(nilai_DO, 0.75),
max(nilai_DO), diff(range(nilai_DO)),
sd(nilai_DO) / mean(nilai_DO) * 100, skewness(nilai_DO),
kurtosis(nilai_DO) - 3, shapiro.test(nilai_DO)$p.value), 3)
)
knitr::kable(statistik_DO, caption = "Statistik deskriptif DO (mg/L)")| Statistik | Nilai |
|---|---|
| Jumlah titik | 78.000 |
| Mean | 5.345 |
| Median | 5.420 |
| Standar deviasi | 0.682 |
| Minimum | 2.665 |
| Kuartil 1 | 4.977 |
| Kuartil 3 | 5.680 |
| Maksimum | 7.140 |
| Rentang | 4.475 |
| CV (%) | 12.752 |
| Skewness | -0.667 |
| Kurtosis (excess) | 2.200 |
| Shapiro-Wilk p | 0.024 |
p_hist <- ggplot(titik_aceh, aes(x = DO)) +
geom_histogram(bins = 10, fill = "steelblue", color = "white") +
labs(title = "Distribusi Oksigen Terlarut", x = "DO (mg/L)", y = "Frekuensi") +
theme_minimal()
print(p_hist)
p_box <- ggplot(titik_aceh, aes(y = DO)) +
geom_boxplot(fill = "lightblue") +
labs(title = "Boxplot Oksigen Terlarut", x = NULL, y = "DO (mg/L)") +
theme_minimal() +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
print(p_box)
# Titik pencilan (di luar 1,5 x IQR)
batas_iqr <- quantile(nilai_DO, c(0.25, 0.75)) + c(-1, 1) * 1.5 * IQR(nilai_DO)
is_pencilan <- nilai_DO < batas_iqr[1] | nilai_DO > batas_iqr[2]
pencilan <- titik_aceh[is_pencilan, ]
# Nama sungai dan lokasi sampling dari titik pencilan (untuk diperiksa)
data_analisis %>%
semi_join(st_drop_geometry(pencilan), by = c("lon", "lat")) %>%
select(Nama_Sungai, Lokasi_Sampling, Tahap, DO_num, lon, lat)## Nama_Sungai
## 1 Danau Aneuk Laot
## 2 Danau Aneuk Laot
## 3 Danau Aneuk Laot
## Lokasi_Sampling
## 1 Titik V - Outlet Danau (Air Alur) Desa Pasiran Kec. Suka Karya
## 2 Titik II - Tengah Danau Lingkungan Putro Ijo Gampong Danau Aneuk Laot Kec. Suka Karya
## 3 Titik V - Outlet Danau (Air Alur) Desa Pasiran Kec. Suka Karya
## Tahap DO_num lon lat
## 1 I 3.27 95.32467 5.881306
## 2 II 7.14 95.33136 5.868639
## 3 II 2.06 95.32467 5.881306Hasil analisis terhadap 78 lokasi pengamatan menunjukkan bahwa rata-rata konsentrasi oksigen terlarut (Dissolved Oxygen/DO) adalah 5,35 mg/L, sedangkan median sebesar 5,42 mg/L. Kedekatan kedua nilai tersebut menunjukkan bahwa pusat distribusi data berada pada kisaran yang relatif serupa. Sebanyak 50% data berada pada rentang kuartil pertama sebesar 4,98 mg/L hingga kuartil ketiga sebesar 5,68 mg/L.
Dari segi penyebaran, standar deviasi sebesar 0,68 mg/L dan koefisien variasi sebesar 12,8% menunjukkan bahwa variasi nilai DO relatif kecil dibandingkan dengan rata-ratanya. Namun, rentang nilai keseluruhan cukup lebar, yaitu dari 2,67 hingga 7,14 mg/L. Hal ini menunjukkan bahwa meskipun sebagian besar lokasi memiliki nilai DO yang berdekatan, terdapat beberapa lokasi dengan nilai yang berbeda cukup jauh dari kelompok utama.
Histogram menunjukkan bahwa nilai DO paling banyak terkonsentrasi pada kisaran 5–6 mg/L. Nilai skewness sebesar −0,67 menunjukkan distribusi yang cenderung menceng ke kiri, yang mengindikasikan adanya sejumlah nilai DO rendah. Sementara itu, kurtosis berlebih sebesar 2,20 menunjukkan bahwa distribusi memiliki ekor yang lebih berat dibandingkan distribusi normal.
Hasil uji Shapiro-Wilk menghasilkan nilai signifikansi sebesar 0,024, lebih kecil daripada 0,05. Dengan demikian, hipotesis bahwa data berasal dari distribusi normal ditolak pada taraf signifikansi 5%. Ketidaknormalan ini dapat berkaitan dengan keberadaan nilai ekstrem, meskipun penyebabnya tidak dapat dipastikan hanya dari hasil pengujian tersebut.
Berdasarkan aturan 1,5 × IQR, terdapat dua nilai pencilan, yaitu 2,665 mg/L dan 7,140 mg/L. Keduanya berkaitan dengan pengamatan di sekitar Danau Aneuk Laot, Sabang. Perbedaan nilai tersebut menunjukkan variasi DO yang cukup besar pada lokasi yang berdekatan. Kondisi ini perlu ditinjau lebih lanjut untuk membedakan variasi lingkungan yang nyata dari kemungkinan kesalahan pencatatan atau pengukuran.
Jika dibandingkan dengan baku mutu DO dalam Lampiran VI PP Nomor 22 Tahun 2021, batas minimum DO untuk air kelas I adalah 6 mg/L, kelas II adalah 4 mg/L, kelas III adalah 3 mg/L, dan kelas IV adalah 1 mg/L.
Rata-rata DO sebesar 5,35 mg/L berada di atas batas minimum kelas II, tetapi di bawah batas minimum kelas I. Namun, nilai rata-rata ini tidak berarti seluruh lokasi memenuhi kelas II. Penilaian harus dilakukan pada setiap lokasi dan mempertimbangkan peruntukan badan air serta parameter mutu air lainnya.
Kesimpulannya, data DO tahun 2020 menunjukkan konsentrasi yang umumnya terkumpul pada kisaran 5–6 mg/L, dengan variasi relatif kecil dibandingkan rata-ratanya, tetapi terdapat pencilan dan distribusi yang tidak normal. Sebelum interpolasi Kriging dilakukan, pencilan perlu diperiksa dan struktur spasial data dievaluasi melalui variogram serta validasi silang. Hasil tersebut kemudian dapat digunakan untuk menggambarkan pola persebaran DO di wilayah Aceh secara lebih representatif.
Peta Sebaran Titik
peta_DO <- ggplot() +
geom_sf(data = batas_aceh, fill = "grey95", color = "grey40", linewidth = 0.4) +
geom_sf(data = titik_aceh, aes(color = DO), size = 2.5, alpha = 0.9) +
scale_color_gradientn(
colors = c("#FFC200", "#FF6B00", "#E31A1C", "#800080", "#00207E"),
name = "DO (mg/L)"
) +
labs(title = "Sebaran Oksigen Terlarut di Aceh",
subtitle = "Data kualitas air tahun 2020",
x = "Bujur Timur", y = "Lintang Utara",
caption = "Sumber: Satu Data Indonesia (data.go.id), Pemerintah Provinsi Aceh") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(peta_DO)
Peta menunjukkan persebaran konsentrasi oksigen terlarut (Dissolved Oxygen/DO) pada titik pengamatan kualitas air di Provinsi Aceh tahun 2020. Warna pada setiap titik merepresentasikan nilai DO dalam satuan mg/L, dengan warna oranye hingga merah menunjukkan konsentrasi relatif rendah, sedangkan warna ungu hingga biru menunjukkan konsentrasi relatif tinggi.
Sebaran stasiun pengamatan tidak merata. Titik pengamatan terkonsentrasi di wilayah utara, barat, dan tengah Aceh. Beberapa titik membentuk kelompok atau deret yang kemungkinan berkaitan dengan pengambilan sampel berulang pada sungai atau lokasi yang sama. Sementara itu, wilayah Aceh bagian selatan dan tenggara serta Pulau Simeulue tidak tampak memiliki titik pengamatan pada peta. Kondisi ini menunjukkan bahwa cakupan data belum mewakili seluruh wilayah Aceh secara merata.
Pola nilai DO menunjukkan variasi antar lokasi. Nilai yang relatif rendah, sekitar 4–5 mg/L, ditunjukkan oleh sejumlah titik berwarna oranye hingga merah, terutama di bagian barat laut dan beberapa lokasi pesisir. Sementara itu, nilai yang relatif tinggi, sekitar 6 mg/L atau lebih, ditunjukkan oleh titik berwarna ungu, antara lain di wilayah utara-tengah, beberapa lokasi sekitar 97° BT pada lintang 3,6–4,6° LU, serta sebagian titik di pesisir barat. Perbedaan warna pada titik yang berdekatan mengindikasikan bahwa konsentrasi DO dapat bervariasi secara lokal.
Kesesuaian dengan Hukum Tobler bersifat parsial. Hukum Tobler menyatakan bahwa objek yang berdekatan cenderung memiliki hubungan yang lebih kuat dibandingkan objek yang berjauhan. Pada peta ini, beberapa kelompok titik yang berdekatan menunjukkan nilai DO yang relatif serupa. Namun, terdapat pula titik berdekatan dengan nilai yang berbeda cukup besar, termasuk pencilan yang telah diidentifikasi sebelumnya serta titik oranye dan ungu di sekitar 96,9° BT dan 4,6° LU. Hal ini menunjukkan bahwa kedekatan geografis saja belum cukup untuk menjelaskan variasi DO. Kondisi lingkungan setempat, karakteristik aliran, dan aktivitas di sekitar perairan juga mungkin berpengaruh, meskipun faktor-faktor tersebut belum diuji secara langsung dalam analisis ini.
Implikasi terhadap interpolasi spasial. Ketidakmerataan titik pengamatan dan adanya perbedaan nilai pada lokasi berdekatan perlu diperhatikan dalam pemodelan interpolasi. Prediksi umumnya lebih dapat dipertanggungjawabkan di sekitar lokasi yang memiliki dukungan data memadai, sedangkan wilayah yang jauh dari titik pengamatan memiliki ketidakpastian lebih besar. Pembatasan peta interpolasi hingga radius 50 km dari titik terdekat dapat digunakan sebagai langkah awal untuk menghindari prediksi yang terlalu jauh dari data, tetapi radius tersebut perlu didukung oleh evaluasi variogram dan validasi silang. Selain itu, kelompok titik yang berdekatan dapat mengurangi jumlah informasi spasial yang benar-benar independen.
Kesimpulan: Sebaran DO di Aceh pada tahun 2020 menunjukkan variasi spasial, baik antarwilayah maupun pada beberapa lokasi yang berdekatan. Distribusi stasiun yang tidak merata membatasi keterwakilan data, sehingga hasil interpolasi perlu ditafsirkan secara hati-hati, khususnya pada wilayah yang minim pengamatan. Peta ini dapat menjadi dasar untuk mengidentifikasi lokasi yang memerlukan pemantauan tambahan, tetapi belum cukup untuk menyimpulkan kualitas air seluruh Provinsi Aceh.
Persiapan Interpolasi
# Proyeksi lokal dalam meter
crs_lokal <- st_crs("+proj=aeqd +lat_0=4.5 +lon_0=96.5 +datum=WGS84 +units=m +no_defs")
titik_lokal <- st_transform(titik_aceh, crs_lokal)
batas_lokal <- st_transform(batas_aceh, crs_lokal)
# Koordinat dalam km, dipakai sebagai tren pada Universal Kriging
xy <- st_coordinates(titik_lokal)
titik_lokal$x_km <- xy[, 1] / 1000
titik_lokal$y_km <- xy[, 2] / 1000
data_sp <- as(titik_lokal, "Spatial")
observasi <- data_sp$DO
# Grid prediksi
ukuran_sel <- 3000
grid <- st_make_grid(st_union(batas_lokal), cellsize = ukuran_sel, what = "centers")
grid_sf <- st_sf(geometry = grid) %>% st_filter(st_union(batas_lokal))
# Hanya prediksi di sekitar titik pengamatan (hindari ekstrapolasi jauh)
batas_jarak <- 50000 # meter
jarak_min <- apply(st_distance(grid_sf, titik_lokal), 1, min)
grid_sf <- grid_sf[as.numeric(jarak_min) <= batas_jarak, ]
xy_grid <- st_coordinates(grid_sf)
grid_sf$x_km <- xy_grid[, 1] / 1000
grid_sf$y_km <- xy_grid[, 2] / 1000
grid_sp <- as(grid_sf, "Spatial")Interpolasi spasial dilakukan menggunakan 78 titik pengamatan konsentrasi oksigen terlarut (DO) yang berada di dalam batas Provinsi Aceh. Prediksi dibuat pada grid dengan ukuran sel 3 km × 3 km, dengan jumlah sel prediksi sebanyak nrow(grid_sf) sesuai hasil penyaringan grid di dalam wilayah Aceh dan radius maksimum 50 km dari titik pengamatan terdekat. Pembatasan ini bertujuan agar prediksi lebih terfokus pada wilayah yang memiliki dukungan data pengamatan. Jumlah sel grid menunjukkan banyaknya lokasi prediksi, bukan jumlah titik pengukuran aktual.
Interpolasi IDW
# Uji semua kandidat power dalam satu proses agar pemilihan konsisten
powers_semua <- c(0.25, 0.50, 0.75, 1, 1.5, 2, 2.5, 3, 4)
rmse_semua <- sapply(powers_semua, function(pw) {
cv <- krige.cv(DO ~ 1, data_sp, nfold = n,
set = list(idp = pw), verbose = FALSE)
sqrt(mean(cv$residual^2, na.rm = TRUE))
})
tabel_power <- data.frame(power = powers_semua,
RMSE_LOOCV = round(rmse_semua, 4))
knitr::kable(tabel_power,
caption = "Perbandingan RMSE LOOCV untuk Pemilihan Power IDW")| power | RMSE_LOOCV |
|---|---|
| 0.25 | 0.6755 |
| 0.50 | 0.6738 |
| 0.75 | 0.6821 |
| 1.00 | 0.6869 |
| 1.50 | 0.6912 |
| 2.00 | 0.6994 |
| 2.50 | 0.7090 |
| 3.00 | 0.7176 |
| 4.00 | 0.7296 |
p_best <- powers_semua[which.min(rmse_semua)]
rmse_idw_cv_best <- min(rmse_semua, na.rm = TRUE)
cat("Power IDW terbaik:", p_best, "\n")## Power IDW terbaik: 0.5cat("RMSE LOOCV terbaik:", round(rmse_idw_cv_best, 4), "mg/L\n")## RMSE LOOCV terbaik: 0.6738 mg/LPemilihan parameter IDW dilakukan dengan membandingkan sembilan kandidat nilai power menggunakan LOOCV. Parameter yang dipilih adalah nilai power dengan RMSE terkecil dari seluruh kandidat, bukan hanya dari rentang awal. Nilai parameter dan RMSE terbaik ditampilkan pada tabel hasil di atas. Pemilihan ini kemudian digunakan untuk membuat peta IDW dan menghitung ulang perbandingan model.
hasil_idw <- gstat::idw(DO ~ 1, locations = data_sp, newdata = grid_sp,
idp = p_best, debug.level = 0)
hasil_idw_sf <- st_as_sf(hasil_idw)Peta IDW pada bagian ini dibuat menggunakan nilai power yang memberikan RMSE LOOCV terendah dari seluruh kandidat yang diuji. Hasil tabel pemilihan power menjadi acuan, dan nilai tersebut digunakan kembali dalam prediksi grid serta evaluasi LOOCV akhir. Pembobotan jarak yang lebih merata dapat menghasilkan prediksi yang lebih halus, tetapi hasil ini saja belum cukup untuk menyimpulkan bahwa struktur spasial DO lemah; interpretasi tetap perlu mempertimbangkan variogram, validasi silang, dan sebaran titik pengamatan.
Ordinary Kriging
n_peringatan <- 0
# Fit tiga model variogram dengan nilai awal yang realistis;
# peringatan "tidak konvergen" dihitung agar dapat dilaporkan
fit_variogram_3 <- function(v, nilai) {
vr <- var(nilai)
fits <- lapply(c("Sph", "Exp", "Gau"), function(m) {
withCallingHandlers(
fit.variogram(v, vgm(psill = 0.5 * vr, model = m,
range = 30000, nugget = 0.5 * vr)),
warning = function(w) {
n_peringatan <<- n_peringatan + 1
invokeRestart("muffleWarning")
}
)
})
names(fits) <- c("Sph", "Exp", "Gau")
fits
}v_ok <- variogram(DO ~ 1, data_sp, cutoff = 100000, width = 10000)
knitr::kable(round(v_ok[, c("np", "dist", "gamma")], 3),
caption = "Semivariogram empiris DO (jarak dalam meter)")| np | dist | gamma |
|---|---|---|
| 125 | 4473.50 | 0.498 |
| 40 | 14422.54 | 0.246 |
| 70 | 25686.78 | 0.254 |
| 91 | 34789.42 | 0.628 |
| 77 | 45241.92 | 0.225 |
| 107 | 54777.81 | 0.417 |
| 158 | 65700.03 | 0.394 |
| 114 | 74352.13 | 0.283 |
| 109 | 85691.65 | 0.312 |
| 119 | 94863.43 | 0.353 |
plot(v_ok, main = "Semivariogram Empiris DO")
n_peringatan <- 0
fits_ok <- fit_variogram_3(v_ok, observasi)
n_warn_ok <- n_peringatan
rmse_model_ok <- sapply(fits_ok, function(f) {
cv <- krige.cv(DO ~ 1, data_sp, f, nfold = n, verbose = FALSE)
sqrt(mean(cv$residual^2))
})
print(round(rmse_model_ok, 4))## Sph Exp Gau
## 0.6656 0.6593 0.6760nama_ok <- names(which.min(rmse_model_ok))
model_ok <- fits_ok[[nama_ok]]
print(model_ok)## model psill range
## 1 Nug 0.3173826 0.0
## 2 Exp 0.6519559 390438.2plot(v_ok, model_ok, main = paste("Semivariogram OK - model", nama_ok))
hasil_ok <- krige(DO ~ 1, locations = data_sp, newdata = grid_sp,
model = model_ok, debug.level = 0)
hasil_ok_sf <- st_as_sf(hasil_ok)Variogram empiris. Pada jarak terdekat sekitar 4,5 km, nilai semivarian sebesar 0,498, mendekati varians data sebesar 0,465. Semivarian kemudian turun menjadi 0,246 pada jarak 14,4 km, meningkat hingga 0,628 pada jarak 34,8 km, lalu berfluktuasi tanpa pola kenaikan yang konsisten. Hal ini menunjukkan bahwa hubungan spasial DO belum terlihat jelas.
Model terpilih. Model Exponential memiliki RMSE LOOCV terendah, yaitu 0,6593 mg/L, dibandingkan model Spherical sebesar 0,6656 mg/L dan Gaussian sebesar 0,6760 mg/L. Parameter model meliputi nugget 0,317, partial sill 0,652, dan range sekitar 390,44 km.
Catatan tentang model. Garis model variogram tidak mengikuti pola empiris secara konsisten. Nilai range yang jauh melampaui wilayah studi juga menunjukkan bahwa parameter tersebut perlu ditafsirkan secara hati-hati. Dengan demikian, meskipun model Exponential menghasilkan RMSE terendah, struktur spasial DO belum teridentifikasi dengan kuat. Hasil Kriging sebaiknya diperlakukan sebagai estimasi awal dan dievaluasi lebih lanjut sebelum digunakan untuk menarik kesimpulan spasial.
Universal Kriging
v_uk <- variogram(DO ~ x_km + y_km, data_sp, cutoff = 100000, width = 10000)
plot(v_uk, main = "Semivariogram Residual UK")
n_peringatan <- 0
fits_uk <- fit_variogram_3(v_uk, observasi)
n_warn_uk <- n_peringatan
rmse_model_uk <- sapply(fits_uk, function(f) {
cv <- krige.cv(DO ~ x_km + y_km, data_sp, f, nfold = n, verbose = FALSE)
sqrt(mean(cv$residual^2))
})
print(round(rmse_model_uk, 4))## Sph Exp Gau
## 0.6682 0.6623 0.6761nama_uk <- names(which.min(rmse_model_uk))
model_uk <- fits_uk[[nama_uk]]
print(model_uk)## model psill range
## 1 Nug 0.3185835 0.0
## 2 Exp 0.6833149 405356.6plot(v_uk, model_uk, main = paste("Semivariogram UK - model", nama_uk))
hasil_uk <- krige(DO ~ x_km + y_km, locations = data_sp, newdata = grid_sp,
model = model_uk, debug.level = 0)
hasil_uk_sf <- st_as_sf(hasil_uk)Semivariogram residual menunjukkan pola yang berfluktuasi dan tidak meningkat secara konsisten terhadap jarak. Berdasarkan validasi silang LOOCV, model Exponential menghasilkan RMSE terendah sebesar 0,6623 mg/L, dibandingkan model Spherical (0,6682 mg/L) dan Gaussian (0,6761 mg/L). Oleh karena itu, model Exponential dipilih untuk Universal Kriging dengan tren koordinat \(x\) dan \(y\).
Model menghasilkan nugget sebesar 0,3186, partial sill sebesar 0,6833, dan range sekitar 405,36 km. Namun, garis model belum mengikuti fluktuasi semivariogram empiris dengan baik dan nilai range jauh melampaui wilayah studi. Dengan demikian, struktur spasial residual DO belum terlihat kuat, sehingga hasil prediksi UK perlu ditafsirkan secara hati-hati dan dievaluasi melalui validasi silang sebelum digunakan untuk menarik kesimpulan.
Peta Hasil Interpolasi
# Skala warna SAMA untuk IDW, OK, dan UK agar bisa dibandingkan
rentang_pred <- range(c(hasil_idw_sf$var1.pred,
hasil_ok_sf$var1.pred,
hasil_uk_sf$var1.pred))
keterangan <- paste0("Area abu-abu: lebih dari ", batas_jarak / 1000,
" km dari titik pengamatan\nterdekat (tidak diprediksi)")
buat_peta <- function(hasil, judul) {
xy_hasil <- st_coordinates(hasil)
hasil$x <- xy_hasil[, 1]
hasil$y <- xy_hasil[, 2]
ggplot() +
geom_sf(data = batas_lokal, fill = "grey88", color = NA) +
geom_tile(data = hasil, aes(x = x, y = y, fill = var1.pred),
width = ukuran_sel, height = ukuran_sel) +
geom_sf(data = batas_lokal, fill = NA, color = "grey25", linewidth = 0.4) +
geom_sf(data = titik_lokal, shape = 21, fill = "white", color = "black",
size = 2.2, stroke = 0.7) +
scale_fill_gradientn(
colors = c("#FFC200", "#FF6B00", "#E31A1C", "#800080", "#00207E"),
limits = rentang_pred,
name = "DO prediksi\n(mg/L)"
) +
coord_sf(crs = st_crs(batas_lokal), datum = NA, expand = FALSE) +
labs(title = judul, subtitle = "Oksigen terlarut (DO), tahun 2020",
caption = keterangan,
x = "Koordinat X (meter)", y = "Koordinat Y (meter)") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "grey40"),
plot.caption = element_text(color = "grey40", hjust = 0),
plot.caption.position = "plot")
}peta_idw <- buat_peta(hasil_idw_sf, paste0("Interpolasi IDW (p = ", p_best, ") - DO Aceh"))
print(peta_idw)
peta_ok <- buat_peta(hasil_ok_sf, paste0("Ordinary Kriging (", nama_ok, ") - DO Aceh"))
print(peta_ok)
peta_uk <- buat_peta(hasil_uk_sf, paste0("Universal Kriging (", nama_uk, ") - DO Aceh"))
print(peta_uk)
xy_v <- st_coordinates(hasil_ok_sf)
df_v <- data.frame(x = xy_v[, 1], y = xy_v[, 2], varians = hasil_ok_sf$var1.var)
peta_var <- ggplot() +
geom_sf(data = batas_lokal, fill = "grey88", color = NA) +
geom_tile(data = df_v, aes(x = x, y = y, fill = varians),
width = ukuran_sel, height = ukuran_sel) +
geom_sf(data = batas_lokal, fill = NA, color = "grey25", linewidth = 0.4) +
geom_sf(data = titik_lokal, shape = 21, fill = "white", color = "black", size = 2.2) +
scale_fill_gradient(low = "#FFC200", high = "#00207E", name = "Varians") +
coord_sf(crs = st_crs(batas_lokal), datum = NA, expand = FALSE) +
labs(title = "Varians Prediksi Ordinary Kriging - DO Aceh",
caption = keterangan,
x = "Koordinat X (meter)", y = "Koordinat Y (meter)") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", hjust = 0.5),
plot.caption = element_text(color = "grey40", hjust = 0),
plot.caption.position = "plot")
print(peta_var)
Interpretasi peta hasil interpolasi DO di Aceh tahun 2020
- Interpolasi IDW (power = 0.5)
Peta IDW menunjukkan pola prediksi konsentrasi DO pada area yang didukung titik pengamatan. Terdapat beberapa area dengan nilai lebih tinggi dan lebih rendah di sekitar titik pengamatan tertentu. Metode IDW menghasilkan permukaan prediksi berdasarkan pembobotan jarak; hasilnya sensitif terhadap nilai pengamatan dan parameter power terpilih, yaitu 0.5.
- Ordinary Kriging (model Exponential)
Peta Ordinary Kriging menunjukkan pola DO yang relatif lebih tinggi di sebagian wilayah tengah hingga selatan area prediksi, sedangkan beberapa wilayah barat memiliki nilai yang lebih rendah. Model ini menggunakan hubungan spasial yang direpresentasikan oleh variogram Exponential. Namun, karena variogram empiris tidak menunjukkan pola yang konsisten dan nilai range hasil pemodelan jauh melampaui wilayah studi, pola prediksi tersebut perlu ditafsirkan secara hati-hati.
- Universal Kriging (model Exponential)
Peta Universal Kriging menunjukkan pola yang menyerupai Ordinary Kriging, tetapi memperhitungkan tren berdasarkan koordinat X dan Y. Nilai prediksi tampak meningkat pada sebagian wilayah selatan dibandingkan wilayah barat. Perbedaan ini menunjukkan pengaruh komponen tren dalam model, tetapi belum cukup untuk menyimpulkan adanya pola peningkatan DO yang nyata secara geografis.
- Varians prediksi Ordinary Kriging
Peta varians menunjukkan tingkat ketidakpastian prediksi model Ordinary Kriging. Area dengan warna lebih gelap memiliki varians prediksi lebih tinggi, sedangkan area dengan warna lebih terang memiliki varians lebih rendah. Peningkatan varians di beberapa bagian tepi area prediksi menunjukkan bahwa ketidakpastian model dapat berbeda menurut lokasi. Varians ini merupakan ukuran ketidakpastian model, bukan ukuran langsung kualitas air.
Kesimpulan: Ketiga peta memperlihatkan pola prediksi DO yang berbeda, meskipun terdapat kemiripan pada beberapa wilayah. IDW menghasilkan prediksi berdasarkan jarak, sedangkan Ordinary Kriging mempertimbangkan struktur spasial melalui variogram dan Universal Kriging menambahkan komponen tren koordinat. Karena struktur spasial data belum teridentifikasi dengan kuat, pemilihan model sebaiknya didasarkan pada hasil validasi silang, bukan hanya tampilan peta.
Perlu dicatat bahwa area abu-abu menunjukkan wilayah yang tidak diprediksi karena berjarak lebih dari 50 km dari titik pengamatan terdekat. Oleh sebab itu, peta tidak dapat digunakan untuk menyimpulkan kondisi DO di seluruh wilayah Aceh, terutama pada area tanpa dukungan data pengamatan.
Perbandingan Model dengan LOOCV
# IDW
pred_idw <- rep(NA_real_, n)
for (i in seq_len(n)) {
latih <- data_sp[-i, ]
uji <- data_sp[i, ]
pred_idw[i] <- gstat::idw(DO ~ 1, locations = latih, newdata = uji,
idp = p_best, debug.level = 0)$var1.pred
}
# OK dan UK
cv_ok <- krige.cv(DO ~ 1, locations = data_sp, model = model_ok,
nfold = n, verbose = FALSE)
cv_uk <- krige.cv(DO ~ x_km + y_km, locations = data_sp, model = model_uk,
nfold = n, verbose = FALSE)
pred_ok <- cv_ok$var1.pred
pred_uk <- cv_uk$var1.pred
# Baseline: menebak dengan rata-rata (tanpa informasi spasial)
pred_mean <- sapply(seq_len(n), function(i) mean(observasi[-i]))
hitung_metrik <- function(aktual, prediksi, metode) {
valid <- is.finite(aktual) & is.finite(prediksi)
e <- aktual[valid] - prediksi[valid]
data.frame(
Metode = metode,
Jumlah_valid = sum(valid),
RMSE = sqrt(mean(e^2)),
MAE = mean(abs(e)),
R2 = 1 - sum(e^2) / sum((aktual[valid] - mean(aktual[valid]))^2)
)
}
hasil_loocv <- bind_rows(
hitung_metrik(observasi, pred_idw, paste0("IDW (p = ", p_best, ")")),
hitung_metrik(observasi, pred_ok, paste0("Ordinary Kriging (", nama_ok, ")")),
hitung_metrik(observasi, pred_uk, paste0("Universal Kriging (", nama_uk, ")")),
hitung_metrik(observasi, pred_mean, "Baseline: rata-rata global")
) %>%
mutate(RMSE = round(RMSE, 4), MAE = round(MAE, 4), R2 = round(R2, 3)) %>%
arrange(RMSE)
knitr::kable(hasil_loocv, caption = "Hasil validasi silang Leave-One-Out")| Metode | Jumlah_valid | RMSE | MAE | R2 |
|---|---|---|---|---|
| Ordinary Kriging (Exp) | 78 | 0.6593 | 0.4783 | 0.052 |
| Universal Kriging (Exp) | 78 | 0.6623 | 0.4802 | 0.044 |
| IDW (p = 0.5) | 78 | 0.6738 | 0.5070 | 0.010 |
| Baseline: rata-rata global | 78 | 0.6860 | 0.5110 | -0.026 |
par(mfrow = c(1, 3))
for (m in list(list("IDW", pred_idw), list("Ordinary Kriging", pred_ok),
list("Universal Kriging", pred_uk))) {
plot(observasi, m[[2]], pch = 19, col = "#00207E",
xlab = "Observasi", ylab = "Prediksi LOOCV",
main = paste(m[[1]], "- RMSE =", round(sqrt(mean((observasi - m[[2]])^2)), 3)))
abline(0, 1, lty = 2)
}
par(mfrow = c(1, 1))# Apakah OK benar-benar lebih akurat daripada menebak rata-rata?
e2_ok <- (observasi - pred_ok)^2
e2_mean <- (observasi - pred_mean)^2
uji_wil <- suppressWarnings(
wilcox.test(e2_mean, e2_ok, paired = TRUE, alternative = "greater", exact = FALSE)
)
uji_wil##
## Wilcoxon signed rank test with continuity correction
##
## data: e2_mean and e2_ok
## V = 1747, p-value = 0.1524
## alternative hypothesis: true location shift is greater than 0Interpretasi hasil validasi model interpolasi
Tabel LOOCV membandingkan IDW dengan nilai power terpilih, Ordinary Kriging, Universal Kriging, dan baseline rata-rata global. Model dengan RMSE dan MAE lebih kecil memiliki kesalahan prediksi yang lebih rendah. Nilai \(R^2\) yang lebih tinggi menunjukkan prediksi yang lebih baik relatif terhadap variasi nilai observasi, tetapi nilai negatif berarti model tidak lebih baik daripada prediksi menggunakan rata-rata observasi dalam ukuran tersebut.
Penetapan metode terbaik harus mengikuti tabel LOOCV yang dihasilkan setelah seluruh kode dijalankan ulang. Perhatikan juga apakah model interpolasi mengungguli baseline. Jika selisih RMSE kecil atau uji berpasangan tidak signifikan, hasil tidak boleh ditafsirkan sebagai bukti kuat bahwa model tersebut lebih unggul. Hasil validasi merupakan evaluasi pada 78 lokasi yang tersedia dan tidak otomatis menjamin ketepatan prediksi di wilayah Aceh yang tidak terwakili oleh titik pengamatan.
Plot observasi dan prediksi digunakan untuk memeriksa apakah prediksi mengikuti rentang nilai aktual atau cenderung mengumpul di sekitar rata-rata. Nilai ekstrem dapat sulit diprediksi, terutama ketika titik pengamatan jarang atau nilai pada lokasi yang berdekatan sangat berbeda. Hal ini perlu dipertimbangkan bersama hasil variogram dan sebaran titik.
Untuk Ordinary dan Universal Kriging, model variogram pada skrip ini dipilih menggunakan LOOCV, sementara variogram dipasang terlebih dahulu menggunakan seluruh data. Oleh karena itu, evaluasi tersebut dapat sedikit optimistis. Pengujian yang lebih ketat dapat dilakukan dengan memasang ulang variogram di setiap lipatan validasi silang.
Keterbatasan
- Struktur spasial belum teridentifikasi dengan kuat. Variogram empiris berfluktuasi dan model terpilih memiliki range jauh melampaui wilayah studi. Jumlah peringatan saat fitting ditampilkan sebagai informasi tambahan (OK 3, UK 3), tetapi jumlah peringatan saja tidak menjadi bukti tunggal ketidakstabilan model.
- Sebaran titik tidak merata. Aceh bagian selatan dan tenggara serta Pulau Simeulue tidak terwakili, dan banyak titik berkelompok pada sungai yang sama. Hasil hanya andal di sekitar titik pengamatan.
- Pencilan dan kualitas data. Dua titik pencilan yang berdekatan memengaruhi variogram, dan satu lokasi dikeluarkan karena koordinatnya berada di luar batas. Metadata sumber juga belum sepenuhnya memenuhi Prinsip SDI.
- Pemilihan model dan LOOCV. Power IDW dan model variogram dipilih memakai LOOCV yang sama dengan evaluasi akhir, dan variogram difit dari seluruh data, sehingga RMSE sedikit optimistis. Dengan n = 78, selisih kecil antar metode belum tentu bermakna.
- Nilai ber-tanda “<” pada DO diperlakukan sebagai angka batas yang tertera, yang merupakan penyederhanaan.
Kesimpulan
Data oksigen terlarut (Dissolved Oxygen/DO) pada sungai dan danau di Aceh tahun 2020 yang dianalisis terdiri atas 78 lokasi pengamatan. Rata-rata konsentrasi DO sebesar 5,35 mg/L, dengan rentang 2,665–7,14 mg/L. Distribusi data tidak normal berdasarkan uji Shapiro–Wilk (\(p=0,024\)), dengan koefisien variasi sebesar 12,75%. Terdapat dua pencilan berdasarkan aturan 1,5 × IQR yang perlu diperiksa lebih lanjut.
Sebaran titik pengamatan tidak merata dan terkonsentrasi di wilayah utara, barat, dan tengah Aceh. Konsentrasi DO bervariasi antar lokasi, termasuk pada beberapa titik yang berdekatan. Pola tersebut hanya sebagian sesuai dengan Hukum Tobler, yang menyatakan bahwa lokasi berdekatan cenderung memiliki hubungan spasial lebih kuat.
Berdasarkan hasil validasi silang Leave-One-Out Cross-Validation (LOOCV), metode Ordinary Kriging dengan model variogram eksponensial menghasilkan RMSE terendah sebesar 0,6593 dan MAE sebesar 0,4783, dengan nilai R² sebesar 0,052. Universal Kriging menghasilkan RMSE sebesar 0,6623, sedangkan IDW dengan power 0,5 menghasilkan RMSE sebesar 0,6738. Dengan demikian, Ordinary Kriging merupakan metode dengan kinerja prediksi terbaik di antara metode yang dibandingkan berdasarkan nilai RMSE.
Meskipun Ordinary Kriging menghasilkan RMSE terendah, nilai R² yang rendah menunjukkan bahwa kemampuan model dalam menjelaskan variasi konsentrasi DO masih terbatas. Uji Wilcoxon antara Ordinary Kriging dan baseline rata-rata global menghasilkan nilai p sebesar 0,1524, sehingga belum terdapat bukti perbedaan yang signifikan secara statistik pada taraf 5%. Selain itu, variogram tidak menunjukkan struktur spasial yang kuat dan titik pengamatan tersebar tidak merata. Oleh karena itu, peta interpolasi sebaiknya digunakan untuk menggambarkan kecenderungan pola spasial secara umum, bukan sebagai prediksi yang sepenuhnya akurat di seluruh wilayah Aceh.
Penelitian selanjutnya disarankan untuk menambah titik pemantauan, terutama di wilayah Aceh bagian selatan dan timur, serta memeriksa kembali pencilan dan koordinat yang tidak sesuai. Penelitian juga dapat mempertimbangkan variabel tambahan, seperti suhu air, elevasi, dan jarak dari muara, serta menggunakan data dari beberapa tahun pengamatan untuk mengevaluasi konsistensi pola spasial DO.
Daftar Pustaka
Pemerintah Republik Indonesia. (2021). Peraturan Pemerintah Nomor 22 Tahun 2021 tentang Penyelenggaraan Perlindungan dan Pengelolaan Lingkungan Hidup.
Shepard, D. (1968). A two-dimensional interpolation function for irregularly-spaced data. Proceedings of the 1968 ACM National Conference, 517–524.
Tobler, W. R. (1970). A computer movie simulating urban growth in the Detroit region. Economic Geography, 46(2), 234–240.
Moraga, P. (2023). Spatial Statistics for Data Science: Theory and Practice with R. Chapman & Hall/CRC.
Satu Data Indonesia. (t.t.). Parameter Kualitas Air Sungai dan Danau. Pemerintah Provinsi Aceh. https://data.go.id/dataset/dataset/parameter-kualitas-air-sungai-dan-danau.
Pebesma, E. (2018). Simple features for R: Standardized support for spatial vector data. The R Journal, 10(1), 439–446.
Pebesma, E. J. (2004). Multivariable geostatistics in S: The gstat package. Computers & Geosciences, 30(7), 683–691.