Script ini digunakan untuk membuat laporan bulanan IPU (Indeks Pencemar Udara) berdasarkan data PM2.5 dan data meteorologi (arah dan kecepatan angin).
Pastikan beberapa hal berikut sebelum dijalankan:

  1. Unduh data PM2.5 dari DBKU untuk 1 lokasi dengan rentang waktu tanggal terakhir bulan sebelumnya hingga tanggal terakhir bulan analisis (contoh: 31 Juli – 31 Agustus).
  2. Unduh data meteorologis dari BMKGSoft (unsur: arah angin dan kecepatan angin) dengan rentang waktu sama.
  3. Pastikan zona waktu lokasi sudah sesuai.
  4. Periksa pengaturan separator file CSV (sep = "," atau ";").

1. Pemuatan library

rm(list = ls())
library(dplyr)
library(tidyr)
library(openair)

2. Baca data PM2.5

# Jendela pop-up akan muncul. Pilih file CSV data PM2.5
data <- read.csv(file.choose(), header = TRUE, sep = ";", na.string = "NA")

3. Edit format waktu data PM2.5

# Mengambil waktu (HH:MM) dari kolom JAM
time <- substr(as.POSIXct(sprintf("%s", data$JAM), format='%H'), 12, 16) 
# Menggabungkan TANGGAL dan Waktu, lalu mengkonversi ke objek waktu di UTC
date <- as.POSIXct(paste(data$TANGGAL, time), format= "%d/%m/%Y %H:%M", tz = "UTC")

4. Penyesuaian zona waktu

attr(date, 'tzone') <- "Asia/Jakarta"  
# Contoh: WIB = "Asia/Jakarta", WITA = "Asia/Makassar, WIT = "Asia/Jayapura"

5. Gabung data PM2.5

data <- data.frame(date, data$PM25)
colnames(data) <- c("date", "pm2.5")

6. Lengkapi time series PM2.5

data <- data %>%
  mutate(date = as.POSIXct(date)) %>%
  complete(date = seq.POSIXt(min(date), max(date), by="hour"))
head(data)

7. Baca data meteorologi

# Jendela pop-up akan muncul. Pilih file CSV data meteorologi
data_angin <- read.csv(file.choose(), header = TRUE, sep = ",", na.string = "NA")

date2 <- as.POSIXct(data_angin$DATA.TIMESTAMP, format = "%Y-%m-%d %H:%M", tz = "UTC")
attr(date2, 'tzone') <- "Asia/Jakarta"
data_angin <- data.frame(date2, data_angin$WIND.DIR.DEG.DD, data_angin$WIND.SPEED.FF)
colnames(data_angin) <- c("date2", "wd", "ws")

8. Lengkapi time series data meteorologi

data_angin <- data_angin %>%
  mutate(date2 = as.POSIXct(date2)) %>%
  complete(date2 = seq.POSIXt(min(date2), max(date2), by="hour"))

9. Gabungkan data PM2.5 dan data meteorologi

n <- nrow(data)
data <- data.frame(
  date = data$date,
  pm2.5 = data$pm2.5,
  wd = data_angin$wd[1:n], # CATATAN: Penggabungan ini rentan error jika panjang data berbeda
  ws = data_angin$ws[1:n]
)

10. Koreksi zona waktu

data <- tail(data, -17) # WIB = -17, WITA = -16, WIT = -15
data <- head(data, -7)  # WIB = -7, WITA = -8, WIT = -9
head(data)

11. Visualisasi data

plot(
  data$date, data$pm2.5,
  type = "l",
  col = "steelblue",
  main = "Konsentrasi PM2.5 - Stasiun Kemayoran (Januari 2024)",
  xlab = "Tanggal",
  ylab = expression("Konsentrasi PM2.5 ("*mu*"g/m"^3*")")
)

windRose(data, key.footer = "knot")

pollutionRose(data, pollutant = "pm2.5",
              breaks = c(0, 15.5, 55.4, 150.4, 250.4, 500),
              cols = c("#00FF00", "#0033FF", "#FFFF00", "#FF0000", "#000000"))
              
calendarPlot(data, pollutant = "pm2.5",
             main = "Kalender Level Konsentrasi PM2.5 - Staklim Banjarbaru (Desember 2023)",
             breaks = c(0, 15.5, 55.4, 150.4, 250.4, 500),
             labels = c("Baik", "Sedang", "Tidak sehat", "Sangat tidak sehat", "Berbahaya"),
             cols = c("#00FF00", "#0033FF", "#FFFF00", "#FF0000", "#000000"))
             
timeVariation(data, pollutant = "pm2.5", ylab = "PM2.5 (µg/m³)")

12. Summary statistik

nmax <- which.max(data$pm2.5)
nmin <- which.min(data$pm2.5)
max_val <- data$pm2.5[nmax]
min_val <- data$pm2.5[nmin]
mean_val <- mean(data$pm2.5, na.rm = TRUE)
cat("Konsentrasi Maksimum:", max_val, "µg/m³ pada", data$date[nmax], "\n")
cat("Konsentrasi Minimum:", min_val, "µg/m³ pada", data$date[nmin], "\n")
cat("Rata-rata Konsentrasi:", round(mean_val, 2), "µg/m³\n")