Visualisasi Data Hujan
dengan R & ggplot2

# Instalasi (jalankan sekali saja)
install.packages("tidyverse")  # termasuk ggplot2, dplyr, lubridate
install.packages("readxl")     # baca file Excel
install.packages("scales")     # format sumbu
install.packages("patchwork")  # multi-panel layout

# Load package di setiap sesi
library(tidyverse)
library(readxl)
library(scales)
library(patchwork)

library(tidyverse)
library(readxl)

# ── 1. Import sheet data_harian ─────────────────────────────
df_raw <- read_excel(
  "sta_hujan_cemara_bandung_2000_2025.xlsx",
  sheet = "data_harian"
)

# ── 2. Bersihkan & buat kolom turunan ───────────────────────
df <- df_raw |>
  rename(
    tanggal  = Tanggal,
    rr_raw   = RR,
    rr_mm    = `Curah_hujan (mm)`
  ) |>
  mutate(
    rr_mm   = if_else(rr_raw == 8888, NA_real_, rr_raw),  # ganti 8888 → NA
    tahun   = year(tanggal),
    bulan   = month(tanggal, label = TRUE, abbr = TRUE),
    bln_num = month(tanggal),
    hujan   = if_else(rr_mm > 0, 1L, 0L, missing = 0L)  # flag hari hujan
  )

# ── 3. Cek hasil ────────────────────────────────────────────
glimpse(df)
summary(df$rr_mm)

# ── Rekap bulanan: rata-rata curah hujan per bulan (26 tahun) ──
df_bulanan <- df |>
  group_by(tahun, bln_num, bulan) |>
  summarise(
    total_mm   = sum(rr_mm, na.rm = TRUE),
    maks_mm    = max(rr_mm, na.rm = TRUE),
    hari_hujan = sum(hujan, na.rm = TRUE),
    .groups    = "drop"
  )

# ── Rekap tahunan ────────────────────────────────────────────
df_tahunan <- df_bulanan |>
  group_by(tahun) |>
  summarise(
    total_mm   = sum(total_mm),
    maks_mm    = max(maks_mm),
    hari_hujan = sum(hari_hujan)
  )

# ── Rata-rata bulanan (lintas tahun) ─────────────────────────
df_avg_bln <- df_bulanan |>
  group_by(bln_num, bulan) |>
  summarise(avg_mm = mean(total_mm), .groups = "drop")

# ── Versi dasar: 3 baris ────────────────────────────────────
ggplot(df_avg_bln, aes(x = bulan, y = avg_mm)) +
  geom_col(fill = "#378add") +
  labs(title = "Rata-rata Curah Hujan Bulanan",
       y = "mm", x = NULL)

# ── Warna per kategori musim ─────────────────────────────────
warna_bln <- c(
  "Jan"="#185fa5", "Feb"="#185fa5", "Mar"="#185fa5",
  "Apr"="#378add", "Mei"="#378add",
  "Jun"="#85b7eb", "Jul"="#85b7eb", "Agt"="#b5d4f4", "Sep"="#85b7eb",
  "Okt"="#378add", "Nov"="#185fa5", "Des"="#185fa5"
)

p_bulanan <- ggplot(df_avg_bln, aes(x = bulan, y = avg_mm, fill = bulan)) +

  # ① Bar utama
  geom_col(show.legend = FALSE, width = 0.75) +

  # ② Garis rata-rata tahunan
  geom_hline(yintercept = 184, linetype = "dashed",
             color = "#ef9f27", linewidth = 0.8) +

  # ③ Label di atas bar
  geom_text(aes(label = round(avg_mm, 0)),
            vjust = -0.4, size = 3, color = "#495057") +

  # ④ Anotasi garis rata-rata
  annotate("text", x = 11.5, y = 192,
           label = "Rata-rata 184 mm",
           size = 3, color = "#ef9f27", fontface = "bold") +

  # ⑤ Warna & skala
  scale_fill_manual(values = warna_bln) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.12)),
                     labels = label_number(suffix = " mm")) +

  # ⑥ Label
  labs(
    title   = "Rata-rata Curah Hujan Bulanan",
    subtitle= "Sta. Hujan Cemara, Bandung · 2000–2025",
    x = NULL, y = "Curah Hujan (mm)",
    caption = "Sumber: BMKG Bandung"
  ) +

  # ⑦ Tema
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(face = "bold"),
    plot.subtitle = element_text(color = "grey50", size = 10),
    panel.grid.major.x = element_blank()
  )

print(p_bulanan)

# ── Tambahkan kolom kategori tahun ──────────────────────────
df_tahunan <- df_tahunan |>
  mutate(
    kategori = case_when(
      total_mm > 3000 ~ "Ekstrem Basah",
      total_mm < 1800 ~ "Kering",
      .default        = "Normal"
    )
  )

# ── Bar chart tahunan ────────────────────────────────────────
p_tahunan <- ggplot(df_tahunan, aes(x = tahun, y = total_mm, fill = kategori)) +

  geom_col(width = 0.7) +
  geom_hline(yintercept = 2215, linetype = "dashed",
             color = "#1d9e75", linewidth = 0.9) +

  # Label hanya untuk tahun ekstrem
  geom_text(
    data = ~filter(.x, kategori != "Normal"),
    aes(label = comma(total_mm)),
    vjust = -0.4, size = 2.8, fontface = "bold"
  ) +

  scale_fill_manual(
    values = c("Ekstrem Basah" = "#a32d2d",
               "Kering"        = "#ef9f27",
               "Normal"        = "#378add")
  ) +
  scale_y_continuous(labels = comma,
                     expand = expansion(mult = c(0.02, 0.12))) +
  scale_x_continuous(breaks = seq(2000, 2025, 5)) +
  labs(
    title    = "Curah Hujan Tahunan 2000–2025",
    subtitle = "Sta. Hujan Cemara, Bandung | Garis hijau = rata-rata 2.215 mm",
    x = NULL, y = "Curah Hujan (mm)",
    fill = "Kategori", caption = "Sumber: BMKG Bandung"
  ) +
  theme_minimal(base_size = 12) +
  theme(panel.grid.major.x = element_blank(),
        plot.title = element_text(face = "bold"))

print(p_tahunan)

# ── Heatmap curah hujan bulanan × tahunan ───────────────────
p_heatmap <- ggplot(df_bulanan,
  aes(x = factor(bln_num), y = factor(tahun), fill = total_mm)) +

  # ① Kotak warna
  geom_tile(color = "white", linewidth = 0.3) +

  # ② Label nilai di tiap sel
  geom_text(aes(label = round(total_mm, 0)),
            size = 2.5, color = ifelse(df_bulanan$total_mm > 300, "white", "#333")) +

  # ③ Skala warna: putih (kering) → biru tua (basah)
  scale_fill_gradient2(
    low      = "#f0f7ff",
    mid      = "#378add",
    high     = "#0a2540",
    midpoint = 200,
    name     = "mm"
  ) +

  # ④ Label sumbu X: nama bulan
  scale_x_discrete(labels = c("Jan","Feb","Mar","Apr","Mei","Jun",
                                "Jul","Agt","Sep","Okt","Nov","Des")) +

  labs(
    title    = "Heatmap Curah Hujan Bulanan (mm)",
    subtitle = "Sta. Hujan Cemara, Bandung · 2000–2025",
    x = "Bulan", y = "Tahun",
    caption  = "Sumber: BMKG Bandung"
  ) +
  theme_minimal(base_size = 11) +
  theme(
    plot.title    = element_text(face = "bold"),
    axis.text.y   = element_text(size = 8),
    panel.grid    = element_blank(),
    legend.position = "right"
  )

print(p_heatmap)

library(patchwork)

# ── Chart hari hujan tahunan ─────────────────────────────────
p_harihujan <- ggplot(df_tahunan, aes(x = tahun, y = hari_hujan)) +
  geom_line(color = "#639922", linewidth = 1.1) +
  geom_point(color = "#639922", size = 2.5, fill = "white", shape = 21) +
  geom_hline(yintercept = mean(df_tahunan$hari_hujan),
             linetype = "dashed", color = "#639922", alpha = 0.5) +
  scale_x_continuous(breaks = seq(2000, 2025, 5)) +
  labs(title = "Hari Hujan per Tahun", x = NULL, y = "Hari") +
  theme_minimal(base_size = 11)

# ── Gabungkan: tahunan (atas) / hari hujan (bawah) ──────────
figure_final <- p_tahunan / p_harihujan +
  plot_annotation(
    title   = "Ringkasan Hidrologi Sta. Hujan Cemara, Bandung",
    subtitle= "Periode Pengamatan 2000–2025",
    caption = "Sumber: BMKG Bandung · Analisis: irpanchumaedi.com",
    theme   = theme(plot.title = element_text(face = "bold", size = 14))
  )

print(figure_final)

# ── PNG resolusi tinggi (untuk laporan/presentasi) ───────────
ggsave(
  filename = "curah_hujan_bulanan_cemara.png",
  plot     = p_bulanan,
  width    = 10,        # lebar dalam inci
  height   = 6,         # tinggi dalam inci
  dpi      = 300,       # resolusi (300 dpi = print quality)
  bg       = "white"    # background putih
)

# ── PDF vektoral (untuk jurnal/laporan teknis) ───────────────
ggsave(
  filename = "heatmap_cemara.pdf",
  plot     = p_heatmap,
  width    = 14, height = 8,
  device   = "pdf"
)

# ── Multi-panel figure ───────────────────────────────────────
ggsave(
  filename = "ringkasan_hidrologi_cemara.png",
  plot     = figure_final,
  width    = 12, height = 9, dpi = 300
)

# ── Boxplot curah hujan harian per bulan (hanya hari hujan) ─
p_boxplot <- df |>
  filter(rr_mm > 0) |>                       # hanya hari hujan

  ggplot(aes(x = bulan, y = rr_mm, fill = bulan)) +

  # ① Boxplot utama
  geom_boxplot(
    show.legend  = FALSE,
    outlier.shape  = 21,
    outlier.fill   = "#e24b4a",
    outlier.color  = "white",
    outlier.size   = 1.8,
    outlier.alpha  = 0.7,
    width          = 0.6
  ) +

  # ② Titik rata-rata (diamond)
  stat_summary(fun = mean, geom = "point",
               shape = 23, size = 3,
               fill = "#ef9f27", color = "white") +

  # ③ Garis threshold hujan lebat (≥50mm/hari, SNI)
  geom_hline(yintercept = 50, linetype = "dashed",
             color = "#e24b4a", linewidth = 0.7) +
  annotate("text", x = 11.6, y = 53,
           label = "≥50 mm = Hujan Lebat (BMG)",
           size = 3, color = "#e24b4a") +

  scale_fill_manual(values = setNames(
    colorRampPalette(c("#b5d4f4", "#185fa5"))(12),
    levels(df$bulan)
  )) +
  scale_y_continuous(labels = label_number(suffix = " mm")) +
  labs(
    title    = "Distribusi Curah Hujan Harian per Bulan",
    subtitle = "Sta. Cemara, Bandung 2000–2025 | ◆ = rata-rata | ● merah = outlier ekstrem",
    x = NULL, y = "Curah Hujan (mm/hari)",
    caption  = "Hanya hari dengan curah hujan > 0 mm | Sumber: BMKG Bandung"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"),
        panel.grid.major.x = element_blank())

print(p_boxplot)

# ── 26 mini chart bulanan — satu baris facet_wrap ───────────
p_facet <- ggplot(df_bulanan,
  aes(x = bln_num, y = total_mm, fill = total_mm)) +

  geom_col(show.legend = FALSE, width = 0.85) +

  # ← INI YANG TIDAK BISA DILAKUKAN DI EXCEL
  facet_wrap(~ tahun, ncol = 6) +

  scale_fill_gradient(low = "#b5d4f4", high = "#0a2540") +
  scale_x_continuous(breaks = c(1,6,12),
                     labels = c("Jan","Jun","Des")) +
  scale_y_continuous(labels = label_number(suffix="mm", scale=0.001, accuracy=0.1)) +
  labs(
    title    = "Curah Hujan Bulanan per Tahun — Sta. Cemara Bandung",
    subtitle = "26 chart sekaligus | Warna lebih gelap = curah hujan lebih tinggi",
    x = NULL, y = NULL,
    caption  = "Sumber: BMKG Bandung"
  ) +
  theme_minimal(base_size = 9) +
  theme(
    strip.text       = element_text(face = "bold", size = 8),
    panel.grid.minor = element_blank(),
    plot.title       = element_text(face = "bold")
  )

# Simpan ukuran besar agar semua panel terbaca
ggsave("facet_bulanan_cemara.png", p_facet,
       width = 16, height = 12, dpi = 300)

# ── Tren curah hujan tahunan + CI 95% ───────────────────────
p_trend <- ggplot(df_tahunan, aes(x = tahun, y = total_mm)) +

  # ① Bar abu-abu sebagai latar
  geom_col(fill = "#85b7eb", alpha = 0.5, width = 0.7) +

  # ② Garis tren + pita CI 95% ← tidak bisa di Excel
  geom_smooth(
    method    = "lm",          # linear model
    color     = "#e24b4a",
    fill      = "#e24b4a",
    alpha     = 0.15,           # transparansi pita CI
    linewidth = 1.2,
    se        = TRUE            # tampilkan CI
  ) +

  # ③ Hitung dan tampilkan nilai tren
  annotate("text",
    x = 2002, y = 3500,
    label = paste("Tren:",
      round(coef(lm(total_mm ~ tahun, df_tahunan))[2], 1),
      "mm/tahun"),
    size = 4, color = "#e24b4a", fontface = "bold", hjust = 0
  ) +

  scale_x_continuous(breaks = seq(2000, 2025, 5)) +
  scale_y_continuous(labels = comma,
                     limits = c(0, 4200)) +
  labs(
    title    = "Tren Curah Hujan Tahunan 2000–2025",
    subtitle = "Pita merah = Confidence Interval 95% | Area pita = zona ketidakpastian",
    x = NULL, y = "Curah Hujan Tahunan (mm)",
    caption  = "Metode: Linear Model (lm) | Sumber: BMKG Bandung"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

print(p_trend)

library(ggridges)   # install.packages("ggridges")
library(forcats)   # sudah termasuk dalam tidyverse

# ── Ridgeline: distribusi intensitas hujan harian per bulan ─
p_ridge <- df |>
  filter(rr_mm > 0) |>
  mutate(bulan = fct_rev(bulan)) |>   # Jan di atas

  ggplot(aes(x = rr_mm, y = bulan,
             fill = after_stat(x))) +   # warna = intensitas

  # ① Kurva distribusi bergradien
  geom_density_ridges_gradient(
    scale          = 2.5,    # tingkat tumpang tindih
    rel_min_height = 0.01,
    color          = "white",
    linewidth      = 0.4
  ) +

  # ② Gradien warna: biru muda (ringan) → biru tua (lebat)
  scale_fill_gradientn(
    colors = c("#e6f1fb", "#378add", "#0a2540", "#e24b4a"),
    values = c(0, 0.3, 0.7, 1),
    guide  = "none"
  ) +

  # ③ Garis referensi: 50mm (hujan lebat)
  geom_vline(xintercept = 50, linetype = "dashed",
             color = "#ef9f27", linewidth = 0.8) +

  scale_x_continuous(
    limits = c(0, 150),
    labels = label_number(suffix = " mm")
  ) +
  labs(
    title    = "Distribusi Intensitas Hujan Harian per Bulan",
    subtitle = "Sta. Cemara, Bandung 2000–2025 | Garis oranye = 50mm (Hujan Lebat)",
    x = "Intensitas (mm/hari)", y = NULL,
    caption  = "Hanya hari hujan (RR > 0) | Sumber: BMKG Bandung"
  ) +
  theme_ridges(grid = FALSE) +
  theme(plot.title = element_text(face = "bold"))

print(p_ridge)

# ── 1. Buat fungsi yang bisa dipakai ulang ──────────────────
buat_chart_hujan <- function(df, nama_stasiun, warna = "#378add") {

  df_avg <- df |>
    filter(rr_mm != 8888) |>
    group_by(tahun, bln_num) |>
    summarise(total = sum(rr_mm, na.rm = TRUE), .groups = "drop") |>
    group_by(bln_num) |>
    summarise(avg = mean(total))

  ggplot(df_avg, aes(x = bln_num, y = avg)) +
    geom_col(fill = warna, width = 0.75) +
    geom_hline(yintercept = mean(df_avg$avg),
               linetype = "dashed", color = "#ef9f27") +
    scale_x_continuous(breaks = 1:12,
      labels = c("Jan","Feb","Mar","Apr","Mei","Jun",
                 "Jul","Agt","Sep","Okt","Nov","Des")) +
    labs(title = paste("Rata-rata Bulanan —", nama_stasiun),
         x = NULL, y = "mm") +
    theme_minimal(base_size = 11) +
    theme(panel.grid.major.x = element_blank())
}

# ── 2. Daftar semua stasiun ──────────────────────────────────
stasiun_list <- tibble(
  nama  = c("Cemara", "Lembang", "Padalarang", "Cimahi",
             "Margahayu", "Baleendah"),
  file  = paste0("data/", c("cemara","lembang","padalarang",
                              "cimahi","margahayu","baleendah"), ".xlsx"),
  warna = c("#185fa5","#388e3c","#ba7517",
             "#7b1fa2","#c62828","#00838f")
)

# ── 3. Loop: baca data + buat + simpan chart ─────────────────
dir.create("output_charts", showWarnings = FALSE)

stasiun_list |>
  pwalk(function(nama, file, warna) {

    df_sta <- read_excel(file, sheet = "data_harian")

    p <- buat_chart_hujan(df_sta, nama, warna)

    ggsave(
      filename = paste0("output_charts/", tolower(nama), ".png"),
      plot     = p,
      width = 10, height = 6, dpi = 300, bg = "white"
    )
    message("✓ Selesai: ", nama)
  })

# Output: 6 file PNG tersimpan otomatis di folder output_charts/