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

  1. Menyajikan data numerik lingkungan bertipe titik (DO) di Provinsi Aceh beserta sumbernya.
  2. Meringkas karakteristik data melalui statistik deskriptif.
  3. Memvisualisasikan sebaran spasial DO dalam bentuk peta.
  4. 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:

  1. Koordinat berformat derajat-menit-detik (DMS) dikonversi menjadi derajat desimal. Format yang tidak lazim, seperti 0950, diperbaiki menjadi 095.
  2. Nilai DO yang bertanda < dibaca sebagai angka batas yang tertera (penyederhanaan).
  3. 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.
  4. 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 34
dms_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.881
summary(data_air$lon)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   95.32   95.72   96.41   96.43   96.98   97.90
data_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                 2
nilai_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 \(\{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 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.881306

Hasil 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

powers <- c(1, 1.5, 2, 2.5, 3, 4)
rmse_p <- sapply(powers, function(pw) {
  cv <- krige.cv(DO ~ 1, data_sp, nfold = n, set = list(idp = pw), verbose = FALSE)
  sqrt(mean(cv$residual^2))
})
tabel_power <- data.frame(power = powers, RMSE_LOOCV = round(rmse_p, 4))
knitr::kable(tabel_power, caption = "RMSE LOOCV IDW untuk berbagai nilai power")
RMSE LOOCV IDW untuk berbagai nilai power
power RMSE_LOOCV
1.0 0.6869
1.5 0.6912
2.0 0.6994
2.5 0.7090
3.0 0.7176
4.0 0.7296 
p_best <- powers[which.min(rmse_p)]
cat("Power IDW terbaik:", p_best, "\n")
## Power IDW terbaik: 1

Berdasarkan hasil pengujian menggunakan metode Leave-One-Out Cross-Validation (LOOCV), diperoleh nilai RMSE yang berbeda untuk setiap nilai power pada metode Inverse Distance Weighting (IDW).

Nilai RMSE terendah diperoleh pada power = 1, yaitu sebesar 0,6869 mg/L. Sementara itu, nilai RMSE tertinggi diperoleh pada power = 4, yaitu sebesar 0,7296 mg/L.

Hasil tersebut menunjukkan bahwa power = 1 merupakan parameter terbaik di antara nilai power yang diuji, karena menghasilkan kesalahan prediksi paling kecil. Semakin besar nilai power yang digunakan, RMSE cenderung meningkat. Hal ini menunjukkan bahwa pada data DO di Aceh tahun 2020, pembobotan jarak yang lebih rendah memberikan hasil prediksi yang lebih baik dibandingkan pembobotan yang lebih tinggi dalam pengujian ini.

Secara matematis, IDW dengan power = 1 memberikan bobot berdasarkan kebalikan jarak, sehingga pengaruh titik pengamatan terdekat tidak terlalu dominan dibandingkan pada power yang lebih besar.

powers_kecil <- c(0.25, 0.5, 0.75)
rmse_p_kecil <- sapply(powers_kecil, function(pw) {
  cv <- krige.cv(DO ~ 1, data_sp, nfold = n, set = list(idp = pw), verbose = FALSE)
  sqrt(mean(cv$residual^2))
})
knitr::kable(data.frame(power = powers_kecil, RMSE_LOOCV = round(rmse_p_kecil, 4)),
             caption = "RMSE LOOCV IDW pada power yang lebih kecil (pemeriksaan tambahan)")
RMSE LOOCV IDW pada power yang lebih kecil (pemeriksaan tambahan)
power RMSE_LOOCV
0.25 0.6755
0.50 0.6738
0.75 0.6821

Berdasarkan hasil pemeriksaan tambahan menggunakan metode Leave-One-Out Cross-Validation (LOOCV), diperoleh nilai RMSE sebesar 0,6755 mg/L pada power = 0,25, sebesar 0,6738 mg/L pada power = 0,50, dan sebesar 0,6821 mg/L pada power = 0,75.

Nilai RMSE terendah diperoleh pada power = 0,50, yaitu sebesar 0,6738 mg/L. Nilai ini lebih rendah dibandingkan RMSE pada power = 1, yaitu 0,6869 mg/L, sehingga menunjukkan bahwa penurunan nilai power dapat meningkatkan ketepatan prediksi pada data yang dianalisis.

Dengan demikian, power = 0,50 menjadi parameter IDW terbaik sementara dari seluruh nilai yang telah diuji, yaitu 0,25; 0,50; 0,75; 1; 1,5; 2; 2,5; 3; dan 4. Hasil ini menunjukkan bahwa pembobotan jarak yang lebih landai menghasilkan kesalahan prediksi lebih kecil pada data DO di Aceh tahun 2020.

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)

Berdasarkan hasil Leave-One-Out Cross-Validation (LOOCV), nilai RMSE IDW meningkat dari 0,6869 mg/L pada power = 1 menjadi 0,7296 mg/L pada power = 4. Hasil ini menunjukkan bahwa peningkatan nilai power pada rentang tersebut cenderung menurunkan ketepatan prediksi.

Namun, pemeriksaan tambahan pada nilai power yang lebih kecil menghasilkan RMSE minimum sebesar 0,6738 mg/L pada power = 0,50, lebih rendah dibandingkan RMSE minimum pada pengujian awal. Dengan demikian, power = 0,50 merupakan parameter terbaik sementara dari seluruh nilai yang telah diuji, bukan power = 1.

Kecenderungan bahwa nilai power yang lebih rendah menghasilkan RMSE lebih kecil menunjukkan bahwa pembobotan jarak yang lebih merata memberikan hasil prediksi lebih baik pada data ini. Hal tersebut dapat mengindikasikan bahwa pengaruh jarak terhadap kemiripan nilai DO relatif lemah, tetapi belum cukup untuk menyimpulkan bahwa struktur spasial DO memang lemah. Kesimpulan tersebut perlu didukung oleh analisis autokorelasi spasial, variogram, atau evaluasi model tambahan.

Kesimpulan: berdasarkan hasil LOOCV, IDW dengan power = 0,50 memberikan kesalahan prediksi terendah, yaitu 0,6738 mg/L. Parameter ini dapat digunakan untuk interpolasi selanjutnya, dengan tetap mempertimbangkan hasil validasi dan karakteristik spasial data DO.

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)")
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.6760
nama_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.2
plot(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.6761
nama_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.6
plot(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

  1. Interpolasi IDW (power = 1)

Peta IDW menunjukkan bahwa sebagian besar wilayah yang diprediksi memiliki konsentrasi DO sekitar 5,2–5,6 mg/L. Terdapat beberapa area dengan nilai lebih tinggi dan lebih rendah di sekitar titik pengamatan tertentu. Metode IDW menghasilkan permukaan prediksi berdasarkan pembobotan jarak, tetapi hasilnya sensitif terhadap nilai pengamatan dan parameter power yang digunakan.

  1. 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.

  1. 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.

  1. 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")
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
Baseline: rata-rata global 78 0.6860 0.5110 -0.026
IDW (p = 1) 78 0.6869 0.5006 -0.029 
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 0

Interpretasi hasil validasi model interpolasi

  • Model terbaik. Berdasarkan RMSE dan MAE terendah serta \(R^2\) tertinggi di antara ketiga metode interpolasi, Ordinary Kriging (Exponential) merupakan model terbaik, dengan RMSE 0,6593 mg/L. Urutan metode berdasarkan RMSE dari yang terkecil adalah Ordinary Kriging (0,6593), Universal Kriging (0,6623), baseline rata-rata global (0,6860), dan IDW dengan \(p=1\) (0,6869).

  • Keunggulan model masih kecil. Dibandingkan baseline, Ordinary Kriging menurunkan RMSE sekitar 3,9%, sedangkan Universal Kriging sekitar 3,5%. Nilai \(R^2\) Ordinary Kriging sebesar 0,052 menunjukkan kemampuan menjelaskan variasi data masih terbatas. Uji Wilcoxon menghasilkan \(p=0,1524\), sehingga keunggulan Ordinary Kriging terhadap baseline belum signifikan secara statistik.

  • Kinerja IDW. IDW menghasilkan RMSE 0,6869 mg/L, sedikit lebih tinggi daripada baseline sebesar 0,6860 mg/L, dengan \(R^2=-0,029\). Artinya, IDW dengan \(p=1\) pada pengujian ini tidak memberikan perbaikan prediksi dibandingkan penggunaan rata-rata global.

  • Perbandingan OK dan UK. Universal Kriging memiliki RMSE 0,6623 mg/L, sedikit lebih tinggi daripada Ordinary Kriging sebesar 0,6593 mg/L. Selisih 0,0030 mg/L menunjukkan bahwa penambahan tren linear berdasarkan koordinat belum memberikan peningkatan akurasi.

  • Plot observasi dan prediksi. Prediksi ketiga metode cenderung mengumpul di sekitar rata-rata, sedangkan nilai observasi berkisar antara 2,665 hingga 7,14 mg/L. Nilai ekstrem sulit diprediksi secara akurat. Sebagai contoh, DO sebesar 2,06 mg/L diprediksi jauh lebih tinggi oleh seluruh metode, sedangkan nilai 7,14 mg/L cenderung diprediksi lebih rendah. Hal ini menunjukkan bahwa model belum mampu menangkap variasi ekstrem dengan baik.

  • Penyebab akurasi terbatas. Data DO berasal dari sungai dan danau dengan karakteristik lingkungan berbeda. Variasi debit, suhu, beban pencemar, dan waktu pengambilan sampel dapat memengaruhi konsentrasi DO. Karena itu, kedekatan geografis saja belum tentu cukup untuk menjelaskan kemiripan nilai DO. Kondisi ini sejalan dengan pola variogram yang belum menunjukkan struktur spasial yang kuat.

Kesimpulan: Ordinary Kriging (Exponential) merupakan metode terbaik di antara tiga metode interpolasi yang diuji berdasarkan RMSE dan MAE. Namun, peningkatan akurasinya terhadap baseline masih kecil dan belum signifikan secara statistik. Oleh karena itu, hasil interpolasi perlu diperlakukan sebagai estimasi awal, bukan gambaran pasti kondisi DO di seluruh wilayah Aceh.

Keterbatasan

  1. Struktur spasial lemah dan fit variogram tidak stabil (peringatan tidak konvergen: OK 3 dari 3 model, UK 3 dari 3 model), sehingga parameter variogram tidak dapat ditafsirkan secara fisik.
  2. 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.
  3. 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.
  4. 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.
  5. Nilai ber-tanda “<” pada DO diperlakukan sebagai angka batas yang tertera, yang merupakan penyederhanaan.

Kesimpulan

  1. 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.

  2. 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.

  3. Berdasarkan validasi silang LOOCV, Ordinary Kriging dengan model Exponential merupakan metode interpolasi terbaik di antara tiga metode yang dibandingkan, dengan RMSE 0,6593 mg/L dan MAE 0,4783 mg/L. Universal Kriging menghasilkan RMSE 0,6623 mg/L, sedangkan IDW dengan power = 1 menghasilkan RMSE 0,6869 mg/L.

  4. Keunggulan Ordinary Kriging terhadap baseline rata-rata global masih terbatas. Penurunan RMSE hanya sekitar 3,9%, dengan nilai \(R^2=0,052\). Uji Wilcoxon berpasangan juga tidak menunjukkan perbedaan yang signifikan (\(p=0,1524\)). Hasil ini menunjukkan bahwa kemampuan model dalam memprediksi variasi DO berdasarkan lokasi masih terbatas, sehingga peta interpolasi lebih tepat digunakan untuk menggambarkan kecenderungan spasial secara umum daripada sebagai prediksi yang presisi.

  5. Penelitian selanjutnya disarankan untuk menambah titik pemantauan di wilayah Aceh bagian selatan dan timur, memeriksa kembali pencilan serta koordinat yang tidak sesuai, dan mempertimbangkan variabel tambahan seperti suhu air, elevasi, atau jarak dari muara. Penggunaan data dari beberapa tahun pengamatan juga dapat membantu mengevaluasi konsistensi pola spasial DO. Selain itu, model interpolasi perlu dievaluasi kembali setelah parameter IDW terbaik diperbarui.

Daftar Pustaka

  1. Pemerintah Republik Indonesia. (2021). Peraturan Pemerintah Nomor 22 Tahun 2021 tentang Penyelenggaraan Perlindungan dan Pengelolaan Lingkungan Hidup.

  2. Shepard, D. (1968). A two-dimensional interpolation function for irregularly-spaced data. Proceedings of the 1968 ACM National Conference, 517–524.

  3. Tobler, W. R. (1970). A computer movie simulating urban growth in the Detroit region. Economic Geography, 46(2), 234–240.

  4. Moraga, P. (2023). Spatial Statistics for Data Science: Theory and Practice with R. Chapman & Hall/CRC.

  5. 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.

  6. Pebesma, E. (2018). Simple features for R: Standardized support for spatial vector data. The R Journal, 10(1), 439–446.

  7. Pebesma, E. J. (2004). Multivariable geostatistics in S: The gstat package. Computers & Geosciences, 30(7), 683–691.