Skip to contents

Objective

By the end of this tutorial you will be able to bring the observed climate dimension into the pipeline using data from the National Institute of Meteorology (INMET) weather stations. Specifically, you will learn to:

  1. Import the hourly series of a station with sus_climate_inmet();
  2. Fill gaps in the series with machine learning models (XGBoost) per station, with sus_climate_fill_inmet();
  3. Visualize the gap-filling result (observed vs imputed) with sus_climate_plot_fill();
  4. Compute the 30-year climatological normals with sus_climate_normals(), which serve as a climatological baseline;
  5. Align the climate series with the health series using epidemiologically grounded temporal strategies (same day, lag, moving window…) with sus_climate_aggregate().

This is the seventh stage of the pipeline (stage = "climate"). It connects the already-aggregated health series (Tutorial 4) to observed climate, producing the exposure × outcome base that will feed the gridded climate stage (Tutorial 8) and modelling (Tutorial 9).

Scientific context

Environmental epidemiology relates climate exposures to health outcomes. In Brazil, INMET’s automatic stations provide hourly observations of temperature, humidity, precipitation, pressure, and wind, with growing national coverage since the 2000s. Unlike gridded reanalysis products (Tutorial 8), station data measure the environment in situ, with high temporal resolution but with limited spatial representativeness around the station (Saldanha et al., 2019; Geirinhas et al., 2018).

Three methodological challenges run through this stage:

  • Quality and gaps. Sensors fail; series have missing values and physical inconsistencies. Gap-filling with models such as XGBoost (Chen & Guestrin, 2016) recovers gaps from temporal patterns (lags, daily and seasonal cycles), preserving the series structure without introducing the crude biases of linear interpolation.
  • Climatological baseline. The climatological normals (30-year averages) define the reference climate of each location — essential for contextualizing anomalies and defining the reference point of dose-exposure response models (Gasparrini, 2011).
  • Causal temporal alignment. The effect of heat or cold on health is rarely instantaneous: there are lags (delayed effects) and accumulation (cumulative exposure). The choice of temporal strategy must reflect the hypothesized biological mechanism (Bhaskaran et al., 2013; Gasparrini et al., 2010), and must never use future climate to explain a present event (anticipation bias).

Guiding case. We investigate temperature and health outcomes in older adults (60+) in the São Paulo Metropolitan Region (RMSP), 2014–2019. In this stage, we import the series of the RMSP reference station — São Paulo – Mirante de Santana (code A701) —, fill its gaps, compute the normals, and align temperature with the daily health series inherited from Tutorial 4. The resulting base (health × climate) is saved for Tutorial 8.

Data used

Source Function Role in the case
Hourly INMET series, station A701 (RMSP) sus_climate_inmet() Observed climate exposure
A701 series with gaps filled (XGBoost) sus_climate_fill_inmet() Continuous temperature series
Climatological normals 1991–2020 (SP) sus_climate_normals() Climatological baseline
Daily health series from T4 (dados/caso_serie.rds) sus_climate_aggregate() Outcome to be aligned

Chained input: dados/caso_serie.rds (health series from Tutorial 4). Chained output: dados/caso_clima_estacao.rds (health + aligned climate), which Tutorial 8 will load. All data are public and real (INMET).

Full pipeline

The blocks below reproduce what the script scripts/tutorial_07_clima_estacoes.R executes. Load the package:

library(climasus4r)

Import the station series (A701)

sus_climate_inmet() downloads INMET’s hourly data, standardizes the column names to canonical names (e.g., tair_dry_bulb_c, rh_mean_porc, rainfall_mm) and applies physical quality control. We filter by the RMSP reference station via station_code.

clima_horario <- sus_climate_inmet(
  years        = 2014:2019,
  uf           = "SP",
  station_code = "A701",    # São Paulo - Mirante de Santana
  use_cache    = TRUE,       # cache em dois níveis (memória + disco/Parquet)
  parallel     = TRUE,       # paraleliza entre anos e entre arquivos do ano
  workers      = 4,
  lang         = "pt",
  verbose      = TRUE
)

The return is a climasus_df at stage = "climate", type = "inmet", with hourly frequency and dates in UTC.

Fill gaps with XGBoost

sus_climate_fill_inmet() fits an XGBoost model per station and per variable, automatically creating temporal features (hour, day, month, cyclical transforms), lags (1, 2, 3, 12, 24, 168 periods) and rolling statistics. We first assess accuracy (evaluation mode), then actually fill the gaps (production mode):

# Modo AVALIAÇÃO: cria lacunas artificiais (MCAR) e devolve métricas por estação
aval <- sus_climate_fill_inmet(
  df                = clima_horario,
  target_var        = "tair_dry_bulb_c",
  quality_threshold = 0.4,    # exclui estações com > 40% de faltantes
  run_evaluation    = TRUE,
  gap_percentage    = 0.2,    # 20% de lacunas artificiais
  parallel          = TRUE,
  workers           = 4,
  lang              = "pt"
)
aval$tair_dry_bulb_c$metrics   # RMSE, MAE, R², sMAPE, slope_bias, n_gaps
# Modo PRODUÇÃO: preenche as lacunas reais e adiciona a flag is_imputed_<var>
clima_preenchido <- sus_climate_fill_inmet(
  df                = clima_horario,
  target_var        = "tair_dry_bulb_c",
  quality_threshold = 0.4,
  run_evaluation    = FALSE,
  keep_features     = FALSE,   # retorna só as colunas originais + flag
  parallel          = TRUE,
  workers           = 4,
  lang              = "pt"
)

To impute all 17 canonical INMET variables at once, use target_var = "all".

Visualize the gap-filling

sus_climate_plot_fill() automatically detects the mode (production vs evaluation) and draws the observed series with the imputed points highlighted. In static mode (interactive = FALSE) it returns a ggplot2 object:

sus_climate_plot_fill(
  df_filled   = clima_preenchido,
  df_original = clima_horario,
  target_var  = "tair_dry_bulb_c",
  interactive = FALSE,        # FALSE = ggplot2; TRUE = plotly interativo
  output_type = "plot",       # "plot", "table", "metrics" ou "all"
  lang        = "pt"
)

Compute the climatological normals

sus_climate_normals() downloads INMET’s 30-year normals. We choose the most recent period and the temperature variables:

normais <- sus_climate_normals(
  period     = "1991-2020",                       # ou "1961-1990", "1981-2010"
  target_var = c("t_max", "t_min", "t_mean_comp"),
  use_cache  = TRUE,
  lang       = "pt"
)

# Para explorar os códigos de variável disponíveis antes:
# sus_climate_normals_meta(period = "1991-2020", lang = "pt")

Align climate and health (temporal strategies)

sus_climate_aggregate() links each health record to the nearest station (Euclidean distance) and applies the chosen temporal strategy. We compare same-day exposure (exact) with the 7-day moving window (moving_window, recommended for cumulative exposure):

# Exposição aguda (mesmo dia)
saude_exact <- sus_climate_aggregate(
  health_data       = saude_df,             # série de saúde do T4 (com geometria)
  climate_data      = clima_preenchido,
  climate_var       = "tair_dry_bulb_c",
  temporal_strategy = "exact",
  lang              = "pt"
)

# Exposição cumulativa (média móvel retroativa de 7 dias: t-6 .. t)
saude_movel <- sus_climate_aggregate(
  health_data       = saude_df,
  climate_data      = clima_preenchido,
  climate_var       = "tair_dry_bulb_c",
  temporal_strategy = "moving_window",
  window_days       = 7,
  lang              = "pt"
)

The resulting series (health + aligned temperature) is saved to dados/caso_clima_estacao.rds for Tutorial 8.

Function characteristics

sus_climate_inmet()

Imports, standardizes, and quality-controls INMET’s hourly data.

Argument Essential? Description
years recommended Year(s), 2000–current year (e.g., 2014:2019); NULL uses the last 2
uf recommended State(s) ("SP", c("RJ","MG")); NULL downloads all 27 states
station_code optional Filters stations by code ("A701"); NULL = all in the state
use_cache / cache_dir recommended Two-level cache (memory + Parquet/Zstd)
parallel / workers optional Parallelizes across years and across CSVs within a year
lang / verbose optional Language ("pt"/"en"/"es") and progress

Return: climasus_df (stage = "climate", type = "inmet"), hourly frequency, dates in UTC, canonical columns. Automatic QC: validation of physical ranges, dew point consistency (Magnus formula), and zeroed nighttime solar radiation.

sus_climate_fill_inmet()

Imputes missing values with XGBoost, one model per station per variable.

Argument Essential? Description
df yes Output of sus_climate_inmet()
target_var yes Column(s) to impute, or "all" (the 17 canonical variables)
quality_threshold recommended Max. missing proportion per station (default 0.4)
run_evaluation optional TRUE creates artificial gaps and returns metrics
gap_percentage optional Proportion of artificial gaps in evaluation (default 0.2)
keep_features optional TRUE keeps engineered features (lags/rolling means)
parallel / workers optional Processes stations in parallel (furrr)
datetime_col / station_col optional Column names; auto-detected if NULL

Return (production): the same tibble with target_var imputed and flag is_imputed_<var>. Return (evaluation): named list with $data and $metrics (RMSE, MAE, R², sMAPE, slope_bias, n_gaps) per station.

sus_climate_plot_fill()

Visualizes the gap-filling; automatically detects production vs evaluation.

Argument Essential? Description
df_filled yes Output of sus_climate_fill_inmet()
df_original production only Raw series (from sus_climate_inmet())
target_var yes Variable(s) to visualize
interactive optional FALSE → ggplot2; TRUE → plotly
output_type optional "plot", "table", "metrics" or "all"
save_plot optional Path to save (.png or .html)

Return: the plot, the table, the metrics, or a list with all four components. Evaluation mode draws a 2×2 diagnostic panel (observed vs predicted, residual histogram, residuals over time and vs observed).

sus_climate_normals()

Downloads INMET’s 30-year climatological normals.

Argument Essential? Description
period recommended "1961-1990", "1981-2010" or "1991-2020" (default)
target_var optional Variable codes (c("t_max","t_min")); NULL = all
use_cache / cache_dir recommended Disk cache of the catalog and files
lang / verbose optional Language and progress

Return: climasus_df (type = "normals"), long format with codigo, nome_estacao, uf, mes, decada, valor, var_code and multilingual labels. Companion function: sus_climate_normals_meta() lists the available codes.

sus_climate_aggregate()

Links climate to health via the nearest station and applies the temporal strategy.

Argument Essential? Description
health_data yes climasus_df with date, code_muni and geometry (sf)
climate_data yes Output of sus_climate_fill_inmet()
climate_var recommended Climate variable(s); "all" (default)
temporal_strategy yes One of 10 strategies (see below)
window_days conditional Window size (moving, degree-days, threshold, weighted)
lag_days conditional Lag(s) for discrete_lag / max for distributed_lag
offset_days conditional c(W1, W2) for offset_window
temp_base / gdd_temp_var conditional Thermal base and variable for degree_days
threshold_value / threshold_direction conditional Threshold for *_exceedance
climate_region optional "auto" (default), "tropical", "subtropical", "temperate"
min_obs optional Minimum proportion of valid observations in the window (default 0.7)

Temporal strategies (temporal_strategy): exact, discrete_lag, moving_window, offset_window, distributed_lag, degree_days, seasonal, threshold_exceedance, cold_wave_exceedance, weighted_window.

Return: climasus_df with the health data and the integrated climate variables as new columns.

Formulas

Backward-looking moving average (moving_window, window of WW days). The cumulative exposure of day tt is the average of the observations from t−W+1t-W+1 to tt — it never uses the future (avoids anticipation bias):

x‾t(W)=1W∑ℓ=0W−1xt−ℓ \bar{x}_t^{(W)} = \frac{1}{W} \sum_{\ell=0}^{W-1} x_{t-\ell}

Discrete lag (discrete_lag). Exposure is the value observed exactly LL days before the event:

xt(L)=xt−L x_t^{(L)} = x_{t-L}

Degree-days (degree_days, base TbT_b). Accumulation of the temperature excess over the base across the window (thermal stress):

GDDt=∑ℓ=0W−1max⁡(xt−ℓ−Tb,0) \mathrm{GDD}_t = \sum_{\ell=0}^{W-1} \max\!\left(x_{t-\ell} - T_b,\; 0\right)

Threshold exceedance (threshold_exceedance). Number of days in the window in which temperature exceeds the threshold τ\tau (heat wave proxy):

Nt=∑ℓ=0W−1𝟙(xt−ℓ>τ) N_t = \sum_{\ell=0}^{W-1} \mathbb{1}\!\left(x_{t-\ell} > \tau\right)

Illustrated results

Artefato indisponível — rode Rscript scripts/tutorial_07_clima_estacoes.R (a partir da raiz do pacote) para gerá-lo.

Figure 1. Daily mean temperature observed at station A701 (São Paulo – Mirante de Santana), 2014–2019. The winter seasonality (lower temperatures between June and August) is clear and will be the axis of alignment with mortality in older adults.

Artefato indisponível — rode Rscript scripts/tutorial_07_clima_estacoes.R (a partir da raiz do pacote) para gerá-lo.

Artefato indisponível — rode Rscript scripts/tutorial_07_clima_estacoes.R (a partir da raiz do pacote) para gerá-lo.

Figure 2. Result of the gap-filling (sus_climate_plot_fill()): the observed series with the values imputed by XGBoost highlighted. The imputation preserves the daily and seasonal cycles, rather than “flattening” the missing stretches.

Artefato indisponível — rode Rscript scripts/tutorial_07_clima_estacoes.R (a partir da raiz do pacote) para gerá-lo.

Figura 3 indisponível: as normais climatológicas exigem acesso ao portal do INMET. Reexecute o script de geração com internet.

Figure 3. Annual cycle of the temperature climatological normals (1991–2020) for the state of São Paulo. The normals provide the climatological baseline against which anomalies and extremes are measured in the following tutorials.

Artefato indisponível — rode Rscript scripts/tutorial_07_clima_estacoes.R (a partir da raiz do pacote) para gerá-lo.

Figure 4. Temperature exposure by temporal strategy in the winter of 2014 (RMSP): the same-day value (exact) oscillates more, while the 7-day moving average (moving_window) smooths the signal and captures cumulative exposure — a choice that reflects the biological mechanism of interest in the final model.

References

  • Bhaskaran, K., Gasparrini, A., Hajat, S., Smeeth, L., & Armstrong, B. (2013). Time series regression studies in environmental epidemiology. International Journal of Epidemiology, 42(4), 1187–1195.
  • Chen, T., & Guestrin, C. (2016). XGBoost: A scalable tree boosting system. Proceedings of the 22nd ACM SIGKDD, 785–794.
  • Gasparrini, A. (2011). Distributed lag linear and non-linear models in R: the package dlnm. Journal of Statistical Software, 43(8), 1–20.
  • Gasparrini, A., Armstrong, B., & Kenward, M. G. (2010). Distributed lag non-linear models. Statistics in Medicine, 29(21), 2224–2234.
  • Geirinhas, J. L., et al. (2018). Climatic and synoptic characterization of heat waves in Brazil. International Journal of Climatology, 38(4), 1760–1776.
  • Saldanha, R. de F., Bastos, R. R., & Barcellos, C. (2019). Microdatasus: pacote para download e pré-processamento de microdados do DATASUS. Cadernos de Saúde Pública, 35(9), e00032419.
  • INMET — Instituto Nacional de Meteorologia. Dados históricos. Available at: https://portal.inmet.gov.br/dadoshistoricos.