ANALISIS SPASIAL HOTSPOT KEBAKARAN HUTAN DAN LAHAN Wilayah studi: Kalimantan, Indonesia Variabel: Fire Radiative Power (FRP)Metode: IDW, Ordinary Kriging, Universal Kriging, LOOCV

1. LOAD PACKAGES

library(sf)

## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)
library(gstat)
library(sp)

2. MEMBACA DAN MEMBERSIHKAN DATA

raw <- read.delim(
  "C:\\Users\\user\\Downloads\\indonesia2023.csv",
  sep = ";",
  header = TRUE,
  stringsAsFactors = FALSE,
  colClasses = "character"
)

# Ubah koma desimal pada koordinat menjadi titik desimal
raw$latitude <- as.numeric(
  gsub(",", ".", trimws(raw$latitude), fixed = TRUE)
)

raw$longitude <- as.numeric(
  gsub(",", ".", trimws(raw$longitude), fixed = TRUE)
)

# FRP pada file asli menggunakan titik desimal
raw$frp <- as.numeric(raw$frp)
raw$confidence <- as.numeric(raw$confidence)

# Konversi tanggal
raw$acq_date <- as.Date(raw$acq_date, format = "%d/%m/%Y")

# Hapus baris dengan data utama yang tidak valid
data_bersih <- raw |>
  filter(
    is.finite(latitude),
    is.finite(longitude),
    is.finite(frp),
    !is.na(acq_date),
    latitude >= -90, latitude <= 90,
    longitude >= -180, longitude <= 180,
    frp >= 0
  )

cat("Jumlah baris raw:", nrow(raw), "\n")

## Jumlah baris raw: 52156

cat("Jumlah baris valid:", nrow(data_bersih), "\n")

## Jumlah baris valid: 52156

# Periksa hasil parsing dan statistik deskriptif
print(head(data_bersih))

##   latitude longitude   acq_date confidence  frp
## 1  -0.9296  101.7016 2023-01-01         21 16.2
## 2  -1.3048  101.8272 2023-01-01         66 27.5
## 3  -2.5064  120.3362 2023-01-02         66  6.5
## 4  -2.5153  120.3375 2023-01-02         38  4.0
## 5  -2.3944  120.8498 2023-01-02         54 10.6
## 6  -2.7547  121.7294 2023-01-02         93 50.1

print(summary(data_bersih[, c(
  "latitude", "longitude", "confidence", "frp"
)]))

##     latitude        longitude          confidence          frp         
##  Min.   :-9.999   Min.   :  0.0101   Min.   :  0.00   Min.   :   0.00  
##  1st Qu.:-6.468   1st Qu.:105.3466   1st Qu.: 50.00   1st Qu.:   9.30  
##  Median :-2.978   Median :112.7319   Median : 66.00   Median :  15.40  
##  Mean   :-3.468   Mean   :105.1973   Mean   : 63.57   Mean   :  25.94  
##  3rd Qu.:-1.305   3rd Qu.:119.8307   3rd Qu.: 79.00   3rd Qu.:  28.70  
##  Max.   : 5.571   Max.   :140.9771   Max.   :100.00   Max.   :1804.90

3. MEMILIH WILAYAH STUDI: KALIMANTAN

# Batas koordinat ini merupakan penyaringan awal,
# bukan batas administrasi resmi Pulau Kalimantan.
data_kalimantan <- data_bersih |>
  filter(
    latitude >= -4.5,
    latitude <= 4.5,
    longitude >= 108,
    longitude <= 119
  )

cat("Jumlah hotspot dalam wilayah studi:",
    nrow(data_kalimantan), "\n")

## Jumlah hotspot dalam wilayah studi: 19166

# Peta awal titik hotspot
ggplot(data_kalimantan, aes(
  x = longitude, y = latitude, color = frp
)) +
  geom_point(alpha = 0.5, size = 0.8) +
  scale_color_viridis_c(name = "FRP (MW)") +
  coord_equal() +
  labs(
    title = "Sebaran Hotspot Kebakaran di Kalimantan",
    subtitle = "Data tahun 2023",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal()

4. MENYIAPKAN DATA UNTUK INTERPOLASI

# LOOCV Kriging pada puluhan ribu titik sangat berat.
# Ambil sampel acak maksimal 200 titik untuk analisis awal.
# Hasil evaluasi hanya mewakili sampel yang digunakan.
set.seed(123)

data_sampel <- data_kalimantan |>
  distinct(latitude, longitude, .keep_all = TRUE)

if (nrow(data_sampel) < 15) {
  stop("Titik valid terlalu sedikit untuk analisis Kriging.")
}

if (nrow(data_sampel) > 200) {
  data_sampel <- data_sampel |>
    slice_sample(n = 200)
}

# Bentuk objek spasial dari koordinat geografis
titik_wgs84 <- st_as_sf(
  data_sampel,
  coords = c("longitude", "latitude"),
  crs = 4326,
  remove = FALSE
)

# Proyeksi ke koordinat planar untuk perhitungan jarak
# EPSG:3857 digunakan sebagai pendekatan awal.
titik <- st_transform(titik_wgs84, 3857)

# Simpan koordinat proyeksi sebagai kolom untuk Universal Kriging
xy <- st_coordinates(titik)
titik$X <- xy[, 1]
titik$Y <- xy[, 2]

cat("Jumlah titik yang dianalisis:", nrow(titik), "\n")

## Jumlah titik yang dianalisis: 200

5. MEMBUAT GRID PREDIKSI

bbox <- st_bbox(titik)

grid_xy <- expand.grid(
  X = seq(bbox["xmin"], bbox["xmax"], length.out = 80),
  Y = seq(bbox["ymin"], bbox["ymax"], length.out = 80)
)

grid <- st_as_sf(
  grid_xy,
  coords = c("X", "Y"),
  crs = st_crs(titik)
)

6. INTERPOLASI IDW

hasil_idw <- gstat::idw(
  frp ~ 1,
  locations = titik,
  newdata = grid,
  idp = 2
)

## Registered S3 method overwritten by 'stars':
##   method                  from
##   st_interpolate_aw.stars sf

## [inverse distance weighted interpolation]

grid$IDW <- hasil_idw$var1.pred

7. MEMBENTUK DAN MEMASANG MODEL VARIOGRAM

# Ordinary Kriging menggunakan rata-rata konstan.
vg_ok_emp <- variogram(frp ~ 1, titik)

vg_ok_fit <- fit.variogram(
  vg_ok_emp,
  model = vgm(
    psill = var(titik$frp, na.rm = TRUE) * 0.8,
    model = "Sph",
    range = max(sp::spDists(
      as(titik, "Spatial"),
      longlat = FALSE
    )) / 3,
    nugget = var(titik$frp, na.rm = TRUE) * 0.2
  )
)

## Warning in fit.variogram(vg_ok_emp, model = vgm(psill = var(titik$frp, na.rm =
## TRUE) * : No convergence after 200 iterations: try different initial values?

print(vg_ok_fit)

##   model    psill   range
## 1   Nug 2009.911       0
## 2   Sph 4779.647 2137964

plot(vg_ok_emp, vg_ok_fit, main = "Variogram Ordinary Kriging")

# Ordinary Kriging
hasil_ok <- krige(
  frp ~ 1,
  locations = titik,
  newdata = grid,
  model = vg_ok_fit
)

## [using ordinary kriging]

grid$OK <- hasil_ok$var1.pred

8. UNIVERSAL KRIGING

# Memodelkan tren linear berdasarkan koordinat X dan Y.
vg_uk_emp <- variogram(frp ~ X + Y, titik)

vg_uk_fit <- fit.variogram(
  vg_uk_emp,
  model = vgm(
    psill = var(titik$frp, na.rm = TRUE) * 0.8,
    model = "Sph",
    range = max(sp::spDists(
      as(titik, "Spatial"),
      longlat = FALSE
    )) / 3,
    nugget = var(titik$frp, na.rm = TRUE) * 0.2
  )
)

## Warning in fit.variogram(vg_uk_emp, model = vgm(psill = var(titik$frp, na.rm =
## TRUE) * : No convergence after 200 iterations: try different initial values?

print(vg_uk_fit)

##   model    psill   range
## 1   Nug 2011.779       0
## 2   Sph 4720.160 2128614

plot(vg_uk_emp, vg_uk_fit, main = "Variogram Universal Kriging")

# PERBAIKAN UNIVERSAL KRIGING

# Menambahkan koordinat X dan Y ke grid prediksi
xy_grid <- st_coordinates(grid)

grid$X <- xy_grid[, 1]
grid$Y <- xy_grid[, 2]

# Menjalankan Universal Kriging
hasil_uk <- krige(
  frp ~ X + Y,
  locations = titik,
  newdata = grid,
  model = vg_uk_fit
)

## [using universal kriging]

# Menyimpan hasil prediksi
grid$UK <- hasil_uk$var1.pred

# Melihat hasil prediksi
head(grid[, c("X", "Y", "UK")])

## Simple feature collection with 6 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 12024640 ymin: -438018.7 xmax: 12093820 ymax: -438018.7
## Projected CRS: WGS 84 / Pseudo-Mercator
##          X         Y       UK                   geometry
## 1 12024642 -438018.7 21.59676 POINT (12024642 -438018.7)
## 2 12038477 -438018.7 21.94573 POINT (12038477 -438018.7)
## 3 12052312 -438018.7 22.31243 POINT (12052312 -438018.7)
## 4 12066147 -438018.7 22.69483 POINT (12066147 -438018.7)
## 5 12079981 -438018.7 23.09029 POINT (12079981 -438018.7)
## 6 12093816 -438018.7 23.49583 POINT (12093816 -438018.7)

# Membuat peta Universal Kriging
koordinat_grid <- st_coordinates(grid)

data_peta_uk <- data.frame(
  X = koordinat_grid[, 1],
  Y = koordinat_grid[, 2],
  UK = grid$UK
)

ggplot(data_peta_uk, aes(x = X, y = Y, fill = UK)) +
  geom_raster() +
  scale_fill_viridis_c(name = "FRP (MW)") +
  coord_equal() +
  labs(
    title = "Peta Interpolasi Universal Kriging",
    subtitle = "Sebaran FRP Hotspot Kalimantan Tahun 2023",
    x = "Koordinat X",
    y = "Koordinat Y"
  ) +
  theme_minimal()

10. LOOCV IDW

n <- nrow(titik)
aktual <- titik$frp
pred_idw <- rep(NA_real_, n)

for (i in seq_len(n)) {
  latih <- titik[-i, ]
  uji <- titik[i, ]

  hasil <- gstat::idw(
    frp ~ 1,
    locations = latih,
    newdata = uji,
    idp = 2
  )

  pred_idw[i] <- hasil$var1.pred
}

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]

11. LOOCV ORDINARY KRIGING

cv_ok <- krige.cv(
  frp ~ 1,
  locations = titik,
  model = vg_ok_fit,
  nfold = n
)

pred_ok <- cv_ok$var1.pred

12. LOOCV UNIVERSAL KRIGING

cv_uk <- krige.cv(
  frp ~ X + Y,
  locations = titik,
  model = vg_uk_fit,
  nfold = n
)

pred_uk <- cv_uk$var1.pred

13. EVALUASI RMSE DAN MAE

metrik <- function(aktual, prediksi) {
  valid <- is.finite(aktual) & is.finite(prediksi)

  c(
    RMSE = sqrt(mean((aktual[valid] - prediksi[valid])^2)),
    MAE = mean(abs(aktual[valid] - prediksi[valid])),
    N = sum(valid)
  )
}

evaluasi <- rbind(
  IDW = metrik(aktual, pred_idw),
  Ordinary_Kriging = metrik(aktual, pred_ok),
  Universal_Kriging = metrik(aktual, pred_uk)
)

print(round(evaluasi, 4))

##                      RMSE     MAE   N
## IDW               57.3387 22.3901 200
## Ordinary_Kriging  57.1862 21.7909 200
## Universal_Kriging 57.2115 21.7911 200

# Metode terbaik berdasarkan RMSE
rmse <- evaluasi[, "RMSE"]
if (all(is.finite(rmse))) {
  cat(
    "\nMetode dengan RMSE terkecil:",
    rownames(evaluasi)[which.min(rmse)],
    "\n"
  )
}

## 
## Metode dengan RMSE terkecil: Ordinary_Kriging

UTS LINGKUNGAN

Tazkia Aulia - 3338240001

1. LOAD PACKAGES

2. MEMBACA DAN MEMBERSIHKAN DATA

3. MEMILIH WILAYAH STUDI: KALIMANTAN

4. MENYIAPKAN DATA UNTUK INTERPOLASI

5. MEMBUAT GRID PREDIKSI

6. INTERPOLASI IDW

7. MEMBENTUK DAN MEMASANG MODEL VARIOGRAM

8. UNIVERSAL KRIGING

10. LOOCV IDW

11. LOOCV ORDINARY KRIGING

12. LOOCV UNIVERSAL KRIGING

13. EVALUASI RMSE DAN MAE