Bioclimatic Indices for Public Health

About this tutorial This document demonstrates the full climasus4r pipeline for generating bioclimatic indices from hourly INMET data: download → gap-filling → computation of 15 indicators → visualization and regional categorization of thermal risk.

Theoretical Foundations

Why bioclimatic indices?

Air temperature (T) alone does not adequately represent human thermal stress. The thermoregulatory system responds to a complex energy balance that integrates:

Air temperature (T, °C)
Relative humidity (RH, %): affects the efficiency of sweat evaporation.
Wind speed (WS, m/s): modifies the convective boundary layer.
Solar radiation (SR, W/m²): adds radiant heat load.

Indices such as WBGT, UTCI, and PET translate this interaction into a single scalar value with direct physiological meaning, and are standard tools in occupational health, heat wave surveillance, and urban planning.

📊 Bioclimatic and Thermal Stress Indicators

Code	Name	Input Variables	Domain of Validity	Main Reference
`wbgt`	Wet-Bulb Globe Temperature	T, RH, SR, WS	General	Liljegren et al. (2008)
`hi`	Heat Index	T, RH	T ≥ 26.7 °C	Rothfusz (1990); NWS
`thi`	Temperature-Humidity Index	T, RH	General	Thom (1959)
`wcet`	Wind Chill Equivalent Temp. (CA)	T, WS	T ≤ 0 °C	Environment Canada (2001)
`wct`	Wind Chill Temp. (NWS)	T, WS	T ≤ 10 °C	NWS (2001)
`et`	Effective Temperature	T, RH, WS	General	Missenard (1937) apud Auliciems (2007)
`utci`	Universal Thermal Climate Index	T, RH, WS, SR	-60 to +50 °C	Bröde et al. (2012)
`pet`	Physiological Equivalent Temp.	T, RH, WS, SR	General	Matzarakis et al. (1999)
`cdd`	Cooling Degree Days	T	T > T_base	ASHRAE (2021)
`hdd`	Heating Degree Days	T	T < T_base	ASHRAE (2021)
`gdd`	Growing Degree Days	T	T_base < T < T_upper	McMaster & Wilhelm (1997)
`diurnal_range`	Diurnal Thermal Amplitude	T (min, max)	General	-
`vapor_pressure`	Actual Vapor Pressure	T, RH	General	Magnus-Tetens (1930)
`heat_stress_risk`	Thermal Stress Risk	WBGT, UTCI	Categorical	OSHA/ACGIH (2017)
`koppen_humidity`	Humidity Class (Köppen)	RH	Categorical	Köppen (1936)

Mathematical Formulas

1. WBGT (Wet-Bulb Globe Temperature)

Model implemented for outdoor environments under direct solar load.

WBGT_out = 0.7 × T_nwb + 0.2 × T_g + 0.1 × T_db

Where:

T_db = Dry-bulb temperature (°C).
T_nwb = Natural wet-bulb temperature (°C). For numerical robustness, the average of two estimates is used:
- Liljegren: T_nwb_lilj = T * atan(0.16 * sqrt(e_a + 0.1)) + 3.0
- Bernard & Pourmoghani: T_nwb_bp = T + 0.33 * (RH/100) * exp(0.0514 * T) - 4.0
- T_nwb = (T_nwb_lilj + T_nwb_bp) / 2
- Where e_a is the actual vapor pressure (kPa), calculated from T and RH.
T_g = Globe temperature (°C), which models the radiant load:
- T_g = T + 0.0144 * max(G_0, 0)^0.6 / WS^0.2 - 2.0
- G_0 is the solar irradiance in W/m², converted from SR (kJ/m²) via SR / 3.6.

Weighting: The coefficients 0.7, 0.2, and 0.1 represent the relative contribution of evaporative cooling (wet bulb), radiant load (globe), and sensible heat (dry bulb) to thermal stress in outdoor workers.

2. Heat Index (HI)

Implementation of the polynomial regression model from the National Weather Service (NWS), with adjustments for extreme humidity conditions.

HI_°F = −42.379 + 2.04901523·T + 10.14333127·RH − 0.22475541·T·RH − 0.00683783·T² − 0.05481717·RH² + 0.00122874·T²·RH + 0.00085282·T·RH² − 0.00000199·T²·RH²

The final value is converted to Celsius: HI = (HI°F - 32) * 5/9.

Validity Threshold Adapted by Biome: * Standard: T ≥ 26.7 °C. * Tropical regions (Amazon, Atlantic Forest): T ≥ 24.0 °C, due to the acclimatization of local populations.

Post-Calculation Adjustments: * If RH < 13% and 26.7 ≤ T ≤ 44.4 °C: HI = HI - ((13 - RH)/4) * sqrt((17 - |T - 95|) / 17). * If RH > 85% and 26.7 ≤ T ≤ 30.6 °C: HI = HI + ((RH - 85)/10) * ((87 - T)/5).

3. Universal Thermal Climate Index (UTCI)

Polynomial approximation of the Fiala reference model, optimized for large-scale operational applications.

Mean Radiant Temperature (T_mrt): T_mrt = T + 0.07 * G_0 - 1.5 (approximation for shortwave radiation).
Vapor Pressure (VP, kPa): VP = 0.6108 * exp(17.27 * T / (T + 237.3)) * (RH / 100) @Tetens1930.
Wind Speed at 10 m (WS_10m): WS_10m = max(WS * 1.5, 0.5) (conversion from the standard 2 m height to 10 m).

UTCI (°C) = T + 0.048 * ((VP1000) - 1000)/100 - 0.29 ln(WS_10m + 0.1) + 0.016 * ΔT_mrt + 0.010 * (ΔT_mrt²/10) - 6e-5 * ΔT_mrt² * WS_10m + δ_regional

Where ΔT_mrt = T_mrt - T. The term δ<sub>regional</sub> is an empirical correction calibrated for Brazilian biomes.

4. Degree Days (CDD, HDD, GDD)

CDD (Cooling Degree Days): CDD = max(0, T - T_base_cdd). Regional base: 18–20 °C.
HDD (Heating Degree Days): HDD = max(0, T_base_hdd - T). Regional base: 14–18 °C.
GDD (Growing Degree Days): GDD = max(0, min(T, T_upper) - T_base_gdd). Regional base: 8–15 °C. Regional upper limit: 30–35 °C.

Full Pipeline

Step 1 · INMET Import

library(climasus4r)
library(dplyr)
# ── Pipeline para 4 biomas representativos ───────────────────────────────────
# Amazônia (AM), Caatinga (CE), Mata Atlântica (SP), Pampa (RS)

df_am <- sus_climate_inmet(
  years   = 2023,
  uf      = "AM",          # Amazônia
  parallel = TRUE,
  workers  = 4,
  verbose  = TRUE
)

df_ce <- sus_climate_inmet(years = 2023, uf = "CE")   # Caatinga
df_sp <- sus_climate_inmet(years = 2023, uf = "SP")   # Mata Atlântica
df_rs <- sus_climate_inmet(years = 2023, uf = "RS")   # Pampa

# Inspecionar estrutura
head(df_am)

Step 2 · Gap-Filling

# ── XGBoost gap-filling — todas as variáveis ─────────────────────────────────
df_am_filled <- df_am |>
  sus_climate_fill_inmet(
    target_var        = "all",
    quality_threshold = 0.4,
    parallel          = TRUE,
    workers           = 4,
    verbose           = TRUE
  )

df_ce_filled <- df_ce |> sus_climate_fill_inmet(target_var = "all")
df_sp_filled <- df_sp |> sus_climate_fill_inmet(target_var = "all")
df_rs_filled <- df_rs |> sus_climate_fill_inmet(target_var = "all")

Step 3 · Computing the Indicators

# ── Carregar função ───────────────────────────────────────────────────────────
# Em produção: já disponível via library(climasus4r)
# Aqui: source do arquivo do pacote
# source("sus_climate_compute_indicators.R")  # descomente se necessário

# ── Configuração de indicadores por bioma ────────────────────────────────────

# Amazônia — adaptação tropical, foco em estresse por calor e umidade
ind_am <- df_am_filled |>
  sus_climate_compute_indicators(
    region           = "amazon",
    confidence_flags = TRUE,
    verify_physics   = TRUE,
    verbose          = FALSE
  )

# Caatinga — semiárido, limiares rebaixados para RH
ind_ce <- df_ce_filled |>
  sus_climate_compute_indicators(
    region           = "caatinga",
    confidence_flags = TRUE,
    verbose          = FALSE
  )

# Mata Atlântica — subtropical úmido
ind_sp <- df_sp_filled |>
  sus_climate_compute_indicators(
    region           = "atlantic_forest",
    confidence_flags = TRUE,
    verbose          = FALSE
  )

# Pampa — subtropical temperado, wind chill relevante no inverno
ind_rs <- df_rs_filled |>
  sus_climate_compute_indicators(
    region           = "pampa",
    confidence_flags = TRUE,
    verbose          = FALSE
  )

# ── Combinar todos os biomas ──────────────────────────────────────────────────
ind_all <- bind_rows(
  mutate(ind_am, bioma = "Amazônia",      regiao = "Norte"),
  mutate(ind_ce, bioma = "Caatinga",      regiao = "Nordeste"),
  mutate(ind_sp, bioma = "Mata Atlântica",regiao = "Sudeste"),
  mutate(ind_rs, bioma = "Pampa",         regiao = "Sul")
)

glimpse(ind_all |> select(date, bioma, wbgt_c, hi_c, thi_c, utci_c,
                           pet_c, heat_stress_risk, koppen_humidity))

Visualizations

4.1 Overview of the Indicators

# ── Distribuições dos principais índices por bioma ───────────────────────────
library(tidyr)
library(ggplot2)
library(ggridges)
library(patchwork)
library(scales)


bioma_colors <- c(
  "Amazônia"       = "#1a7a4a",
  "Caatinga"       = "#c0392b",
  "Mata Atlântica" = "#2471a3",
  "Pampa"          = "#7d3c98"
)

ind_long <- ind_all |>
  select(date, bioma, wbgt_c, hi_c, thi_c, utci_c, pet_c, et_c) |>
  pivot_longer(
    cols      = c(wbgt_c, hi_c, thi_c, utci_c, pet_c, et_c),
    names_to  = "indice",
    values_to = "valor"
  ) |>
  filter(!is.na(valor)) |>
  mutate(indice = recode(indice,
    wbgt_c = "WBGT (°C)",
    hi_c   = "Heat Index (°C)",
    thi_c  = "THI (°C)",
    utci_c = "UTCI (°C)",
    pet_c  = "PET (°C)",
    et_c   = "ET (°C)"
  ))

ggplot(ind_long, aes(x = valor, y = bioma, fill = bioma)) +
  ggridges::geom_density_ridges(
    alpha = 0.75, scale = 0.9,
    quantile_lines = TRUE, quantiles = c(0.25, 0.5, 0.75),
    quantile_fun = function(x, ...) quantile(x, ...)
  ) +
  scale_fill_manual(values = bioma_colors, guide = "none") +
  facet_wrap(~indice, scales = "free_x", ncol = 3) +
  labs(
    title    = "Distribuição dos Índices Bioclimáticos por Bioma — 2023",
    subtitle = "Linhas internas = quartis (Q1, mediana, Q3)",
    x = "Valor do índice", y = NULL,
    caption  = "Fonte: INMET / climasus4r | Dados horários 2023"
  ) +
  theme(strip.background = element_rect(fill = "#eaf4fb", colour = "grey80"))

Note: the block above uses ggridges; install it with install.packages("ggridges") if needed.

# Versão alternativa sem ggridges
ggplot(ind_long |> filter(!is.na(valor)),
       aes(x = valor, colour = bioma, fill = bioma)) +
  geom_density(alpha = 0.35, linewidth = 0.9) +
  scale_fill_manual(values  = bioma_colors, name = "Bioma") +
  scale_colour_manual(values = bioma_colors, name = "Bioma") +
  facet_wrap(~indice, scales = "free", ncol = 3) +
  labs(
    title   = "Distribuição dos Índices Bioclimáticos — 2023",
    x = "Valor do índice (°C equiv.)", y = "Densidade",
    caption = "Fonte: INMET / climasus4r"
  )

4.2 Seasonal Time Series

# ── Média diária de WBGT por bioma ───────────────────────────────────────────
diario <- ind_all |>
  mutate(data_d = as.Date(date)) |>
  group_by(bioma, data_d) |>
  summarise(
    wbgt_med  = mean(wbgt_c,  na.rm = TRUE),
    utci_med  = mean(utci_c,  na.rm = TRUE),    .groups   = "drop"
  )

# Bandas de risco WBGT (OSHA/ACGIH)
bandas_wbgt <- tibble(
  ymin  = c(-Inf, 20, 25, 28, 30, 33),
  ymax  = c(20,   25, 28, 30, 33, Inf),
  risco = factor(
    c("Nenhum","Baixo","Moderado","Alto","Muito Alto","Extremo"),
    levels = c("Nenhum","Baixo","Moderado","Alto","Muito Alto","Extremo")
  )
)
cores_risco <- c(
  "Nenhum"    = "#d5f5e3",
  "Baixo"     = "#a9dfbf",
  "Moderado"  = "#f9e79f",
  "Alto"      = "#f0b27a",
  "Muito Alto"= "#ec7063",
  "Extremo"   = "#922b21"
)

p1 <- ggplot() +
  geom_rect(
    data = bandas_wbgt,
    aes(xmin = as.Date("2022-12-31"), xmax = as.Date("2024-01-01"),
        ymin = ymin, ymax = ymax, fill = risco),
    alpha = 0.25
  ) +
  scale_fill_manual(values = cores_risco, name = "Risco WBGT") +
  #new_scale_fill() +
  geom_line(data = diario,
            aes(x = data_d, y = wbgt_med, colour = bioma), linewidth = 0.8) +
  scale_colour_manual(values = bioma_colors, name = "Bioma") +
  scale_x_date(date_labels = "%b", date_breaks = "1 month") +
  labs(
    title    = "WBGT Médio Diário — Biomas Brasileiros 2023",
    subtitle = "Bandas de risco: OSHA/ACGIH (adaptadas por bioma)",
    x = NULL, y = "WBGT (°C)",
    caption  = "Zonas de risco: Baixo = 20–25°C | Moderado = 25–28°C | Alto = 28–30°C | Muito Alto = 30–33°C | Extremo > 33°C"
  )

# UTCI e WBGT sobrepostos para Caatinga
p2 <- diario |>
  filter(bioma == "Caatinga") |>
  pivot_longer(c(wbgt_med, utci_med),
               names_to = "var", values_to = "val") |>
  mutate(var = recode(var,
    wbgt_med = "WBGT", utci_med = "UTCI")) |>
  ggplot(aes(data_d, val, colour = var)) +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = c("WBGT"="#c0392b","UTCI"="#8e44ad","T ar"="#2980b9"),
                      name = "Índice") +
  scale_x_date(date_labels = "%b", date_breaks = "2 months") +
  labs(
    title    = "Caatinga — WBGT vs UTCI vs T ar",
    subtitle = "WBGT e UTCI consistentemente acima da temperatura seca",
    x = NULL, y = "°C"
  )

p1 / p2

4.4 Regional Risk Categories

# ── Distribuição mensal das categorias de risco por bioma ────────────────────
risco_mensal <- ind_all |>
  mutate(mes = month(date, label = TRUE, abbr = TRUE, locale = "pt_BR.UTF-8")) |>
  group_by(bioma, mes) |>
  count(heat_stress_risk) |>
  mutate(pct = n / sum(n) * 100) |>
  ungroup() |>
  mutate(
    heat_stress_risk = factor(heat_stress_risk,
      levels = c("None","Low","Moderate","High","Very High","Extreme")),
    risco_label = recode(heat_stress_risk,
      "None"      = "Nenhum",
      "Low"       = "Baixo",
      "Moderate"  = "Moderado",
      "High"      = "Alto",
      "Very High" = "Muito Alto",
      "Extreme"   = "Extremo"
    ),
    risco_label = factor(risco_label,
      levels = c("Nenhum","Baixo","Moderado","Alto","Muito Alto","Extremo"))
  )

ggplot(risco_mensal, aes(mes, pct, fill = risco_label)) +
  geom_col(position = "stack", width = 0.85, colour = "white", linewidth = 0.2) +
  scale_fill_manual(
    values = c(
      "Nenhum"    = "#abebc6",
      "Baixo"     = "#82e0aa",
      "Moderado"  = "#f9e79f",
      "Alto"      = "#f0b27a",
      "Muito Alto"= "#ec7063",
      "Extremo"   = "#7b241c"
    ),
    name = "Risco\ntérmico"
  ) +
  scale_y_continuous(labels = percent_format(scale = 1)) +
  facet_wrap(~bioma, ncol = 2) +
  labs(
    title    = "Composição Mensal do Risco de Estresse Térmico (WBGT-base)",
    subtitle = "Adaptado por bioma: limiares OSHA/ACGIH + correção regional",
    x = NULL, y = "% horas por categoria",
    caption  = paste("Caatinga: pico de risco Extremo em set–out (seca)",
                     "| Amazônia: risco Alto-Muito Alto persistente",
                     "| Pampa: dominado por Nenhum no inverno",
                     sep = "\n")
  ) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

4.5 Agricultural Degree Days

# ── Acumulado de GDD, CDD e HDD por bioma ────────────────────────────────────
grausdia <- ind_all |>
  mutate(data_d = as.Date(date)) |>
  group_by(bioma, data_d) |>
  summarise(
    cdd = sum(cdd_c, na.rm = TRUE),
    hdd = sum(hdd_c, na.rm = TRUE),
    gdd = sum(gdd_c, na.rm = TRUE),
    .groups = "drop"
  ) |>
  group_by(bioma) |>
  arrange(data_d) |>
  mutate(
    cdd_acum = cumsum(cdd),
    hdd_acum = cumsum(hdd),
    gdd_acum = cumsum(gdd)
  ) |>
  ungroup()

p_gdd1 <- grausdia |>
  select(bioma, data_d, cdd_acum, hdd_acum) |>
  pivot_longer(c(cdd_acum, hdd_acum), names_to = "tipo", values_to = "val") |>
  mutate(tipo = recode(tipo,
    cdd_acum = "CDD (resfriamento)",
    hdd_acum = "HDD (aquecimento)"
  )) |>
  ggplot(aes(data_d, val, colour = bioma, linetype = tipo)) +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = bioma_colors, name = "Bioma") +
  scale_linetype_manual(values = c("solid","dashed"), name = NULL) +
  scale_x_date(date_labels = "%b") +
  labs(
    title = "Graus-dia Acumulados — CDD e HDD",
    subtitle = "Base regional (ASHRAE adaptado): 18–20°C para CDD | 14–18°C para HDD",
    x = NULL, y = "Graus-dia acumulados"
  )

p_gdd2 <- grausdia |>
  ggplot(aes(data_d, gdd_acum, colour = bioma)) +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = bioma_colors, name = "Bioma") +
  scale_x_date(date_labels = "%b") +
  labs(
    title = "GDD Acumulados — Potencial Agrícola",
    subtitle = "Base 8–15°C (regional) | Limite superior 30–35°C",
    x = NULL, y = "GDD acumulados"
  )

p_gdd1 | p_gdd2

4.7 Thermal Amplitude and Vapor Pressure

vapor_dtr <- ind_all |>
  mutate(
    mes = month(date),
    trimestre = case_when(
      mes %in% 1:3  ~ "T1 (Jan–Mar)",
      mes %in% 4:6  ~ "T2 (Abr–Jun)",
      mes %in% 7:9  ~ "T3 (Jul–Set)",
      TRUE          ~ "T4 (Out–Dez)"
    )
  ) |>
  filter(!is.na(diurnal_range_c), !is.na(vapor_pressure_kpa))

p_vap <- ggplot(vapor_dtr,
       aes(diurnal_range_c, vapor_pressure_kpa, colour = bioma)) +
  geom_point(alpha = 0.08, size = 0.6) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 1.1) +
  scale_colour_manual(values = bioma_colors, name = "Bioma") +
  scale_y_continuous(labels = label_number(suffix = " kPa")) +
  labs(
    title    = "Amplitude Térmica Diurna × Pressão de Vapor Atual",
    subtitle = "Caatinga: alta amplitude + baixa pressão de vapor (aridez)\nAmazônia: baixa amplitude + alta pressão (carga evaporativa máxima)",
    x = "DTR — Amplitude Diurna (°C)", y = "Pressão de vapor ea (kPa)"
  ) +
  facet_wrap(~bioma, ncol = 2) +
  theme(legend.position = "none")

p_vap

4.8 Categorical Classifications

# ── Köppen Humidity e Heat Stress Risk — proporção anual por bioma ────────────

koppen_prop <- ind_all |>
  count(bioma, koppen_humidity) |>
  group_by(bioma) |>
  mutate(pct = n / sum(n) * 100) |>
  ungroup()

risk_anual <- ind_all |>
  count(bioma, heat_stress_risk) |>
  group_by(bioma) |>
  mutate(
    pct = n / sum(n) * 100,
    heat_stress_risk = factor(heat_stress_risk,
      levels = c("None","Low","Moderate","High","Very High","Extreme"))
  ) |>
  ungroup()

p_kopp <- ggplot(koppen_prop, aes(bioma, pct, fill = koppen_humidity)) +
  geom_col(width = 0.7, colour = "white", linewidth = 0.3) +
  scale_fill_manual(
    values = c("Arid"="#f0b27a","Semi-arid"="#f9e79f",
               "Humid"="#82e0aa","Perhumid"="#2ecc71"),
    name = "Classe Köppen"
  ) +
  scale_y_continuous(labels = percent_format(scale = 1)) +
  coord_flip() +
  labs(
    title = "Classificação Köppen de Umidade",
    subtitle = "% horas anuais por classe",
    x = NULL, y = "% horas"
  )

p_risk <- ggplot(risk_anual, aes(bioma, pct, fill = heat_stress_risk)) +
  geom_col(width = 0.7, colour = "white", linewidth = 0.3) +
  scale_fill_manual(
    values = c(
      "None"      = "#abebc6",
      "Low"       = "#82e0aa",
      "Moderate"  = "#f9e79f",
      "High"      = "#f0b27a",
      "Very High" = "#ec7063",
      "Extreme"   = "#7b241c"
    ),
    name = "Risco Térmico\n(WBGT-base)"
  ) +
  scale_y_continuous(labels = percent_format(scale = 1)) +
  coord_flip() +
  labs(
    title = "Risco de Estresse Térmico (OSHA/ACGIH)",
    subtitle = "% horas anuais por nível de risco",
    x = NULL, y = "% horas"
  )

p_kopp | p_risk

Effect of the Regional Threshold on HI

# ── Demonstração do efeito do limiar regional no HI ──────────────────────────
t_seq  <- seq(22, 38, 0.2)
rh_val <- 65

calc_hi_vec <- function(airT, rh, min_temp) {
  tf <- (airT * 9/5) + 32
  hi_f <- -42.379 + 2.04901523*tf + 10.14333127*rh - 0.22475541*tf*rh -
           6.83783e-3*tf^2 - 5.481717e-2*rh^2 + 1.22874e-3*tf^2*rh +
           8.5282e-4*tf*rh^2 - 1.99e-6*tf^2*rh^2
  hi_c <- (hi_f - 32) * 5 / 9
  ifelse(airT < min_temp, NA, hi_c)
}

df_hi_reg <- tibble(
  airT = rep(t_seq, 3),
  bioma = rep(c("Sul/Pampa (26°C)","Padrão (26.7°C)","Amazônia (24°C)"), each = length(t_seq)),
  hi    = c(
    calc_hi_vec(t_seq, rh_val, 26.0),
    calc_hi_vec(t_seq, rh_val, 26.7),
    calc_hi_vec(t_seq, rh_val, 24.0)
  )
)

ggplot(df_hi_reg, aes(airT, hi, colour = bioma)) +
  geom_line(linewidth = 1.3, na.rm = TRUE) +
  geom_vline(xintercept = c(24.0, 26.0, 26.7),
             linetype = "dotted", colour = c("#1a7a4a","#7d3c98","#2471a3"),
             linewidth = 0.8) +
  annotate("text", x = c(24.2, 26.2, 26.9), y = 22,
           label = c("Amazônia\n24°C","Sul\n26°C","Padrão\n26.7°C"),
           hjust = 0, size = 3, colour = c("#1a7a4a","#7d3c98","#2471a3")) +
  scale_colour_manual(
    values = c("Sul/Pampa (26°C)"="#7d3c98",
               "Padrão (26.7°C)"="#2471a3",
               "Amazônia (24°C)"="#1a7a4a"),
    name = "Limiar regional"
  ) +
  labs(
    title    = "Impacto do Limiar Regional no Heat Index (RH = 65%)",
    subtitle = "Populações tropicais aclimatadas percebem desconforto a T mais baixas",
    x = "Temperatura do ar (°C)", y = "Heat Index (°C)",
    caption  = "Rothfusz (1990) — limiar padrão 26.7°C | adaptado: 24–26°C para trópicos"
  )

Uncertainty & Confidence Flags

Uncertainty Estimation (Monte Carlo)

# ── Incerteza para São Paulo — 200 simulações ─────────────────────────────────
ind_sp_unc <- df_sp_filled |>
  sus_climate_compute_indicators(
    indicators          = c("wbgt", "hi", "utci"),
    region              = "atlantic_forest",
    compute_uncertainty = TRUE,      # ← 200 simulações MC
    confidence_flags    = TRUE,
    verbose             = TRUE
  )

# Colunas geradas: wbgt_c, wbgt_c_ci_low, wbgt_c_ci_high
# Interpretação: IC 95% propagando incertezas T±0.5°C, RH±3%, WS±0.3 m/s, SR±10%

ind_sp_unc |>
  filter(!is.na(wbgt_c_ci_low)) |>
  mutate(ci_width = wbgt_c_ci_high - wbgt_c_ci_low) |>
  summary()

Confidence Flags

# ── Visualizar os flags gerados ───────────────────────────────────────────────
flags_summary <- ind_all |>
  select(bioma, date, wbgt_c,
         wbgt_c_flag_extreme, wbgt_c_flag_high, wbgt_c_flag_low) |>
  group_by(bioma) |>
  summarise(
    pct_extreme = 100 * mean(wbgt_c_flag_extreme, na.rm = TRUE),
    pct_high    = 100 * mean(wbgt_c_flag_high,    na.rm = TRUE),
    pct_low     = 100 * mean(wbgt_c_flag_low,     na.rm = TRUE),
    .groups     = "drop"
  ) |>
  pivot_longer(-bioma, names_to = "flag", values_to = "pct") |>
  mutate(
    flag = recode(flag,
      pct_extreme = "Flag: Extremo (> 31°C)",
      pct_high    = "Flag: Alto (> 28°C)",
      pct_low     = "Flag: Frio (< 15°C)"
    )
  )

ggplot(flags_summary, aes(bioma, pct, fill = flag)) +
  geom_col(position = "dodge", width = 0.7,
           colour = "white", linewidth = 0.3) +
  scale_fill_manual(
    values = c("Flag: Extremo (> 31°C)"="#922b21",
               "Flag: Alto (> 28°C)"   ="#e67e22",
               "Flag: Frio (< 15°C)"  ="#2980b9"),
    name = NULL
  ) +
  scale_y_continuous(labels = percent_format(scale = 1)) +
  labs(
    title    = "Frequência dos Flags de Confiança WBGT por Bioma",
    subtitle = "Flags: _flag_extreme | _flag_high | _flag_low — gerados automaticamente",
    x = NULL, y = "% horas com flag ativo",
    caption  = paste("Limiares adaptados por bioma:",
                     "Amazônia: extremo>31°C | Caatinga: extremo>33°C | Sul: frio<15°C",
                     sep = "\n")
  )

Advanced Usage Patterns

Full Integration into the SUS Pipeline

# ── Pipeline de produção para análise epidemiológica ─────────────────────────

library(climasus4r)

# 1. Download multi-UF, multi-ano
df_raw <- sus_climate_inmet(
  years    = 2020:2023,
  uf       = c("AM","PA","CE","BA","SP","MG","RS"),
  parallel = TRUE,
  workers  = 6
)

# 2. Preenchimento XGBoost
df_filled <- df_raw |>
  sus_climate_fill_inmet(
    target_var        = "all",
    quality_threshold = 0.40,
    parallel          = TRUE,
    workers           = 6
  )

# 3. Indicadores com adaptação regional automática
df_indicators <- df_filled |>
  sus_climate_compute_indicators(
    region              = "auto",       # detecta bioma por UF/latitude
    confidence_flags    = TRUE,
    verify_physics      = TRUE,
    compute_uncertainty = FALSE,        # TRUE para análise de sensibilidade
    keep_source_vars    = FALSE
  )

# 4. Resumo mensal por estação — para linkage com dados de mortalidade SIM
resumo_mensal <- df_indicators |>
  mutate(
    ano_mes = format(date, "%Y-%m"),
    station = station_code
  ) |>
  group_by(station, ano_mes) |>
  summarise(
    across(c(wbgt_c, hi_c, utci_c, pet_c, thi_c),
           list(med = ~mean(.x, na.rm=TRUE),
                p90 = ~quantile(.x, 0.90, na.rm=TRUE),
                p95 = ~quantile(.x, 0.95, na.rm=TRUE)),
           .names = "{.col}_{.fn}"),
    horas_risco_alto    = sum(heat_stress_risk %in% c("High","Very High","Extreme"),
                              na.rm = TRUE),
    horas_risco_extremo = sum(heat_stress_risk == "Extreme", na.rm = TRUE),
    cdd_acum            = sum(cdd_c, na.rm = TRUE),
    gdd_acum            = sum(gdd_c, na.rm = TRUE),
    .groups             = "drop"
  )

Custom Thresholds for Research

# ── Thresholds customizados para estudo específico ───────────────────────────
ind_custom <- df_am_filled |>
  sus_climate_compute_indicators(
    region = "amazon",
    custom_thresholds = list(
      hi_min_temp = 22.0,   # limiar ainda mais conservador
      hi_min_rh   = 35,
      cdd_base    = 22,     # base climática local calibrada
      gdd_base    = 18,
      gdd_upper   = 33
    ),
    confidence_flags = TRUE
  )

Selected Indicators Only

# ── Cálculo seletivo — mais eficiente para análises focadas ──────────────────

# Só indicadores de saúde pública ocupacional
ind_ocupacional <- df_ce_filled |>
  sus_climate_compute_indicators(
    indicators       = c("wbgt", "utci", "heat_stress_risk"),
    region           = "caatinga",
    confidence_flags = TRUE,
    verbose          = FALSE
  )

# Só graus-dia para modelagem energética/agrícola
ind_energia <- df_rs_filled |>
  sus_climate_compute_indicators(
    indicators = c("cdd", "hdd", "gdd", "diurnal_range"),
    region     = "pampa",
    verbose    = FALSE
  )

# Só indicadores que não dependem de radiação solar (dados incompletos)
ind_sem_sr <- df_sp_filled |>
  sus_climate_compute_indicators(
    indicators = c("hi", "thi", "et", "vapor_pressure", "koppen_humidity"),
    region     = "atlantic_forest",
    verbose    = FALSE
  )

cat("Colunas geradas (ocupacional):", names(ind_ocupacional), "\n")

Risk Category Reference Table

🌡️ Thermal Risk Classification by Bioclimatic Indices

Index	Category	Threshold (°C)	Physiological Meaning	Source
WBGT	None	< 20	No perceptible thermal risk	OSHA/ACGIH
WBGT	Low	20–25	Mild stress during intense activity	OSHA/ACGIH
WBGT	Moderate	25–28	Rest breaks recommended for workers	OSHA/ACGIH
WBGT	High	28–30	Real risk of thermal disorders	OSHA/ACGIH
WBGT	Very High	30–33	Severe reduction in working time	OSHA/ACGIH
WBGT	Extreme	> 33	Risk of heat collapse	OSHA/ACGIH
UTCI	Extreme cold	< −27	Imminent hypothermia	Bröde et al. (2012)
UTCI	Neutral	9–26	Thermal comfort zone	Bröde et al. (2012)
UTCI	Mild heat	26–32	Mild heat stress	Bröde et al. (2012)
UTCI	Strong heat	32–38	Moderate heat stress	Bröde et al. (2012)
UTCI	Extreme heat	> 46	Extreme thermal stress	Bröde et al. (2012)
PET	Extreme cold	< 4	Thermal discomfort from cold	Matzarakis et al. (1999)
PET	Comfort	18–23	Neutral thermal zone	Matzarakis et al. (1999)
PET	Mild heat	23–29	Mild heat discomfort	Matzarakis et al. (1999)
PET	Extreme heat	> 41	Very strong thermal stress	Matzarakis et al. (1999)
HI	Caution	27–32	Fatigue possible with prolonged exposure	NWS
HI	Extreme danger	> 54	Heat stroke highly likely	NWS

Scientific Considerations and Limitations

⚠️ Important Limitations to Report

UTCI and PET are approximations — the full UTCI (Fiala) requires a polynomial with 210+ terms; this implementation may diverge by 3–5°C under extreme conditions (Bröde et al., 2012). For critical epidemiological studies, use reticulate + pythermalcomfort for the full UTCI.
Nighttime solar radiation = 0 — NA values in sr_kj_m2 during the night are imputed as 0 (no solar load). Nighttime WBGT reflects only temperature, humidity, and wind — physically correct.
PET and clothing — the simplified model assumes average seasonal clothing (clo 0.5–0.9); populations with very different habits (e.g., miners, rural workers) require specific calibration.
Regional adaptations are empirical — the offsets (+1.5°C Amazon, −0.5°C Caatinga, etc.) reflect bibliographic consensus; local validation against globe thermometer measurements is recommended for publications.
INMET data — automatic stations have maintenance issues; XGBoost gap-filling reduces gaps but does not eliminate measurement uncertainty. Use compute_uncertainty = TRUE to propagate this uncertainty.

References

Alduchov, O. A., & Eskridge, R. E. (1996). Improved Magnus form approximation of saturation vapor pressure. Journal of Applied Meteorology, 35, 601–609.
Bernard, T. E., & Pourmoghani, M. (1999). Prediction of workplace wet bulb globe temperature. AIHAJ, 60(1), 32–37.
Bröde, P., Fiala, D., Błażejczyk, K., Holmér, I., Jendritzky, G., Kampmann, B., et al. (2012). Deriving the operational procedure for the Universal Thermal Climate Index (UTCI). International Journal of Biometeorology, 56(4), 481–494.
International Organization for Standardization. (1989). ISO 7243: Hot environments—Estimation of the heat stress on working man based on the WBGT-index.
Joint Action Group for Temperature Indices. (2001). Development of the wind chill temperature index.
Liljegren, J. C., Carhart, R. A., Lawday, P., Tschopp, S., & Sharp, R. (2008). Modeling the wet bulb globe temperature using standard meteorological measurements. Journal of Occupational and Environmental Hygiene, 5(10), 645–655.
Magnus, G. (1844). Versuche über die spannkräfte des wasserdampfs.
Matzarakis, A., Mayer, H., & Iziomon, M. G. (1999). Applications of a universal thermal index: Physiological equivalent temperature. International Journal of Biometeorology, 43, 76–84.
Missenard, A. (1937). Température sensible.
McMaster, G. S., & Wilhelm, W. W. (1997). Growing degree-days: One equation, two interpretations. Agricultural and Forest Meteorology, 87, 291–300.
Missenard, A. (1937). Température sensible.
National Weather Service. (n.d.). Heat index equation. https://www.wpc.ncep.noaa.gov/html/heatindex_equation.shtml
Occupational Safety and Health Administration. (2017). Protecting workers from the effects of heat (OSHA 3154).
Rothfusz, L. P. (1990). The heat index equation (Technical Attachment SR/SSD 90-23). National Weather Service.
Tetens, O. (1930). Über einige meteorologische Begriffe. Zeitschrift für Geophysik, 6, 297–309.
Thom, E. C. (1959). The discomfort index. Weatherwise, 12(2), 57–59.

🌲 Have feedback or suggestions?

Do you have an idea for improvement or found an error? Click the button below to open a new issue on GitHub and share your suggestions directly with us.

Open issue on GitHub

Theory and INMET → Gap-filling → Thermal Indicators Pipeline | `sus_climate_compute_indicators()`

Max Anjos

João Gobo

Thauã Menezes

17 October, 2026

Theoretical Foundations

Why bioclimatic indices?

📊 Bioclimatic and Thermal Stress Indicators

Mathematical Formulas

1. WBGT (Wet-Bulb Globe Temperature)

2. Heat Index (HI)

3. Universal Thermal Climate Index (UTCI)

4. Degree Days (CDD, HDD, GDD)

Full Pipeline

Step 1 · INMET Import

Step 2 · Gap-Filling

Step 3 · Computing the Indicators

Visualizations

4.1 Overview of the Indicators

4.2 Seasonal Time Series

4.4 Regional Risk Categories

4.5 Agricultural Degree Days

4.7 Thermal Amplitude and Vapor Pressure

4.8 Categorical Classifications

Effect of the Regional Threshold on HI

Uncertainty & Confidence Flags

Uncertainty Estimation (Monte Carlo)

Confidence Flags

Advanced Usage Patterns

Full Integration into the SUS Pipeline

Custom Thresholds for Research

Selected Indicators Only

Risk Category Reference Table

🌡️ Thermal Risk Classification by Bioclimatic Indices

Scientific Considerations and Limitations

References

Thom, E. C. (1959). The discomfort index. Weatherwise, 12(2), 57–59.

🌲 Have feedback or suggestions?

Bioclimatic Indices for Public Health

Theory and INMET → Gap-filling → Thermal Indicators Pipeline | sus_climate_compute_indicators()

Max Anjos

João Gobo

Thauã Menezes

17 October, 2026

Theoretical Foundations

Why bioclimatic indices?

📊 Bioclimatic and Thermal Stress Indicators

Mathematical Formulas

1. WBGT (Wet-Bulb Globe Temperature)

2. Heat Index (HI)

3. Universal Thermal Climate Index (UTCI)

4. Degree Days (CDD, HDD, GDD)

Full Pipeline

Step 1 · INMET Import

Step 2 · Gap-Filling

Step 3 · Computing the Indicators

Visualizations

4.1 Overview of the Indicators

4.2 Seasonal Time Series

4.4 Regional Risk Categories

4.5 Agricultural Degree Days

4.7 Thermal Amplitude and Vapor Pressure

4.8 Categorical Classifications

Effect of the Regional Threshold on HI

Uncertainty & Confidence Flags

Uncertainty Estimation (Monte Carlo)

Confidence Flags

Advanced Usage Patterns

Full Integration into the SUS Pipeline

Custom Thresholds for Research

Selected Indicators Only

Risk Category Reference Table

🌡️ Thermal Risk Classification by Bioclimatic Indices

Scientific Considerations and Limitations

References

Thom, E. C. (1959). The discomfort index. Weatherwise, 12(2), 57–59.

🌲 Have feedback or suggestions?

Theory and INMET → Gap-filling → Thermal Indicators Pipeline | `sus_climate_compute_indicators()`