Eksplorasi Data Lingkungan: Kebakaran Hutan dan Lahan (Karhutla) serta Curah Hujan di Indonesia
Mata Kuliah Statistika Lingkungan (STL623301) — Statistika UNTIRTA
Tazkia Aulia (3338240001)
12 September 2026
1 Pendahuluan
1.1 Tentang Studi Kasus
Project ini menerapkan teknik eksplorasi data lingkungan yang dibahas pada Pertemuan 2 statistik deskriptif, visualisasi, deteksi outlier, dan autokorelasi (temporal & spasial), pada studi kasus luas kebakaran hutan dan lahan (karhutla) dan curah hujan di 34 provinsi Indonesia.
Dataset utama (Data_Project) berisi 190
baris (34 dan atau 38 provinsi × 5 tahun, 2022–2026), memuat
rata-rata curah hujan tahunan (rfh, mm), luas karhutla
tahunan (luas_karhutla_ha), serta koordinat
(lat, long) tiap provinsi. Data disusun dari
dua sumber mentah pada file yang sama: sheet Data Karhutla
(rekap luas karhutla per provinsi per tahun) dan sheet
Curah Hujan 2022-2026 (data curah hujan dekadal/10-harian
per wilayah administrasi yang digabungkan dengan koordinat pada sheet
Bujur lintang.
1.2 Rumusan Pertanyaan Analisis
- Bagaimana karakteristik statistik deskriptif luas karhutla dan curah hujan secara nasional maupun per provinsi?
- Bagaimana pola distribusi, tren tahunan, dan sebaran spasial karhutla dapat digambarkan melalui visualisasi?
- Apakah terdapat provinsi/tahun dengan luas karhutla yang tergolong outlier?
- Apakah terdapat indikasi autokorelasi temporal (antarwaktu) pada curah hujan dan autokorelasi spasial (antarprovinsi berdekatan) pada karhutla?
Keempat pertanyaan ini dijawab bertahap melalui delapan langkah analisis di bawah ini.
2 Persiapan
2.1 Memuat Package
library(readxl) # membaca file Excel
library(dplyr) # manipulasi data
library(tidyr) # merapikan data
library(ggplot2) # visualisasi
library(sf) # objek spasial
library(tmap) # peta tematik
library(spdep) # autokorelasi spasial (Moran's I)
library(rnaturalearth) # garis batas negara Indonesia (basemap, tanpa perlu shapefile manual)2.2 Langkah 1 — Import dan Struktur Data
# Input data
data_karhutla <- read_excel("C:/Users/user/OneDrive - untirta.ac.id/Datacoba1 (1).xlsx", sheet = "Data_Project")
str(data_karhutla)## tibble [190 × 8] (S3: tbl_df/tbl/data.frame)
## $ kode_provinsi : chr [1:190] "11" "11" "11" "11" ...
## $ provinsi : chr [1:190] "Aceh" "Aceh" "Aceh" "Aceh" ...
## $ tahun : num [1:190] 2022 2023 2024 2025 2026 ...
## $ rfh : chr [1:190] "87.065249777777808" "76.088721000000007" "72.041313486111122" "84.924350111111096" ...
## $ rfq : chr [1:190] "119.20335880555558" "101.73461622222224" "98.849901000000003" "110.10899866666668" ...
## $ luas_karhutla_ha: num [1:190] 3716 1937 7257 8925 4070 ...
## $ lat : chr [1:190] "4.3685499999999999" "4.3685499999999999" "4.3685499999999999" "4.3685499999999999" ...
## $ long : chr [1:190] "97.025300000000001" "97.025300000000001" "97.025300000000001" "97.025300000000001" ...data_karhutla <- data_karhutla %>%
dplyr::mutate(dplyr::across(c(rfh, rfq, lat, long), as.numeric))
str(data_karhutla)## tibble [190 × 8] (S3: tbl_df/tbl/data.frame)
## $ kode_provinsi : chr [1:190] "11" "11" "11" "11" ...
## $ provinsi : chr [1:190] "Aceh" "Aceh" "Aceh" "Aceh" ...
## $ tahun : num [1:190] 2022 2023 2024 2025 2026 ...
## $ rfh : num [1:190] 87.1 76.1 72 84.9 67.1 ...
## $ rfq : num [1:190] 119.2 101.7 98.8 110.1 103.9 ...
## $ luas_karhutla_ha: num [1:190] 3716 1937 7257 8925 4070 ...
## $ lat : num [1:190] 4.37 4.37 4.37 4.37 4.37 ...
## $ long : num [1:190] 97 97 97 97 97 ...head(data_karhutla, 190)## # A tibble: 190 × 8
## kode_provinsi provinsi tahun rfh rfq luas_karhutla_ha lat long
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 11 Aceh 2022 87.1 119. 3716 4.37 97.0
## 2 11 Aceh 2023 76.1 102. 1937. 4.37 97.0
## 3 11 Aceh 2024 72.0 98.8 7257. 4.37 97.0
## 4 11 Aceh 2025 84.9 110. 8925. 4.37 97.0
## 5 11 Aceh 2026 67.1 104. 4070. 4.37 97.0
## 6 51 Bali 2022 76.3 153. 0 -8.24 115.
## 7 51 Bali 2023 46.3 81.7 2011. -8.24 115.
## 8 51 Bali 2024 59.2 98.6 316. -8.24 115.
## 9 51 Bali 2025 72.0 137. 93.1 -8.24 115.
## 10 51 Bali 2026 64.6 102. 108. -8.24 115.
## # ℹ 180 more rowsdim(data_karhutla)## [1] 190 8# Memeriksa data hilang per kolom
colSums(is.na(data_karhutla))## kode_provinsi provinsi tahun rfh
## 0 0 0 20
## rfq luas_karhutla_ha lat long
## 20 0 20 20data_karhutla %>%
dplyr::filter(is.na(lat)) %>%
dplyr::distinct(provinsi)## # A tibble: 4 × 1
## provinsi
## <chr>
## 1 Papua Barat
## 2 Papua Barat Daya
## 3 Papua Pegunungan
## 4 Papua Selatan2.2.1 Interpretasi
Dataset berisi 190 baris dan 8 kolom, sesuai rancangan 38 provinsi ×
5 tahun (2022–2026). Kolom luas_karhutla_ha dan
rfh sudah bertipe numerik dan provinsi bertipe
karakter.
Dataset tidak sepenuhnya lengkap: kolom rfh,
rfq, lat, dan long masing-masing
memiliki 20 nilai hilang (20 = 4 provinsi × 5 tahun), seluruhnya berasal
dari empat provinsi hasil pemekaran Papua 2022 (Papua Barat, Papua Barat
Daya, Papua Pegunungan, Papua Selatan) yang belum tersedia data curah
hujan maupun koordinatnya pada sumber data ini. Kemungkinan karena belum
memperbarui pembagian wilayah administratifnya. Sebaliknya,
luas_karhutla_ha tetap lengkap untuk keempat provinsi
tersebut karena berasal dari sumber data yang berbeda (sheet Data
Karhutla). Ketidaklengkapan ini ditangani dengan mengecualikan baris
terkait hanya dari analisis yang membutuhkan rfh atau
koordinat, bukan menghapusnya dari keseluruhan dataset
3 Statistik Deskriptif
3.1 Langkah 2 — Statistik Deskriptif
# Statistik deskriptif keseluruhan — luas karhutla
summary(data_karhutla$luas_karhutla_ha)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 384.4 2404.6 12129.1 8288.4 190394.6sd(data_karhutla$luas_karhutla_ha)## [1] 28932.59IQR(data_karhutla$luas_karhutla_ha)## [1] 7903.957# Statistik deskriptif keseluruhan — curah hujan (rfh)
summary(data_karhutla$rfh)## Min. 1st Qu. Median Mean 3rd Qu. Max. NAs
## 42.21 64.20 78.28 78.92 94.23 122.16 20sd(data_karhutla$rfh, na.rm = TRUE)## [1] 18.59752Rata-rata luas karhutla nasional adalah 12.129 ha per provinsi per tahun, jauh di atas mediannya (2.406 ha), dengan simpangan baku yang sangat besar (28.932,59 ha) hampir 2,4 kali nilai rata-ratanya sendiri. Selisih mean–median yang lebar ini adalah indikasi kuat distribusi right-skewed: sebagian kecil provinsi/tahun dengan karhutla ekstrem menarik rata-rata jauh di atas kondisi yang tipikal, sehingga median lebih representatif sebagai ukuran pemusatan untuk data ini.
Rata-rata curah hujan sebesar 78,92, sedangkan mediannya 78,28. Nilai yang hampir sama menunjukkan bahwa data curah hujan relatif simetris. Simpangan baku sebesar 18,60 menunjukkan adanya variasi curah hujan yang cukup terlihat dari rata-ratanya. Sebagian besar data berada pada rentang 64,20–94,23, dengan 20 data yang memiliki nilai NA.
# Statistik deskriptif per provinsi (diurutkan dari rerata karhutla tertinggi)
rekap_provinsi <- data_karhutla %>%
group_by(provinsi) %>%
summarise(
n = n(),
rerata = mean(luas_karhutla_ha),
median = median(luas_karhutla_ha),
sd = sd(luas_karhutla_ha),
min = min(luas_karhutla_ha),
maks = max(luas_karhutla_ha)
) %>%
arrange(desc(rerata))
rekap_provinsi## # A tibble: 38 × 7
## provinsi n rerata median sd min maks
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Nusa Tenggara Timur 5 84848. 102537. 40929. 19565. 118409.
## 2 Kalimantan Barat 5 44571. 26703. 38139. 21836 111848.
## 3 Kalimantan Selatan 5 43147. 8020. 82421. 429 190395.
## 4 Papua Selatan 5 42388. 25839. 61711. 0 150813.
## 5 Kalimantan Tengah 5 38270. 7634. 71439. 1554 165896.
## 6 Nusa Tenggara Barat 5 37743. 30567 19197. 17406. 66716.
## 7 Sumatera Selatan 5 31819. 5940. 56290. 1926. 132083.
## 8 Maluku 5 22194. 14954 13720. 13405. 45999.
## 9 Jawa Timur 5 18868. 10835. 18345. 2380 49498.
## 10 Kalimantan Timur 5 14945. 9573. 16214. 373 39494.
## # ℹ 28 more rowsRekap per provinsi menunjukkan adanya variasi spasial yang cukup jelas dalam luas kebakaran hutan dan lahan (karhutla).
Nusa Tenggara Timur memiliki rerata karhutla tertinggi (84.848 ha), diikuti Kalimantan Barat (44.571 ha), Kalimantan Selatan (43.147 ha), Papua Selatan (42.388 ha), dan Kalimantan Tengah (38.270 ha). Beberapa provinsi tersebut juga memiliki nilai maksimum yang cukup tinggi, seperti Kalimantan Selatan yang mencapai 190.395 ha.
Sementara itu, provinsi lainnya memiliki rerata karhutla yang lebih rendah, dengan Jawa Timur sebesar 18.868 ha dan Kalimantan Timur sebesar 14.945 ha. Perbedaan antara rerata dan median pada beberapa provinsi, terutama Kalimantan Selatan, Kalimantan Tengah, dan Sumatera Selatan, juga menunjukkan adanya variasi luas karhutla yang cukup besar antar tahun.
# Rerata per tahun — melihat tahun mana yang paling parah secara nasional
rekap_tahun <- data_karhutla %>%
group_by(tahun) %>%
summarise(
rerata_karhutla = mean(luas_karhutla_ha),
rerata_hujan = mean(rfh, na.rm = TRUE)
)
rekap_tahun## # A tibble: 5 × 3
## tahun rerata_karhutla rerata_hujan
## <dbl> <dbl> <dbl>
## 1 2022 5392. 90.8
## 2 2023 30558. 66.5
## 3 2024 9916. 79.6
## 4 2025 9464. 92.7
## 5 2026 5316. 65.03.1.1 Interpretasi
Rekap per tahun menunjukkan 2023 sebagai tahun terparah secara nasional (rerata karhutla 30.558 ha, lebih dari 3 kali tahun lainnya), justru bersamaan dengan curah hujan rata-rata terendah kedua (66,5 mm).
4 Visualisasi Data
4.1 Langkah 3 — Histogram dan Boxplot
# Histogram distribusi luas karhutla keseluruhan
ggplot(data_karhutla, aes(x = luas_karhutla_ha)) +
geom_histogram(bins = 25, fill = "#147D64", color = "white") +
labs(title = "Distribusi Luas Karhutla — Indonesia 2022–2026",
x = "Luas Karhutla (ha)", y = "Frekuensi") +
theme_minimal()
4.1.1 Interpretasi
Berdasarkan histogram, distribusi luas karhutla terlihat sangat menceng ke kanan (positively skewed). Sebagian besar observasi berada pada luas di bawah 20.000 ha, sedangkan hanya sedikit observasi yang mencapai lebih dari 50.000 ha, bahkan ada yang mendekati 190.000 ha. Hal ini menunjukkan bahwa sebagian besar kejadian karhutla relatif kecil, tetapi terdapat beberapa kejadian dengan luas yang sangat besar sehingga membentuk ekor panjang di sisi kanan distribusi.
# Boxplot luas karhutla per provinsi (diurutkan berdasarkan median)
ggplot(data_karhutla, aes(x = reorder(provinsi, luas_karhutla_ha, median),
y = luas_karhutla_ha)) +
geom_boxplot(fill = "#FFC200") +
coord_flip() +
labs(title = "Sebaran Luas Karhutla per Provinsi (2022–2026)",
x = "Provinsi", y = "Luas Karhutla (ha)") +
theme_minimal()
4.1.2 Interpretasi
Boxplot per provinsi (diurutkan berdasarkan median) menunjukkan:
Nusa Tenggara Timur memiliki median luas karhutla paling tinggi, disusul Nusa Tenggara Barat, Kalimantan Barat, dan Papua Selatan.
Kalimantan Selatan, Kalimantan Tengah, dan Sumatera Selatan memiliki kotak dan whisker yang relatif lebar, menunjukkan variasi luas karhutla yang cukup besar antartahun. Titik yang terpisah jauh di atas whisker menunjukkan adanya nilai ekstrem.
Sebaliknya, sebagian besar provinsi dengan luas karhutla lebih rendah memiliki kotak yang lebih sempit, menandakan sebaran karhutla yang lebih kecil dan relatif terkonsentrasi.
4.2 Langkah 4 — Diagram Deret Waktu
# Fokus pada 8 provinsi dengan rerata karhutla tertinggi agar grafik tetap terbaca
provinsi_top8 <- rekap_provinsi$provinsi[1:8]
data_karhutla %>%
dplyr::filter(provinsi %in% provinsi_top8) %>%
ggplot(aes(x = tahun, y = luas_karhutla_ha, color = provinsi)) +
geom_line(linewidth = 1) +
geom_point(size = 2) +
scale_x_continuous(breaks = 2022:2026) +
labs(title = "Tren Tahunan Luas Karhutla — 8 Provinsi Tertinggi",
x = "Tahun", y = "Luas Karhutla (ha)", color = "Provinsi") +
theme_minimal()
4.2.1 Interpretasi
Grafik deret waktu pada 8 provinsi ber-karhutla tertinggi memperlihatkan lonjakan tajam dan serentak pada 2023 di hampir semua provinsi tersebut paling ekstrem di Kalimantan Selatan, Kalimantan Tengah, dan Papua Selatan — sebelum kembali menurun tajam pada 2024–2026 mendekati level 2022.
# Scatter plot: rerata curah hujan vs rerata karhutla per tahun (agregat nasional)
ggplot(rekap_tahun, aes(x = rerata_hujan, y = rerata_karhutla)) +
geom_point(color = "blue", size = 3) +
geom_text(aes(label = tahun), vjust = -1) +
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(title = "Rerata Curah Hujan vs Rerata Luas Karhutla per Tahun (Nasional)",
x = "Rerata Curah Hujan (mm)", y = "Rerata Luas Karhutla (ha)") +
theme_minimal()
cor(rekap_tahun$rerata_hujan, rekap_tahun$rerata_karhutla)## [1] -0.46042314.2.2 Interpretasi Scatter Plot
Scatter plot menunjukkan kecenderungan hubungan negatif antara rerata curah hujan dan rerata luas karhutla pada tingkat nasional selama periode 2022–2026. Hal ini terlihat dari garis regresi yang memiliki arah menurun, yang mengindikasikan bahwa tahun dengan rerata curah hujan lebih tinggi cenderung memiliki rerata luas karhutla yang lebih rendah.
Variasi antar tahun terlihat cukup besar, terutama pada tahun 2023 yang memiliki rerata luas karhutla paling tinggi meskipun rerata curah hujannya relatif rendah. Sementara itu, tahun 2026 menunjukkan curah hujan paling rendah, tetapi luas karhutlanya tidak setinggi tahun 2023. Hal ini menunjukkan bahwa curah hujan bukan satu-satunya faktor yang berkaitan dengan luas karhutla
Pita kepercayaan yang lebar dan jumlah data yang hanya lima titik menunjukkan bahwa hubungan tersebut masih bersifat indikatif, sehingga diperlukan data yang lebih banyak untuk memperoleh hasil yang lebih kuat.
4.2.3 Interpretasi Korelasi
Nilai korelasi r = -0,46 menunjukkan hubungan negatif sedang antara curah hujan dan luas karhutla. Artinya, peningkatan curah hujan cenderung diikuti penurunan luas karhutla.
4.3 Langkah 5 Peta Tematik
# Garis batas Indonesia diambil otomatis (tidak perlu shapefile provinsi manual)
indo <- ne_countries(country = "Indonesia", scale = "medium", returnclass = "sf")
# Titik rerata karhutla per provinsi (2022-2026), provinsi tanpa koordinat otomatis dikecualikan
titik_provinsi <- data_karhutla %>%
dplyr::filter(!is.na(lat), !is.na(long)) %>%
tidyr::drop_na(lat, long) %>% # jaring pengaman kedua
dplyr::group_by(provinsi, lat, long) %>%
dplyr::summarise(rerata_karhutla = mean(luas_karhutla_ha), .groups = "drop")
# Sanity check: pastikan tidak ada NA tersisa sebelum dibuat objek spasial
stopifnot(sum(is.na(titik_provinsi$lat)) == 0, sum(is.na(titik_provinsi$long)) == 0)
titik_provinsi <- titik_provinsi %>%
st_as_sf(coords = c("long", "lat"), crs = 4326, remove = FALSE)
tm_shape(indo) +
tm_borders(col = "grey50") +
tm_shape(titik_provinsi) +
tm_symbols(size = "rerata_karhutla", col = "rerata_karhutla",
palette = "YlOrRd", scale = 2,
title.col = "Rerata Karhutla (ha)",
title.size = "Rerata Karhutla (ha)") +
tm_layout(title = "Sebaran Rerata Luas Karhutla per Provinsi (2022-2026)",
legend.outside = TRUE)
4.4 Langkah 5 Peta Tematik
## Menggunakan shapefile provinsi asli (indonesia.shp/.shx/.dbf/.prj),
peta_provinsi <- sf::st_read("C:/Users/user/OneDrive - untirta.ac.id/Documents/Semester 5/shp indo", quiet = TRUE)
# Cek nama provinsi di shapefile vs di data
sort(peta_provinsi$name)## [1] "Aceh" "Bali" "Bangka-Belitung"
## [4] "Banten" "Bengkulu" "Gorontalo"
## [7] "Jakarta Raya" "Jambi" "Jawa Barat"
## [10] "Jawa Tengah" "Jawa Timur" "Kalimantan Barat"
## [13] "Kalimantan Selatan" "Kalimantan Tengah" "Kalimantan Timur"
## [16] "Kepulauan Riau" "Lampung" "Maluku"
## [19] "Maluku Utara" "North Kalimantan" "Nusa Tenggara Barat"
## [22] "Nusa Tenggara Timur" "Papua" "Papua Barat"
## [25] "Riau" "Sulawesi Barat" "Sulawesi Selatan"
## [28] "Sulawesi Tengah" "Sulawesi Tenggara" "Sulawesi Utara"
## [31] "Sumatera Barat" "Sumatera Selatan" "Sumatera Utara"
## [34] "Yogyakarta"# Penyesuaian nama provinsi
# Shapefile ini masih memakai batas 34 provinsi (SEBELUM pemekaran Papua 2022),
# jadi 5 provinsi baru di data kita (Papua, Papua Tengah, Papua Pegunungan,
# Papua Selatan, Papua Barat Daya) tidak punya poligon sendiri.
# Solusi: agregasi data provinsi baru ke wilayah induknya agar tetap bisa dipetakan.
crosswalk_papua <- c(
"Papua" = "Papua",
"Papua Tengah" = "Papua",
"Papua Pegunungan" = "Papua",
"Papua Selatan" = "Papua",
"Papua Barat" = "Papua Barat",
"Papua Barat Daya" = "Papua Barat"
)
titik_provinsi_peta <- data_karhutla %>%
dplyr::mutate(
provinsi_peta = dplyr::case_when(
provinsi %in% names(crosswalk_papua) ~ crosswalk_papua[provinsi],
provinsi == "DKI Jakarta" ~ "Jakarta Raya", # beda penulisan di shapefile
provinsi == "Bangka Belitung" ~ "Bangka-Belitung", # beda penulisan di shapefile
provinsi == "Kalimantan Utara" ~ "North Kalimantan", # beda penulisan di shapefile
TRUE ~ provinsi
)
) %>%
dplyr::group_by(provinsi_peta) %>%
dplyr::summarise(rerata_karhutla = mean(luas_karhutla_ha), .groups = "drop")
# Gabungkan data dengan poligon provinsi
peta_karhutla <- peta_provinsi %>%
dplyr::left_join(titik_provinsi_peta, by = c("name" = "provinsi_peta"))
# Peta choropleth: seluruh wilayah provinsi diwarnai
tm_shape(peta_karhutla) +
tm_polygons(col = "rerata_karhutla",
palette = "YlOrRd", style = "jenks",
title = "Rerata Karhutla (ha)",
border.col = "grey40", colorNA = "grey90",
textNA = "Data tidak tersedia") +
tm_layout(title = "Sebaran Rerata Luas Karhutla per Provinsi (2022-2026)",
legend.outside = TRUE)
4.4.1 Interpretasi
Peta tematik memperlihatkan pola spasial yang cukup jelas: simbol besar dan berwarna gelap terkonsentrasi di Nusa Tenggara Timur dan Kalimantan sedangkan Pulau Jawa hampir seluruhnya bersimbol kecil dan terang. Pola “provinsi berdekatan cenderung mirip” ini menjadi petunjuk visual awal adanya autokorelasi spasial positif yang akan diuji formal pada Langkah 8.
5 Deteksi Outlier
5.1 Langkah 6 — Deteksi Outlier
# Metode IQR
Q1 <- quantile(data_karhutla$luas_karhutla_ha, 0.25)
Q3 <- quantile(data_karhutla$luas_karhutla_ha, 0.75)
IQR_val <- Q3 - Q1
batas_bawah <- Q1 - 1.5 * IQR_val
batas_atas <- Q3 + 1.5 * IQR_val
cat("Batas bawah:", round(batas_bawah, 1), " | Batas atas:", round(batas_atas, 1), "\n")## Batas bawah: -11471.5 | Batas atas: 20144.3outlier_iqr <- data_karhutla %>%
dplyr::filter(luas_karhutla_ha < batas_bawah | luas_karhutla_ha > batas_atas) %>%
dplyr::arrange(desc(luas_karhutla_ha))
outlier_iqr %>% select(provinsi, tahun, luas_karhutla_ha)## # A tibble: 27 × 3
## provinsi tahun luas_karhutla_ha
## <chr> <dbl> <dbl>
## 1 Kalimantan Selatan 2023 190395.
## 2 Kalimantan Tengah 2023 165896.
## 3 Papua Selatan 2023 150813.
## 4 Sumatera Selatan 2023 132083.
## 5 Nusa Tenggara Timur 2025 118409.
## 6 Nusa Tenggara Timur 2024 113090.
## 7 Kalimantan Barat 2023 111848.
## 8 Nusa Tenggara Timur 2023 102537.
## 9 Nusa Tenggara Timur 2022 70637
## 10 Nusa Tenggara Barat 2023 66716.
## # ℹ 17 more rows# Metode Z-score
data_karhutla$z_karhutla <- as.numeric(scale(data_karhutla$luas_karhutla_ha))
outlier_z <- data_karhutla %>%
dplyr::filter(abs(z_karhutla) > 3) %>%
dplyr::arrange(desc(z_karhutla))
outlier_z %>% select(provinsi, tahun, luas_karhutla_ha, z_karhutla)## # A tibble: 8 × 4
## provinsi tahun luas_karhutla_ha z_karhutla
## <chr> <dbl> <dbl> <dbl>
## 1 Kalimantan Selatan 2023 190395. 6.16
## 2 Kalimantan Tengah 2023 165896. 5.31
## 3 Papua Selatan 2023 150813. 4.79
## 4 Sumatera Selatan 2023 132083. 4.15
## 5 Nusa Tenggara Timur 2025 118409. 3.67
## 6 Nusa Tenggara Timur 2024 113090. 3.49
## 7 Kalimantan Barat 2023 111848. 3.45
## 8 Nusa Tenggara Timur 2023 102537. 3.125.1.1 Interpretasi
Dengan aturan IQR, batas atas kewajaran adalah 20.144,3 ha; sebanyak 27 dari 190 observasi (≈14,2%) melampauinya. Proporsi ini terlihat besar untuk standar IQR konvensional, namun wajar mengingat data memang sangat right-skewed pada distribusi seperti ini, aturan 1,5×IQR cenderung menandai banyak titik karena IQR-nya sendiri relatif sempit dibanding rentang datanya. Metode Z-score yang lebih konservatif (|z| > 3) menandai hanya 8 observasi paling ekstrem.
Kedua metode sepakat bahwa observasi paling ekstrem berasal dari tahun 2023, tersebar di beberapa provinsi rawan karhutla (Kalimantan Selatan, Kalimantan Tengah, Papua Selatan, Sumatera Selatan, Nusa Tenggara Timur, Kalimantan Barat).
6 Autokorelasi
6.1 Langkah 7 - Autokorelasi Temporal (ACF)
Deret waktu tahunan karhutla per provinsi hanya
memiliki 5 titik (2022–2026) terlalu pendek untuk ACF yang andal.
Sebagai gantinya, autokorelasi temporal didemonstrasikan pada data
curah hujan dekadal (10-harian) dari sheet
Curah Hujan 2022-2026, yang tersedia untuk 34 provinsi
dengan ±168 titik waktu per provinsi (Januari 2022 – Agustus 2026)
dimana lebih memadai untuk melihat pola autokorelasi/musiman secara
meyakinkan.
# Curah hujan dekadal level provinsi (adm_level == 1), contoh: Nusa Tenggara Timur (PCODE "ID53")
# —provinsi dengan rerata karhutla tertinggi pada rekap_provinsi di Langkah 2
data_hujan_dekadal <- read_excel("C:/Users/user/OneDrive - untirta.ac.id/Datacoba1 (1).xlsx", sheet = "Curah Hujan 2022-2026")
ntt_dekadal <- data_hujan_dekadal %>%
dplyr::filter(adm_level == 1, PCODE == "ID53") %>%
dplyr::arrange(date)
ts_ntt <- ts(ntt_dekadal$rfh, start = c(2022, 1), frequency = 36) # 36 dekade per tahun
acf(ts_ntt, lag.max = 40, main = "ACF Curah Hujan Dekadal — Nusa Tenggara Timur")
6.1.1 Interpretasi
Grafik ACF menunjukkan bahwa pola curah hujan di Nusa Tenggara Timur sangat dipengaruhi oleh musim. Nilai autokorelasi yang tinggi pada lag awal menunjukkan bahwa jika suatu periode memiliki curah hujan tinggi, maka periode-periode setelahnya juga cenderung memiliki curah hujan tinggi. Sebaliknya, pada lag menengah autokorelasi menjadi negatif, yang mencerminkan perbedaan antara musim hujan dan musim kemarau. Setelah mendekati satu tahun, autokorelasi kembali positif, menandakan bahwa pola curah hujan tahun ini cenderung berulang pada periode yang sama di tahun berikutnya. Dengan kata lain, curah hujan di NTT memiliki siklus musiman tahunan yang kuat dan relatif konsisten dari tahun ke tahun.
6.2 Langkah 8 - Autokorelasi Spasial (Moran’s I)
# Matriks pembobot spasial: k = 4 tetangga terdekat, dari rerata karhutla per provinsi (Langkah 5)
# (provinsi tanpa koordinat otomatis dikecualikan — lihat catatan pada Langkah 5)
koordinat <- cbind(titik_provinsi$long, titik_provinsi$lat)
knn4 <- knearneigh(koordinat, k = 4)
wij <- nb2listw(knn2nb(knn4), style = "W")
moran.test(titik_provinsi$rerata_karhutla, listw = wij)##
## Moran I test under randomisation
##
## data: titik_provinsi$rerata_karhutla
## weights: wij
##
## Moran I statistic standard deviate = 1.8432, p-value = 0.03265
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.140049158 -0.030303030 0.0085421346.2.1 Interpretasi
Hasil uji Moran’s I menghasilkan nilai sebesar 0,140 dengan p-value 0,0327. Karena p-value lebih kecil dari 0,05, terdapat autokorelasi spasial positif yang signifikan pada data rerata luas karhutla antarprovinsi di Indonesia.
Autokorelasi spasial positif berarti provinsi yang berdekatan secara geografis cenderung memiliki tingkat karhutla yang mirip. Dengan kata lain, wilayah dengan karhutla tinggi cenderung berada di dekat wilayah lain yang juga memiliki karhutla tinggi, sedangkan wilayah dengan karhutla rendah cenderung berdekatan dengan wilayah yang juga rendah.
Nilai Moran’s I sebesar 0,140 menunjukkan bahwa hubungan spasial tersebut ada, tetapi tidak terlalu kuat. Hal ini mengindikasikan bahwa selain faktor lokasi geografis, masih terdapat faktor-faktor lain yang turut memengaruhi variasi karhutla antarprovinsi, seperti kondisi iklim, tutupan lahan, dan aktivitas manusia. Temuan kuantitatif ini mengonfirmasi pola yang sudah terlihat secara visual pada peta tematik Langkah 5, dan sejalan dengan
Hukum Pertama Tobler
7 Kesimpulan
1. Karakteristik statistik deskriptif. Secara nasional, luas kebakaran hutan dan lahan (karhutla) per provinsi per tahun memiliki rata-rata 12.129 ha, sedangkan mediannya hanya 2.405 ha (sd = 28.933 ha). Perbedaan yang besar antara rata-rata dan median menunjukkan bahwa distribusi data sangat menceng ke kanan (right-skewed), dengan sebagian kecil observasi bernilai sangat tinggi yang menarik rata-rata ke atas. Secara spasial, rerata karhutla tertinggi ditemukan di beberapa provinsi Kalimantan, Sumatera Selatan, Nusa Tenggara Timur, dan Papua Selatan, sedangkan provinsi di Pulau Jawa umumnya memiliki nilai yang lebih rendah.
2. Pola distribusi, tren, dan sebaran spasial. Histogram dan boxplot menunjukkan bahwa sebagian besar observasi memiliki luas karhutla relatif rendah, sementara hanya sedikit observasi yang memiliki nilai sangat besar. Analisis tren tahunan menunjukkan bahwa tahun 2023 merupakan tahun dengan rerata karhutla tertinggi secara nasional, lebih dari tiga kali lipat dibandingkan tahun lainnya. Selain itu, terdapat indikasi hubungan negatif antara curah hujan dan luas karhutla, di mana tahun dengan curah hujan lebih rendah cenderung memiliki luas karhutla yang lebih besar. Peta tematik juga memperlihatkan adanya kecenderungan pengelompokan wilayah dengan tingkat karhutla yang serupa.
3. Provinsi/tahun dengan karhutla outlier. Deteksi outlier menggunakan metode IQR dan Z-score menunjukkan adanya sejumlah observasi ekstrem, terutama pada tahun 2023. Kedua metode secara konsisten mengidentifikasi provinsi-provinsi seperti Kalimantan Selatan, Kalimantan Tengah, Papua Selatan, Sumatera Selatan, Nusa Tenggara Timur, dan Kalimantan Barat sebagai lokasi dengan kejadian karhutla yang sangat tinggi dibandingkan pola umum data.
4. Autokorelasi temporal dan spasial Analisis ACF pada curah hujan dekadal di Nusa Tenggara Timur menunjukkan adanya pola musiman tahunan yang kuat dan konsisten, ditandai oleh pengulangan pola curah hujan dari tahun ke tahun. Sementara itu, hasil uji Moran’s I (I = 0,140; p = 0,033) menunjukkan adanya autokorelasi spasial positif yang signifikan pada rerata karhutla antarprovinsi. Artinya, provinsi yang berdekatan secara geografis cenderung memiliki tingkat karhutla yang mirip. Meskipun demikian, kekuatan autokorelasi spasial tersebut relatif lemah hingga sedang, sehingga faktor-faktor lain seperti iklim, tutupan lahan, dan aktivitas manusia juga berperan dalam membentuk variasi karhutla di Indonesia.
8 Referensi
https://data.humdata.org/m/dataset/idn-rainfall-subnational (Data curah hujan Indonesia untuk 5 tahun terakhir)
https://sipongi.gakkum.kehutanan.go.id/indikasi-luas-kebakaran (Data Luas Kebakaran Hutan dan Lahan (Ha) per Provinsi di Indonesia)
https://github.com/yusufsyaifudin/wilayah-indonesia/blob/master/data/list_of_area/provinces.json (List Lintang Bujur)
https://www.instagram.com/p/Dc01DbFsEjm/?stkn=eG82cm1nOGV0c2Nt (Inspirasi)