Kalibrasi Data GPM dengan Data BMKG

library(ggplot2); library(readxl); library(patchwork)
library(scales);  library(dplyr);   library(tidyr)
library(ggrepel)   # label titik scatter plot

# ── Baca data utama dari Excel ───────────────────────────────
xl  <- "Input_Plot_GPM_NoWMO.xlsx"
df  <- read_excel(xl, sheet = "DATA_UTAMA", skip = 3)
met <- read_excel(xl, sheet = "METRIK_VALIDASI", skip = 3)

# ── Bersihkan nama kolom & konversi tipe data ────────────────
df <- df |>
  rename(raw = `GPM_Raw
(mm)`, m4 = `M4_Baru
(mm)`,
         obs = `Obs_BMKG
(mm)`, fk = `Fk_Baru
(rata²)`) |>
  filter(!is.na(Tahun)) |>
  mutate(across(c(Tahun, raw, m4, obs, fk), as.numeric),
         tahun = as.integer(Tahun),
         err_raw = raw - obs,
         err_m4  = m4  - obs)

# ── Data referensi: HHMT dari stasiun BMKG & piksel GPM ─────
# Asumsi: df sudah memiliki kolom gnd_sy, gnd_ks, gpm_b4

df <- df |>
  mutate(
    # Fk per stasiun = BMKG / GPM_B4 (untuk setiap tahun)
    fk_sy = gnd_sy / gpm_b4,
    fk_ks = gnd_ks / gpm_b4,

    # Fk final = rata-rata stasiun yang tersedia
    # Jika salah satu NA (mis. 2010: gnd_sy = 0), gunakan yang ada
    fk_m4 = case_when(
      is.na(fk_sy) | gnd_sy == 0 ~ fk_ks,
      is.na(fk_ks)             ~ fk_sy,
      .default                 ~ mean(c(fk_sy, fk_ks))
    ),

    # Terapkan Fk ke GPM Raw piksel B4
    gpm_m4 = gpm_b4 * fk_m4
  )

# ── Ringkasan statistik ──────────────────────────────────────
cat("Median Fk:", round(median(df$fk_m4, na.rm = TRUE), 3), "
")
cat("Rentang:", round(min(df$fk_m4), 3), "–",
                round(max(df$fk_m4), 3), "
")

# ── Fungsi hitung metrik validasi ───────────────────────────
hitung_metrik <- function(sim, obs) {
  n     <- length(sim)
  bias  <- sim - obs
  R_val <- cor(obs, sim)
  NSE   <- 1 - sum(bias^2) / sum((obs - mean(obs))^2)
  KGE   <- 1 - sqrt((R_val-1)^2 +
                        (sd(sim)/sd(obs)-1)^2 +
                        (mean(sim)/mean(obs)-1)^2)
  PBIAS <- sum(bias) / sum(obs) * 100
  RMSE  <- sqrt(mean(bias^2))
  MAE   <- mean(abs(bias))
  data.frame(R=R_val, NSE=NSE, KGE=KGE, PBIAS=PBIAS, RMSE=RMSE, MAE=MAE)
}

# ── Bandingkan M0 (Raw) vs M4 ────────────────────────────────
metrik_raw <- hitung_metrik(df$raw, df$obs)
metrik_m4  <- hitung_metrik(df$m4,  df$obs)

bind_rows(
  mutate(metrik_raw, Metode = "GPM Raw"),
  mutate(metrik_m4,  Metode = "GPM M4")
) |> select(Metode, R, NSE, KGE, PBIAS, RMSE, MAE)

# ── Reshape ke long format untuk ggplot ──────────────────────
df_long <- df |>
  select(tahun, raw, m4, obs) |>
  pivot_longer(c(raw, m4, obs), names_to = "seri", values_to = "hhmt") |>
  mutate(seri = factor(seri,
    levels = c("obs", "m4", "raw"),
    labels = c("Obs BMKG", "GPM M4 (Koreksi)", "GPM Raw")))

p1 <- ggplot(df_long, aes(x = tahun, y = hhmt,
       color = seri, linetype = seri, linewidth = seri)) +
  geom_line(na.rm = TRUE) +
  geom_point(data = . |> filter(seri == "Obs BMKG"),
             size = 1.8, shape = 21) +
  geom_label(data = df |> filter(m4 > 200),
             aes(x = tahun, y = m4 + 8,
                 label = paste0(tahun, "
", m4, "mm")),
             inherit.aes = FALSE, size = 2.8, color = "#1A5276") +
  scale_color_manual(
    values = c("Obs BMKG"="#B45309", "GPM M4 (Koreksi)"="#1A5276",
               "GPM Raw"="#B91C1C"), name = NULL) +
  scale_linetype_manual(
    values = c("Obs BMKG"="dashed", "GPM M4 (Koreksi)"="solid",
               "GPM Raw"="solid"), name = NULL) +
  scale_y_continuous(
    limits = c(0, 280),
    labels = function(x) paste0(x, " mm")) +
  labs(
    title    = "Time Series HHMT Tahunan",
    subtitle = "GPM Raw vs GPM M4 vs Observasi BMKG | 1998–2025",
    x = NULL, y = "HHMT (mm)",
    caption  = "Sumber: GPM IMERG V07 (NASA) | BMKG Syamsuddin Noor + BMKG Kalsel"
  ) +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"))

# ── Scatter Plot: M4 vs Obs ───────────────────────────────────
p_scatter <- ggplot(df, aes(x = obs, y = m4)) +
  geom_abline(slope = 1, intercept = 0,
              linetype = "dashed", color = "#94A3B8") +
  geom_point(color = "#1A5276", size = 3, shape = 21) +
  ggrepel::geom_text_repel(aes(label = tahun),
              size = 2.5, color = "#475569") +
  geom_smooth(method = "lm", se = TRUE,
              color = "#1A5276", fill = "#1A527622") +
  coord_equal() +
  labs(title = "Scatter: GPM M4 vs Obs BMKG",
       x = "Observasi BMKG (mm)", y = "GPM M4 (mm)")

# ── Gabung & simpan ──────────────────────────────────────────
fig_final <- p1 / p_scatter +
  plot_annotation(
    title   = "Koreksi Bias GPM IMERG V07 — WPG Banjarmasin",
    caption = "SNI 2415:2016 | Teutschbein & Seibert (2012)")

ggsave("hasil_kalibrasi_gpm.png", fig_final,
       width = 36, height = 30, units = "cm", dpi = 300)