Pendahuluan

Statistika lingkungan merupakan cabang ilmu statistika yang berfokus pada penerapan metode statistik untuk menganalisis fenomena lingkungan, seperti kualitas udara, air, dan tanah. Salah satu pendekatan penting dalam analisis spasial lingkungan adalah interpolasi, yaitu teknik untuk memperkirakan nilai suatu variabel pada lokasi yang tidak terukur berdasarkan nilai di lokasi-lokasi terdekat yang telah diketahui.

Dalam konteks ini, interpolasi digunakan untuk memetakan sebaran Indeks Standar Pencemar Udara (ISPU) di Provinsi DKI Jakarta dengan tujuan mengetahui variasi spasial tingkat pencemaran udara pada periode tertentu. Metode interpolasi yang digunakan meliputi Inverse Distance Weighting (IDW) dan Kriging, yang diimplementasikan menggunakan perangkat lunak R dengan bantuan paket sf, sp, gstat, terra, dan tmap.

Analisis ini diharapkan dapat memberikan gambaran visual mengenai pola sebaran pencemaran udara serta menjadi dasar untuk pengambilan keputusan dalam pengelolaan lingkungan perkotaan.

Analisis

library(sf)
## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
library(sp)
library(gstat)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(terra)
## terra 1.8.70
## 
## Attaching package: 'terra'
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
library(readxl)
library(raster)
## 
## Attaching package: 'raster'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
library(stars)
## Loading required package: abind
library(abind)
library(tmap)
# ===========================
# 1. IMPORT DATA
# ===========================

data <- read_xlsx("/Users/M.Fabian.R.D/Desktop/SEMESTER 5/STATISTIKA LINGKUNGAN/PERTEMUAN 9/data-indeks-standar-pencemar-udara-(ispu)-di-provinsi-dki-jakarta-2023-(1760891606680).xlsx", sheet = 2)   # ganti nama file sesuai milikmu
print(data)
## # A tibble: 5 × 4
##   stasiun_simple mean_ispu_mei2023 latitude longitude
##   <chr>                      <dbl>    <dbl>     <dbl>
## 1 Jagakarsa                   83.0    -6.33      107.
## 2 Kebon Jeruk                 55.1    -6.19      107.
## 3 Kelapa Gading               81.0    -6.16      107.
## 4 Lubang Buaya               106.     -6.29      107.
## 5 bunderan hi                 81.5    -6.19      107.
# Pastikan kolom sesuai: stasiun_simple, mean_ispu_mei2023, latitude, longitude
# ===========================
# 2. KONVERSI KE SPATIAL + SET CRS UTM 48S (EPSG:32748)
# ===========================
# Pastikan tipe datanya numerik
data$latitude <- as.numeric(data$latitude)
data$longitude <- as.numeric(data$longitude)

# Cek apakah ada NA
sum(is.na(data$latitude))
## [1] 0
sum(is.na(data$longitude))
## [1] 0
# Tampilkan baris yang bermasalah (jika ada)
data[!complete.cases(data$latitude, data$longitude), ]
## # A tibble: 0 × 4
## # ℹ 4 variables: stasiun_simple <chr>, mean_ispu_mei2023 <dbl>, latitude <dbl>,
## #   longitude <dbl>
sf_data <- st_as_sf(data, coords = c("longitude","latitude"), crs = 4326) # WGS84
sf_data <- st_transform(sf_data, 32748)  # UTM Zone 48S
sf_data
## Simple feature collection with 5 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 695815.3 ymin: 9299510 xmax: 710852.2 ymax: 9318682
## Projected CRS: WGS 84 / UTM zone 48S
## # A tibble: 5 × 3
##   stasiun_simple mean_ispu_mei2023           geometry
## * <chr>                      <dbl>        <POINT [m]>
## 1 Jagakarsa                   83.0 (701436.3 9299510)
## 2 Kebon Jeruk                 55.1 (695815.3 9315182)
## 3 Kelapa Gading               81.0 (710852.2 9318682)
## 4 Lubang Buaya               106.  (710570.5 9303923)
## 5 bunderan hi                 81.5   (701435 9315115)
# Simpan juga versi SP (untuk gstat)
sp_data <- as(sf_data, "Spatial")
# Tambahkan koordinat ke dalam data
sp_data@data$x <- coordinates(sp_data)[,1]
sp_data@data$y <- coordinates(sp_data)[,2]

# Cek hasil
head(sp_data@data)
##   stasiun_simple mean_ispu_mei2023        x       y
## 1      Jagakarsa          82.96774 701436.3 9299510
## 2    Kebon Jeruk          55.09677 695815.3 9315182
## 3  Kelapa Gading          81.03226 710852.2 9318682
## 4   Lubang Buaya         106.12903 710570.5 9303923
## 5    bunderan hi          81.51613 701435.0 9315115
# ===========================
# 3. BUAT GRID INTERPOLASI (500m)
# ===========================

# --- Gunakan extent dari data spasial
ext <- extent(sp_data)

# --- Buat grid dengan resolusi 500 m
grd <- expand.grid(
  x = seq(ext@xmin, ext@xmax, by = 500),
  y = seq(ext@ymin, ext@ymax, by = 500)
)

# --- Pastikan numeric
grd$x <- as.numeric(grd$x)
grd$y <- as.numeric(grd$y)

# --- Jadikan SpatialPixels
coordinates(grd) <- ~x + y
gridded(grd) <- TRUE
proj4string(grd) <- CRS("+init=epsg:32748")
## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.
# Cek hasil grid
plot(grd, main = "Grid 500 m - DKI Jakarta")
points(sp_data, col = "red", pch = 16)

# ===========================
# 4. IDW (Inverse Distance Weighting)
# ===========================

idw_model <- gstat::idw(mean_ispu_mei2023 ~ 1, sp_data, grd, idp = 2)
## [inverse distance weighted interpolation]
idw_raster <- rast(idw_model)
plot(idw_raster, main = "IDW - ISPU Mei 2023")

# ===========================
# 5. ORDINARY KRIGING (OK)
# ===========================

# Variogram
vgm_OK <- variogram(mean_ispu_mei2023 ~ 1, sp_data)
fit_OK <- fit.variogram(vgm_OK, model = vgm("Sph"))
## Warning in fit.variogram(vgm_OK, model = vgm("Sph")): singular model in
## variogram fit
plot(vgm_OK, fit_OK)

# Kriging
OK_model <- krige(mean_ispu_mei2023 ~ 1, sp_data, grd, model = fit_OK)
## [using ordinary kriging]
OK_raster <- rast(OK_model["var1.pred"])
plot(OK_raster, main = "Ordinary Kriging - ISPU Mei 2023")

# ===========================
# 6. UNIVERSAL KRIGING (UK)
#    menggunakan koordinat sebagai trend (contoh linier)
# ===========================

# ---- Universal Kriging ----
UK_model <- krige(mean_ispu_mei2023 ~ x + y, sp_data, grd, model = fit_OK)
## [using universal kriging]
# Ubah hasil kriging ke format stars (dibutuhkan tmap v4)
UK_stars <- st_as_stars(UK_model)

# Set mode peta interaktif / statis (opsional)
tmap_mode("plot") # pakai "view" kalau mau interaktif
## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`
# Plot peta hasil Universal Kriging
tm_shape(UK_stars) +
  tm_raster(
    col = "var1.pred",
    col.scale = tm_scale(values = "brewer.yl_or_rd"),
    col.legend = tm_legend(title = "ISPU (Prediksi)")
  ) +
  tm_shape(st_as_sf(sp_data)) +
  tm_symbols(col = "black", size = 0.5) +
  tm_title("Universal Kriging - ISPU Mei 2023") +
  tm_layout(
    frame = FALSE,
    legend.outside = TRUE
  )

# ===========================
# 7. CROSS VALIDATION
# ===========================

cv_OK <- krige.cv(mean_ispu_mei2023 ~ 1, sp_data, model = fit_OK, nfold = nrow(sp_data))
cv_OK
## class       : SpatialPointsDataFrame 
## features    : 5 
## extent      : 695815.3, 710852.2, 9299510, 9318682  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs 
## variables   : 6
## names       :        var1.pred,         var1.var,         observed,          residual,           zscore, fold 
## min values  : 75.1532258064516, 872.477887617066, 55.0967741935484, -32.8145161290323, -1.1109353602499,    1 
## max values  : 87.9112903225807, 872.477887617066, 106.129032258064,  30.9758064516129, 1.04868584878836,    5