4  CHIRPS Tie-Breaker

4.1 The Need for a Third Opinion

Chapters 2 and 3 evaluated forecasts against two observation sources:

  • ERA5: Global reanalysis (SEAS5’s “native” validation)
  • ENACTS: Station-blended satellite product (operational standard for Guatemala)

These showed consistent results: SEAS5 performs well for Primera, all models struggle with Postrera. But what if both observation sources happen to favor SEAS5? To break potential ties and increase confidence, we introduce a third independent observation source: CHIRPS.

4.2 Observational Data Sets

This analysis uses three independent observation sources:

CHIRPS (Climate Hazards Group InfraRed Precipitation with Station data) is a quasi-global rainfall dataset that blends satellite imagery with station data. Key characteristics:

Feature CHIRPS ENACTS ERA5
Resolution 0.05° (~5km) 0.05° (~5km) 0.25° (~25km)
Source Satellite + stations Satellite + stations Model reanalysis
Coverage 50°S-50°N Regional Global
Temporal 1981-present Varies 1940-present

CHIRPS and ENACTS both blend satellite and station data, but use different algorithms and station networks. This makes CHIRPS a useful independent check - if forecasts perform well against multiple observation sources, we can be more confident in the skill assessment.

Setup: Libraries and data loading 
library(tidyverse)
library(lubridate)
library(cumulus)
library(yardstick)
library(gghdx)
gghdx()

box::use(../../R/enacts)
box::use(../../R/seas5)

set.seed(42)

# Configuration
BASELINE_START <- 2000
BASELINE_END <- 2024
PRIMERA_MONTHS <- 5:8
POSTRERA_MONTHS <- 9:11
PRIMERA_ISSUED_MONTHS <- c(3, 4, 5)
POSTRERA_ISSUED_MONTHS <- c(6, 7, 8, 9)

# Helper function
calc_rp_threshold <- function(x, rp_target = 4, direction = -1) {
  x <- x[!is.na(x)]
  n <- length(x)
  if (n < 3) return(NA_real_)
  ranks <- rank(x * -direction, ties.method = "average")
  rp <- (n + 1) / ranks
  approx(rp, x, xout = rp_target, rule = 2)$y
}
Load forecast and CHIRPS data 
# Load forecasts
df_insiv <- cumulus::blob_read(
  name = "ds-aa-lac-dry-corridor/data/processed/insivumeh_special/insivumeh_special_models_zonal_seasonal_chiquimula.parquet",
  container = "projects"
)
df_seas5 <- seas5$load_seas5_seasonal()

# Load CHIRPS - extracted via GEE script
df_chirps_raw <- cumulus::blob_read(
  name = "ds-aa-lac-dry-corridor/raw/chirps/2026_cadc_drought_v3_aoi_chirps_monthly_historical.parquet",
  container = "projects"
)

# Wrangle CHIRPS to seasonal totals
df_chirps <- df_chirps_raw |>
  filter(ADM1_NAME == "Chiquimula") |>
  mutate(
    year = year(date),
    month = month(date),
    window = case_when(
      month %in% PRIMERA_MONTHS ~ "primera",
      month %in% POSTRERA_MONTHS ~ "postrera",
      TRUE ~ NA_character_
    )
  ) |>
  filter(!is.na(window), year >= BASELINE_START, year <= BASELINE_END) |>
  group_by(year, window) |>
  summarise(obs_mm = sum(value, na.rm = TRUE), .groups = "drop")

# Combine forecasts
df_fcst_all <- bind_rows(df_insiv, df_seas5) |>
  filter(year >= BASELINE_START, year <= BASELINE_END)

# Filter to operational leadtimes
df_fcst_filtered <- df_fcst_all |>
  mutate(issued_month = month(issued_date)) |>
  filter(
    (window == "primera" & issued_month %in% PRIMERA_ISSUED_MONTHS) |
    (window == "postrera" & issued_month %in% POSTRERA_ISSUED_MONTHS)
  )

# Join forecasts with CHIRPS
df_joined <- df_fcst_filtered |>
  left_join(df_chirps, by = c("year", "window")) |>
  filter(!is.na(obs_mm)) # there are not any NA - any ways

4.3 Metrics Against CHIRPS

NoteWhy Rank-Based Metrics Work Across Different Observation Sources

Because CHIRPS and ENACTS have different climatologies (CHIRPS shows ~300mm higher Primera totals), comparing raw mm errors across sources would be misleading. Instead, we focus on rank-based metrics:

  • Spearman correlation: Evaluates the full ranking across all years. If observations rank years as [driest, 2nd, 3rd, … wettest], how well does the forecast reproduce that ordering? Larger rank errors are penalized more heavily.

  • ROC-AUC: Evaluates separation between drought and non-drought years only. Do drought years consistently get lower forecasts than non-drought years? It ignores ranking within each group—scrambling the order of non-drought years doesn’t affect AUC.

When they diverge: A forecast could perfectly separate drought from non-drought (AUC=1) but scramble rankings within each group (moderate Spearman). Or it could track the overall wet-dry gradient well (high Spearman) but fail to place drought years at the bottom (low AUC).

These metrics allow fair comparison of forecast skill across observation sources with different absolute values.

Spearman & ROC-AUC continue to tell a similar story. Models are skillful in primera with SEAS5 as the dominant competitor. Postrera skill remains low, dubious, and messy. Different models win at different leadtimes with some surprising results of INSIVUMEH provided models showing stronger predictive power at greater leadtimes. SEAS5 remains competitive, but less dominant

Calculate thresholds and metrics 
# Forecast thresholds
df_fcst_thresh <- df_joined |>
  group_by(forecast_source, window, leadtime) |>
  summarise(fcst_thresh = calc_rp_threshold(value), .groups = "drop")

# CHIRPS observation thresholds
df_obs_thresh <- df_chirps |>
  group_by(window) |>
  summarise(obs_thresh = calc_rp_threshold(obs_mm), .groups = "drop")

# Join thresholds
df_analysis <- df_joined |>
  left_join(df_fcst_thresh, by = c("forecast_source", "window", "leadtime")) |>
  left_join(df_obs_thresh, by = "window") |>
  mutate(
    fcst_drought = value < fcst_thresh,
    obs_drought = obs_mm < obs_thresh
  )

# Continuous metrics
df_metrics <- df_analysis |>
  group_by(forecast_source, window, leadtime) |>
  summarise(
    spearman = cor(value, obs_mm, method = "spearman", use = "complete.obs"),
    bias_mm = mean(value - obs_mm, na.rm = TRUE),
    .groups = "drop"
  ) |> 
     mutate(
      window = fct_relevel(window, "primera", "postrera")
    )


# vectorized calc AUC
calc_auc <- function(truth, drought_score) {                                                                                                                                           
  if (length(unique(truth)) < 2) return(NA_real_)                                                                                                                           
  roc_auc_vec(truth, drought_score, event_level = "first")                                                                                                                            
}                                                                                                                                                                                                 

df_roc <- df_analysis |>
  mutate(
    truth = factor(
      obs_drought,
      levels = c(TRUE, FALSE),
      labels = c("drought", "no_drought")
      ),
    # roc_auc is just rank based, magnitude doesnt matter just invert to +drought +score
    drought_score = -value
  )

df_auc <- df_roc |>                                                                                 
  group_by(forecast_source, window, leadtime) |>                                                                                                                                      
  summarise(auc = calc_auc(truth, drought_score), .groups = "drop")

4.3.1 Spearman Correlation

Code 
create_metric_heatmap <- function(df, metric_col, metric_name, title, caption, midpoint = 0) {
  df_plot <- df |>
    group_by(window, leadtime) |>
    mutate(is_best = !!sym(metric_col) == max(!!sym(metric_col), na.rm = TRUE)) |>
    ungroup() 
 

  ggplot(df_plot, aes(x = factor(leadtime), y = forecast_source)) +
    geom_tile(aes(fill = !!sym(metric_col)), color = "white", linewidth = 0.5) +
    geom_tile(
      data = df_plot |> filter(is_best),
      fill = NA, color = "black", linewidth = 1.5
    ) +
    geom_text(aes(label = sprintf("%.2f", !!sym(metric_col))), size = 4, fontface = "bold", color = "black") +
    facet_wrap(~window) +
    scale_fill_gradient2(
      low = "#D73027", mid = "#FFFFBF", high = "#1A9850",
      midpoint = midpoint, name = metric_name
    ) +
    labs(title = title, x = "Leadtime (months)", y = NULL, caption = caption) +
    theme_minimal() +
    theme(legend.position = "right", panel.grid = element_blank(),
          plot.caption = element_text(hjust = 0))
}

create_metric_heatmap(
  df = df_metrics,
  metric_col = "spearman",
  metric_name ="Spearman ρ",
  title = "Spearman Correlation - CHIRPS",
  caption = "All forecasts validated against CHIRPS. Black border = best per leadtime/season."
)

Surprise finding: INSIVUMEH_CESM1 shows the strongest Spearman correlation for Primera at LT1 (0.72) and LT2 (0.52), outperforming SEAS5. This is the opposite of what we saw with ERA5 and ENACTS.

4.3.2 ROC-AUC

Code 
create_metric_heatmap(
  df= df_auc, 
  metric_col = "auc",
  metric_name = "AUC",
  title = "ROC-AUC - CHIRPS",
  caption = "All forecasts validated against CHIRPS. AUC > 0.7 = acceptable skill. Black border = best.",
  midpoint = 0.5
)

4.3.3 F1 Score (Binary)

For completeness, here’s the binary F1 score against CHIRPS - the same metric used in Chapter 2. Story remains consistent.

Code 
# Calculate F1
df_f1 <- df_analysis |>
  mutate(
    truth = factor(obs_drought, levels = c(TRUE, FALSE), labels = c("drought", "no_drought")),
    estimate = factor(fcst_drought, levels = c(TRUE, FALSE), labels = c("drought", "no_drought"))
  ) |>
  group_by(forecast_source, window, leadtime) |>
  summarise(
    f1 = f_meas_vec(truth, estimate, event_level = "first"),
    precision = precision_vec(truth, estimate, event_level = "first"),
    recall = recall_vec(truth, estimate, event_level = "first"),
    .groups = "drop"
  )

# F1 heatmap
df_plot <- df_f1 |>
  group_by(window, leadtime) |>
  mutate(is_best = f1 == max(f1, na.rm = TRUE)) |>
  ungroup() |> 
  mutate(
    window = fct_relevel(window, "primera","postrera")
  )

ggplot(df_plot, aes(x = factor(leadtime), y = forecast_source)) +
  geom_tile(aes(fill = f1), color = "white", linewidth = 0.5) +
  geom_tile(
    data = df_plot |> filter(is_best),
    fill = NA, color = "black", linewidth = 1.5
  ) +
  geom_text(aes(label = sprintf("%.0f%%", f1 * 100)), size = 4, fontface = "bold", color = "black") +
  facet_wrap(~window) +
  scale_fill_gradientn(
    colors = c("#D73027", "#FC8D59", "#FEE08B", "#D9EF8B", "#91CF60", "#1A9850"),
    limits = c(0, 1), labels = scales::percent, name = "F1"
  ) +
  labs(
    title = "F1 Score - CHIRPS",
    x = "Leadtime (months)", y = NULL,
    caption = "All forecasts validated against CHIRPS. Black border = best per leadtime/season."
  ) +
  theme_minimal() +
  theme(legend.position = "right", panel.grid = element_blank(),
        plot.caption = element_text(hjust = 0))

4.4 Conclusions

4.4.1 What We’ve Learned So Far

  1. Primera skill is robust: Multiple models show genuine skill (AUC > 0.7) across all three observation sources. This isn’t an artifact of one particular dataset.

  2. Postrera remains unresolved: No model shows reliable, consistent skill. CCSM4’s apparent better performance at longer leadtimes is suspicious - skill should not improve with leadtime. This inverted pattern is likely noise from the small sample (only ~6 drought events in 25 years).

4.4.2 Open Questions

The tie-breaker analysis raises more questions than it answers for Postrera:

  • Why the inverted skill pattern? CCSM4 showing better skill at LT2-3 than LT1 is backwards - forecast skill should degrade with leadtime. Is this genuine or noise?
  • Which model is actually better? With different metrics favoring different models at different leadtimes, we cannot make a confident recommendation.
  • Is the poor skill a data artifact? Could temporal trends in the observation sources be affecting our skill estimates?

4.4.3 Next Steps

To better understand the poor Postrera skill and the contradictory patterns across models, the next chapter examines temporal drift - whether systematic trends in forecasts or observations might explain some of what we’re seeing.

4.4.4 CHIRPS vs ENACTS Comparison

How different are CHIRPS and ENACTS for Chiquimula? Understanding this helps interpret why forecast skill might differ across observation sources.

Compare CHIRPS and ENACTS observations 
# Load ENACTS for comparison
df_enacts_compare <- enacts$load_enacts_seasonal("chiquimula")

df_compare <- df_chirps |>
  rename(chirps = obs_mm) |>
  left_join(
    df_enacts_compare |> select(year, window, enacts = obs_mm),
    by = c("year", "window")
  ) |>
  filter(!is.na(enacts))

# Correlation
corr_primera <- df_compare |> filter(window == "primera") |>
  summarise(r = cor(chirps, enacts)) |> pull(r)
corr_postrera <- df_compare |> filter(window == "postrera") |>
  summarise(r = cor(chirps, enacts)) |> pull(r)

df_compare |>
  ggplot(aes(x = enacts, y = chirps)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey50") +
  geom_point(alpha = 0.6, size = 3) +
  geom_smooth(method = "lm", se = TRUE, color = "#007CE1", fill = "#007CE1", alpha = 0.2) +
  facet_wrap(~str_to_title(window), scales = "free") +
  labs(
    title = "CHIRPS vs ENACTS: Same Region, Different Estimates",
    subtitle = sprintf("Primera r = %.2f, Postrera r = %.2f", corr_primera, corr_postrera),
    x = "ENACTS (mm)",
    y = "CHIRPS (mm)",
    caption = "Dashed line = 1:1 agreement. CHIRPS consistently higher for Primera."
  )

CHIRPS and ENACTS are correlated but not identical. CHIRPS tends to estimate higher rainfall for Primera (~300mm more on average). This means a forecast calibrated to one source may not match the other perfectly - making the three-source comparison a meaningful robustness check.

4.4.5 Forecast Bias

Bias (mean forecast - observed) is shown for completeness, but is not consequential for our framework since we use rank-based metrics and model-specific thresholds rather than raw precipitation values.

Code 
df_plot <- df_metrics |>
  group_by(window, leadtime) |>
  mutate(is_best = abs(bias_mm) == min(abs(bias_mm), na.rm = TRUE)) |>
  ungroup()

bias_max <- max(abs(df_plot$bias_mm), na.rm = TRUE)

ggplot(df_plot, aes(x = factor(leadtime), y = forecast_source)) +
  geom_tile(aes(fill = bias_mm), color = "white", linewidth = 0.5) +
  geom_tile(
    data = df_plot |> filter(is_best),
    fill = NA, color = "black", linewidth = 1.5
  ) +
  geom_text(aes(label = sprintf("%.0f", bias_mm)), size = 4, fontface = "bold", color = "black") +
  facet_wrap(~window) +
  scale_fill_gradient2(
    low = "#D73027", mid = "#FFFFBF", high = "#1A9850",
    midpoint = 0, limits = c(-bias_max, bias_max), name = "Bias (mm)"
  ) +
  labs(
    title = "Forecast Bias - CHIRPS",
    x = "Leadtime (months)", y = NULL,
    caption = "Units: millimeters. Bias = mean(forecast - observed). Negative = dry bias."
  ) +
  theme_minimal() +
  theme(legend.position = "right", panel.grid = element_blank(),
        plot.caption = element_text(hjust = 0))

Against CHIRPS, all models have dry bias for Primera and wet bias for Postrera. SEAS5 has the largest dry bias for Primera (~280mm under), while INSIVUMEH models are closer (~120mm under).