Predicción de anomalías del balance hídrico mediante Random Forest

Un análisis basado en datos climáticos mensuales (1995–2024) para Madrid - España

Dasaset: Capital Weather Data - 1995-2024

Referencia: Predicting agricultural drought by Endre Harsányi, 2025

1. INTRODUCCIÓN

La sequía es un fenómeno climático complejo que resulta de la interacción entre la precipitación y la demanda atmosférica de agua, afectando directamente la disponibilidad hídrica y la productividad de los cultivos. Su análisis se apoya comúnmente en indicadores como el balance hídrico climático (D = P − PET) y en índices estandarizados de sequía, entre ellos el Índice de Precipitación–Evapotranspiración Estandarizado (SPEI), ampliamente utilizado para caracterizar condiciones secas y húmedas a diferentes escalas temporales.

Sin embargo, la modelación predictiva directa de índices estandarizados como el SPEI presenta limitaciones, especialmente en series temporales de extensión moderada, debido a su normalización y a la alta variabilidad climática inherente. En este estudio se utiliza como variable objetivo la anomalía del balance hídrico mensual (D_anom), calculada a partir del balance hídrico y su climatología mensual, con el fin de reducir la estacionalidad y resaltar condiciones anómalas de déficit o exceso hídrico.

La modelación se realiza mediante un enfoque de regresión Random Forest, implementado con el framework tidymodels, utilizando como variables predictoras la precipitación, temperatura mínima, media, máxima, el balance hídrico mensual, la energía electromagnética, PET y el índice SPEI-3.

El objetivo es evaluar la capacidad del modelo para predecir anomalías del balance hídrico y explorar su utilidad para el análisis y monitoreo de la sequía en regiones mediterráneas.

A partir de este planteamiento, se formula la siguiente pregunta de investigación:

• ¿Es posible modelar condiciones anómalas del balance hídrico a partir algoritmos de aprendizaje automático y variables climáticas?

2. Descripción del dataset y sus variables

Para este estudio, utilizamos el conjunto de datos “Capital Weather Data 1995-2024” (obtenido de Kaggle), que proporciona una serie temporal diaria de variables meteorológicas. Este dataset servirá como base para construir las variables climáticas predictoras necesarias para evaluar el estudio planteado.

El dataset abarca un periodo histórico de casi 30 años (1995-2024), lo cual es ideal para análisis climáticos que requieren series temporales largas para identificar tendencias significativas. La estructura de los datos se alinea con los requisitos del estudio de Harsányi (2025), permitiendo el cálculo del Balance Hídrico Climático (D) y las anomalías asociadas.

El dataset tiene 20 variables de las cuales vamos a usar:

• Date: Fecha de registro de al variable meteorológica.
• Capital: Lugar donde se registro la variable meteorológica.
• Temp_min_c: Hace referencia a la temperatura mínima en grados centígrados.
• Temp_max_x: Hace referencia a la temperatura máxima en grados centígrados.
• Temp_mean_approx: Hace referencia a la tempeeratura media en grados centígrados.
• Precip_mm: Hace referencia a la precipitación en milímetros.
• Shortwave_radiation_MJ_m2: Hace referencia a la energía irradiada en Milijoules por metro cuadrado.

2.1. Variables calculadas

Basándose en la metodología del estudio de Endre Harsányi, se seleccionan y transforman las variables de temperatura máxima, mínima y media, precipitación y radiación de onda corta del conjunto de datos original para construir la variable objetivo D_anom, derivada del cálculo de la evapotranspiración potencial y del balance hídrico climático, incorporando además el índice SPEI-3 como predictor adicional.

• Evapotranspiración Potencial (PET): Se calcula mediante el método de Hargreaves, el cual depende de la temperatura máxima, temperatura mínima y la radiación extraterrestre, proporcionando una estimación robusta de la demanda atmosférica de agua a escala mensual.
• Balance Hídrico Climático (D): Variable obtenida como la diferencia entre la precipitación mensual y la evapotranspiración potencial que representa el superávit o déficit hídrico mensual y constituye la base física para la evaluación de condiciones de sequía.
• SPEI3: Se calculó a partir del balance hídrico acumulado a una escala temporal de tres meses , y estandarizado mediante una distribución log-logística, con el fin de caracterizar condiciones de sequía y humedad a corto plazo.
• Anomalía del Balance Hídrico (D_anom): Se definió como la diferencia entre el balance hídrico mensual y su climatología mensual, permitiendo eliminar la estacionalidad climática y resaltar anomalías de déficit o exceso hídrico, que constituyen la variable objetivo del modelo de aprendizaje automático.

Mientras el SPEI-3 se utiliza como un indicador estandarizado de sequía a corto plazo, la anomalía del balance hídrico permite modelar directamente las desviaciones climáticas mensuales, facilitando su predicción mediante algoritmos de aprendizaje automático.

3. METODOLOGÍA Y RESULTADOS

3.1. Carga y revisión de estructura

Como primer punto, cargamos las librerías que vamos a usar en el estudio.

#Activamos librerías principales
library(tidymodels) #paquete para trabajar aprendizaje automático
library(ranger) # para trabajar con Random Forest
library(RKaggle) #paquete para leer el dataset desde Kaggle
library(tidyverse) #paquete para manipulación de información
library(lubridate) #paquete para trabaajr multitemporal
library(SPEI) #paquete para calcular método de hargreaves

Cargamos el dataset de clima que está alojado directamente en Kaggle. Su fuente oficial es OPEN-METEO. Son series temporales de 1995-2024. Este dataset tiene un total de 2’065.949 observaciones y 20 variables.

# Listar archivos del dataset
kWeatherdata <- get_dataset("wafaaelhusseini/capital-weather-data-1995-2024") #ruta del dataset en el sitio de kaggle.

Creamos un nuevo dataset seleccionando las variables que vamos a utilizar y filtramos por la capital europea Madrid - España. El dataset se redujo a 10.958 observaciones y 7 variables.

#selección de variables y filtrado
kWeatherdata_Madrid <-  kWeatherdata |> 
  select(date, capital, lat, lon, temp_min_c, temp_max_c, temp_mean_c_approx, precip_mm, shortwave_radiation_MJ_m2) |> 
  filter(capital == "Madrid")

Revisamos la estructura del dataset creado para comprender mejor como realizar el procesaminto de los datos. Vemos que la variable date y capital son de tipo character, mientras que, las variables meteorológicas son de tipo numéricas.

#estructura del dataset
str(kWeatherdata_Madrid)

## tibble [10,958 × 9] (S3: tbl_df/tbl/data.frame)
##  $ date                     : chr [1:10958] "1995-01-01" "1995-01-02" "1995-01-03" "1995-01-04" ...
##  $ capital                  : chr [1:10958] "Madrid" "Madrid" "Madrid" "Madrid" ...
##  $ lat                      : num [1:10958] 40.4 40.4 40.4 40.4 40.4 ...
##  $ lon                      : num [1:10958] -3.7 -3.7 -3.7 -3.7 -3.7 ...
##  $ temp_min_c               : num [1:10958] 0.3 -1.4 -1.5 -2 2.6 2.5 2.1 -1.5 0.5 -0.6 ...
##  $ temp_max_c               : num [1:10958] 12 6.7 5.4 5.9 10.8 9.3 8 9.5 10.3 10.5 ...
##  $ temp_mean_c_approx       : num [1:10958] 6.15 2.65 1.95 1.95 6.7 5.9 5.05 4 5.4 4.95 ...
##  $ precip_mm                : num [1:10958] 11.6 0 0 0.1 0 0 0 0 0 0 ...
##  $ shortwave_radiation_MJ_m2: num [1:10958] 8.33 8.86 6.64 4.38 5.71 8.21 6.87 9.1 8.34 8.78 ...

3.2. Análisis exploratorio de datos

Antes de realizar transformaciones al dataset, revisamos ciertos estadísticos para entender su distribución. Se trabajó con la función summary () para un análisis exploratorio general de los datos. Vemos que la temperatura mínimma llego a -11.20 °C y la temperatura máxima fue 41.3 °C. Lo que se traduce a un clima bastante frío en la temporada de invierno y muy caluroso en la temporada de verano. La temperatura media fue de 14.65 °C de clima templado.

Respecto a la precipitación, tenemos valores de 0 mm de lluvia en un día, con una media de 1.171 mm que se considera seco con valores máximos de 80.60 mm por día que podrían ser valores anormales de precipitación.

Por último, la mínima energía solar recibida en magnitud de Megajulios por m2 fue de 0.75 MJ/m2, valores muy bajos típicos de cielos nublados. También se recibieron valores medios de 17.76 MJ/m2 y máximos de 32,35 MJ/m2 que se traduce a días soleados con insolación media.

El dataset no presenta valores nulos para imputar.

#resumen estadístico
summary(kWeatherdata_Madrid)

##      date             capital               lat             lon        
##  Length:10958       Length:10958       Min.   :40.42   Min.   :-3.703  
##  Class :character   Class :character   1st Qu.:40.42   1st Qu.:-3.703  
##  Mode  :character   Mode  :character   Median :40.42   Median :-3.703  
##                                        Mean   :40.42   Mean   :-3.703  
##                                        3rd Qu.:40.42   3rd Qu.:-3.703  
##                                        Max.   :40.42   Max.   :-3.703  
##    temp_min_c        temp_max_c   temp_mean_c_approx   precip_mm     
##  Min.   :-11.200   Min.   : 0.4   Min.   :-2.90      Min.   : 0.000  
##  1st Qu.:  3.400   1st Qu.:12.5   1st Qu.: 8.05      1st Qu.: 0.000  
##  Median :  8.500   Median :19.2   Median :13.70      Median : 0.000  
##  Mean   :  9.092   Mean   :20.2   Mean   :14.65      Mean   : 1.171  
##  3rd Qu.: 15.000   3rd Qu.:27.8   3rd Qu.:21.25      3rd Qu.: 0.300  
##  Max.   : 25.200   Max.   :41.3   Max.   :32.70      Max.   :80.600  
##  shortwave_radiation_MJ_m2
##  Min.   : 0.75            
##  1st Qu.: 9.54            
##  Median :17.13            
##  Mean   :17.16            
##  3rd Qu.:24.78            
##  Max.   :32.35

Para poder realizar análisis multitemporal necesitamos que la variable date sea de tipo date () y con ello poder manipular las fechas correctamente y hacer operaciones temporales con ellas. El ingreso de la fecha puede variar en función del país y los formatos establecidos en cada caso. Por lo que, procedemos a revisar la estructura que contiene la variable date.

#primetras 10 filas
head(kWeatherdata_Madrid, 10)

## # A tibble: 10 × 9
##    date   capital   lat   lon temp_min_c temp_max_c temp_mean_c_approx precip_mm
##    <chr>  <chr>   <dbl> <dbl>      <dbl>      <dbl>              <dbl>     <dbl>
##  1 1995-… Madrid   40.4 -3.70        0.3       12                 6.15      11.6
##  2 1995-… Madrid   40.4 -3.70       -1.4        6.7               2.65       0  
##  3 1995-… Madrid   40.4 -3.70       -1.5        5.4               1.95       0  
##  4 1995-… Madrid   40.4 -3.70       -2          5.9               1.95       0.1
##  5 1995-… Madrid   40.4 -3.70        2.6       10.8               6.7        0  
##  6 1995-… Madrid   40.4 -3.70        2.5        9.3               5.9        0  
##  7 1995-… Madrid   40.4 -3.70        2.1        8                 5.05       0  
##  8 1995-… Madrid   40.4 -3.70       -1.5        9.5               4          0  
##  9 1995-… Madrid   40.4 -3.70        0.5       10.3               5.4        0  
## 10 1995-… Madrid   40.4 -3.70       -0.6       10.5               4.95       0  
## # ℹ 1 more variable: shortwave_radiation_MJ_m2 <dbl>

La variable date tiene la estructura ymd (year, month, day). Por lo que, la transformamos al formato fecha (date) para su manipulación. Adicionalmente, creamos tres variables más para segregar la fecha por año, mes y día.

#Transformamos el tipo de dato a fecha (date)
kWeatherdata_Madrid <- kWeatherdata_Madrid |> 
  mutate(date = ymd(date),
         año  = year(date),
         mes  = month(date, label = FALSE, abbr = FALSE, locale = "es_ES.UTF-8"))

TEMPORALIDAD ANUAL

Agregamos la data de forma anual para poder graficar las diferentes variables climáticas.

#agregación anual de los datos
data_anual <- kWeatherdata_Madrid |>
  group_by(año) |>
  summarise(
    Tmax = max(temp_max_c, na.rm = TRUE),
    Tmin = min(temp_min_c, na.rm = TRUE),
    Tavg = mean(temp_mean_c_approx, na.rm = TRUE),
    Rainfall = sum(precip_mm, na.rm = TRUE),
    Energy = mean(shortwave_radiation_MJ_m2, na.rm = TRUE)) |> 
  ungroup() |> 
  arrange(año)

La precipitación total en Madrid es irregular año tras año, típica de un clima mediterráneo continental con inviernos húmedos y veranos secos. No hay una tendencia lineal clara de aumento o disminución sostenida; en cambio, hay fluctuaciones marcadas, con picos y valles.

Las precipitaciones oscilan entre aproximadamente 200-250 mm (años secos) y 500-600 mm (años muy húmedos). El promedio aproximado parece estar alrededor de 350-400 mm anuales, lo que es consistente con los datos históricos de Madrid (promedio oficial ~450 mm, pero variable).

# Gráfico de precipitación anual
ggplot(data_anual, aes(x = año, y = Rainfall)) +
  geom_line(linewidth = 1, color = "#277da1") +
  geom_point(size = 2) +
  labs(
    title = "Precipitación total anual",
    subtitle = "1995–2024 (MADRID - ESPAÑA)",
    x = "Año",
    y = "Precipitación total (mm)"
  ) +
  scale_x_continuous(breaks = data_anual$año) +
  theme_minimal(base_family = "Tahoma") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16),
    plot.subtitle = element_text(hjust = 0.5, size = 13),
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10),
    axis.text.x = element_text(size = 8, angle = 90),
    axis.text.y = element_text(size = 8)
  )

El gráfico de temperatura máxima anual en Madrid (1995-2024) muestra una clara y marcada tendencia al alza a lo largo de las tres últimas décadas. Desde valores que rondaban los 36-37 °C a mediados de los años 90, las máximas anuales han pasado a situarse habitualmente por encima de los 40 °C a partir de 2012, con un salto especialmente evidente desde 2015 en adelante.

Los años más extremos han sido 2017, 2022 y especialmente 2023-2024, que superan los 43-44 °C de temperatura máxima registrada en algún momento del año. Aunque sigue habiendo cierta variabilidad interanual, la línea de tendencia es inequívocamente ascendente y los valores extremos se han vuelto mucho más frecuentes e intensos, confirmando que Madrid está experimentando un calentamiento significativo y un aumento notable de las olas de calor extremas en los últimos 10-15 años.

# Gráfico de temperatura máxima anual
ggplot(data_anual, aes(x = año, y = Tmax)) +
  geom_line(linewidth = 1, color = "red") +
  geom_point(size = 2) +
  labs(
    title = "Temperatura máxima",
    subtitle = "1995–2024 (MADRID - ESPAÑA)",
    x = "Año",
    y = "Temperatura Máxima (°C)"
  ) +
  scale_x_continuous(breaks = data_anual$año) +
  theme_minimal(base_family = "Tahoma") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16),
    plot.subtitle = element_text(hjust = 0.5, size = 13),
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10),
    axis.text.x = element_text(size = 8, angle = 90),
    axis.text.y = element_text(size = 8)
  )

El gráfico de temperatura mínima anual en Madrid (1995-2024) refleja un aumento muy pronunciado y acelerado del calor nocturno y de invierno en las últimas décadas. Hasta aproximadamente 2010 las mínimas más bajas del año oscilaban mayoritariamente entre -2 °C y +3 °C, con algún año excepcionalmente frío por debajo de -4 °C. Sin embargo, a partir de 2015 se observa una tendencia ascendente clara y, sobre todo, un salto dramático desde 2020.

Este fuerte incremento de las mínimas indica que las noches y madrugadas frías son cada vez más raras, los inviernos se han suavizado notablemente y las olas de calor nocturnas son más intensas y prolongadas, siendo uno de los indicadores más claros y preocupantes del calentamiento en la ciudad de Madrid.

#Gráfico de temperatura mínima anual
ggplot(data_anual, aes(x = año, y = Tmin)) +
  geom_line(linewidth = 1, color = "red") +
  geom_point(size = 2) +
  labs(
    title = "Temperatura mínima",
    subtitle = "1995–2024 (MADRID - ESPAÑA)",
    x = "Año",
    y = "Temperatura Mínima (°C)"
  ) +
  scale_x_continuous(breaks = data_anual$año) +
  theme_minimal(base_family = "Tahoma") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16),
    plot.subtitle = element_text(hjust = 0.5, size = 13),
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10),
    axis.text.x = element_text(size = 8, angle = 90),
    axis.text.y = element_text(size = 8)
  )

El gráfico de radiación solar anual en Madrid (1995-2024) muestra una fuerte tendencia al alza a lo largo de todo el periodo, especialmente marcada desde principios de los años 2000. A mediados de los 90 los valores rondaban los 16,7-17,0 MJ/m² de media anual, mientras que desde 2015 se han estabilizado claramente por encima de los 17,5-17,9 MJ/m², con picos que superan los 18,0 MJ/m² en años como 2003, 2017 o 2022. Esto supone un aumento aproximado de +1,0 a +1,3 MJ/m² (equivalente a un 6-8 % más de radiación anual) en solo tres décadas.

#Gráfico de radiación solar
ggplot(data_anual, aes(x = año, y = Energy)) +
  geom_line(linewidth = 1, color = "yellow") +
  geom_point(size = 2) +
  labs(
    title = "Radiación solar",
    subtitle = "1995–2024 (MADRID - ESPAÑA)",
    x = "Año",
    y = "Radiación Solar (Mj/m2)"
  ) +
  scale_x_continuous(breaks = data_anual$año) +
  theme_minimal(base_family = "Tahoma") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16),
    plot.subtitle = element_text(hjust = 0.5, size = 13),
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10),
    axis.text.x = element_text(size = 8, angle = 90),
    axis.text.y = element_text(size = 8)
  )

TEMPORALIDAD MENSUAL

Agregamos la data de forma mensual para poder graficar las diferentes variables climáticas.

#agregación mensual de los datos
data_mensual <- kWeatherdata_Madrid |>
  group_by(año, mes) |>
  summarise(
    Tmax = max(temp_max_c, na.rm = TRUE),
    Tmin = min(temp_min_c, na.rm = TRUE),
    Tavg = mean(temp_mean_c_approx, na.rm = TRUE),
    Rainfall = sum(precip_mm, na.rm = TRUE),
    Energy = mean(shortwave_radiation_MJ_m2, na.rm = TRUE)) |> 
  ungroup() |> 
  arrange(año, mes)

Precipitación:

En casi todos los años, se observa un ciclo anual estable, donde:

• Los meses con mayor precipitación tienden a concentrarse consistentemente en el mismo periodo del año.
• Los meses secos mantienen un patrón similar a lo largo del tiempo.
• Algunos años muestran picos mucho más altos (eventos lluviosos intensos).
• Otros años presentan niveles más bajos en los meses tradicionalmente húmedos.

Esto indica que la región presenta un régimen climático estacional definido, donde los ciclos de lluvia se repiten anualmente con relativa regularidad. Los eventos extremos que se visualizan en la serie temporal puede deberse a fenómenos climáticos de gran escala como El Niño-Oscilación del Sur (ENSO) u otros.

# Gráfico anual de precipitación 
ggplot(data_mensual, aes(x = mes, y = Rainfall)) +
  geom_line(color = "#277da1", linewidth = 0.8) +
  geom_point(size = 1) +
  facet_wrap(~año, scales = "free_x", nrow = 5, ncol = 6) +
  labs(
    title = "Evolución MENSUAL de la precipitación",
    subtitle = "1995–2024 (MADRID - ESPAÑA)",
    x = "Mes",
    y = "Precipitación (mm)"
  ) +
  scale_x_continuous(
    breaks = 1:12,
    labels = month(1:12, label = TRUE, abbr = TRUE)
  ) +
  theme_minimal(base_family = "Tahoma") +
  theme(
    plot.title = element_text(size = 20, hjust = 0.5),
    plot.subtitle = element_text(size = 14, hjust = 0.5),
    axis.title = element_text(size = 12),
    axis.text.x = element_text(size = 6, angle = 90, vjust = 0.5),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 10, face = "bold")   # título de cada faceta (cada año)
  )

Temperatura Máxima:

En todos los años se observa un patrón estacional cásico del clima mediterráneo continental. - Meses fríos: enero, febrero y diciembre con temperaturas máximas de entre 8°C y 15°C. - Meses cálidos: junio, julio y agosto con temperaturas máximas de entre 28°C y 40°C. - Este patrón se repite con notable regularidad, lo que indica una fuerte estacionalidad términa en Madrid.

# Gráfico anual de temperatura máxima
ggplot(data_mensual, aes(x = mes, y = Tmax)) +
  geom_line(color = "red", linewidth = 0.8) +
  geom_point(size = 1) +
  facet_wrap(~año, scales = "free_x", nrow = 5, ncol = 6) +
  labs(
    title = "Evolución MENSUAL de la temperatura máxima",
    subtitle = "1995–2024 (MADRID - ESPAÑA)",
    x = "Mes",
    y = "Temperatura Máxima (°C)"
  ) +
  scale_x_continuous(
    breaks = 1:12,
    labels = month(1:12, label = TRUE, abbr = TRUE)
  ) +
  theme_minimal(base_family = "Tahoma") +
  theme(
    plot.title = element_text(size = 20, hjust = 0.5),
    plot.subtitle = element_text(size = 14, hjust = 0.5),
    axis.title = element_text(size = 12),
    axis.text.x = element_text(size = 6, angle = 90, vjust = 0.5),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 10, face = "bold")   # título de cada faceta (cada año)
  )

Temperatura Mínima:

El comportamiento anual mínimo es consistente, igual que temperatura máxima, clásico del clima mediterráneo continental. - Meses más fríos: diciembre, enero y febrero con temperaturas mínimas de entre -11°C y 5°C. - Meses cálidos: junio, julio y agosto con temperaturas mínimas de entre 15°C y 23°C. - Este comportamiento refleja el clima mediterráneo de Madrid, con inviernos fríos y veranos muy cálidos.

# Gráfico anual de temperatura mínima
ggplot(data_mensual, aes(x = mes, y = Tmin)) +
  geom_line(color = "red", linewidth = 0.8) +
  geom_point(size = 1) +
  facet_wrap(~año, scales = "free_x", nrow = 5, ncol = 6) +
  labs(
    title = "Evolución MENSUAL de la temperatura mínima",
    subtitle = "1995–2024 (MADRID - ESPAÑA)",
    x = "Mes",
    y = "Temperatura Mínima (°C)"
  ) +
  scale_x_continuous(
    breaks = 1:12,
    labels = month(1:12, label = TRUE, abbr = TRUE)
  ) +
  theme_minimal(base_family = "Tahoma") +
  theme(
    plot.title = element_text(size = 20, hjust = 0.5),
    plot.subtitle = element_text(size = 14, hjust = 0.5),
    axis.title = element_text(size = 12),
    axis.text.x = element_text(size = 6, angle = 90, vjust = 0.5),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 10, face = "bold")   # título de cada faceta (cada año)
  )

Radiación solar:

El gráfico muestra la distribución mensual de la radiación solar global (MJ/m²) durante 30 años. Cada panel corresponde a un año y permite comparar la dinámica estacional, los picos anuales y las posibles variaciones relacionadas con el clima mediterráneo continentalizado de Madrid.

• Mínimos de radiación en los meses de diciembre y enero debido a la menor duración del día.
• Máximos de radiación en los meses de junio y julio que coinciden con el solsticio de verano y mayor insolación diaria.
• Este comportamiento reafirma el patrón estacional del clima de Madrid, donde la radiación solar es determinante para la temperatura y evapotranspiración.

# Gráfico anual de radiación solar
ggplot(data_mensual, aes(x = mes, y = Energy)) +
  geom_line(color = "yellow", linewidth = 0.8) +
  geom_point(size = 1) +
  facet_wrap(~año, scales = "free_x", nrow = 5, ncol = 6) +
  labs(
    title = "Evolución MENSUAL de la radiación solar",
    subtitle = "1995–2024 (MADRID - ESPAÑA)",
    x = "Mes",
    y = "Radiación solar (MJ/m²)"
  ) +
  scale_x_continuous(
    breaks = 1:12,
    labels = month(1:12, label = TRUE, abbr = TRUE)
  ) +
  theme_minimal(base_family = "Tahoma") +
  theme(
    plot.title = element_text(size = 20, hjust = 0.5),
    plot.subtitle = element_text(size = 14, hjust = 0.5),
    axis.title = element_text(size = 12),
    axis.text.x = element_text(size = 6, angle = 90, vjust = 0.5),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 10, face = "bold")   # título de cada faceta (cada año)
  )

3.3. Cáculo de Indicadores Climáticos (PET, D, D_anom, SPEI3)

Siguiendo la ecuación de Hargreaves citada en el estudio, calculamos la Evapotranspiración Potencial (PET) utilizando el algoritmo integrado en R de hargreaves agregado de forma mensual, que utiliza la fórmula siguiente:

\[PET = 0.0023 \cdot R_a \cdot (T_{max} - T_{min})^{0.5} \cdot (T_{avg} + 17.8)\] Posteriormente calculamos el Balance Hídrico Climático (D) en Madrid 1995-2024 de la siguiente forma:

\[D = Rainfall - PET\] Procedemos con el cálculo de la Evapotranspiración potencial y balance hídrico climático con las series temporales agregadas de forma mensual.

# Cálculo de PET (Método Hargreaves), D}
data_indicadores <- kWeatherdata_Madrid |>
  rename(year = año, month = mes) |> 
  group_by(year, month) |> 
  summarise(
    Rainfall = sum(precip_mm, na.rm = TRUE),
    Tmin = mean(temp_min_c, na.rm = TRUE),
    Tmax = mean(temp_max_c, na.rm = TRUE),
    Energy = mean(shortwave_radiation_MJ_m2, nna.rm = TRUE),
    Tmean = mean(temp_mean_c_approx, na.rm = TRUE),
    .groups = "drop") |> 
  mutate(
    PET = as.numeric(
      hargreaves(
        Tmin = Tmin,
        Tmax = Tmax,
        lat = 40.42)), #latitud de Madrid necesaria para la radiación.
    D = Rainfall - PET) # # balance hídrico climático mensual (mm)

## [1] "Calculating reference evapotranspiration using the Hargreaves method. Using latitude (`lat`) to estimate extraterrestrial radiation. Checking for missing values (`NA`): all the data must be complete. Input type is vector. Assuming the data are monthly time series starting in January, all regular (non-leap) years."

Adicionalmente, se calculó la anomalía del balance hídrico mediante la expresión:

\[Danom,t​=Dt​−\overline{D}_m\]

donde \(D_t\) representa el balance hídrico mensual en el tiempo \(t\) y \(\overline{D}_m\) corresponde al valor medio del balance hídrico para el mes \(m\), calculado sobre el período de referencia.

#media del balance hídrico mensual (mm)
data_indicadores <- data_indicadores |>
  group_by(month) |> #agregación mensual
  mutate(D_clim = mean(D, na.rm = TRUE)) |> # media del valor
  ungroup()

#anomalía D
data_indicadores <- data_indicadores |>
  mutate(D_anom = D - D_clim) #anomalia 

El índice SPEI-3 fue calculado a partir del balance hídrico climático mensual (P − PET), utilizando el algoritmo spei() aplicado a una serie temporal mensual, con una ventana de acumulación de tres meses y una distribución log-logística. Se aplicaron diferentes pasos:

1.- Se ordenó la serie temporal

#ordenamos la serie mensual
data_indicadores <- data_indicadores |>
  arrange(year, month)

2.- Se creó la serie temporal

#creamos la serie temporal
D_ts <- ts(
  data_indicadores$D,
  start = c(min(data_indicadores$year), min(data_indicadores$month)),
  frequency = 12
)

3.- Calculamos SPEI3 con la escala temporal 3 meses.

#calculamos SPEI3
spei3 <- spei(
  D_ts,
  scale = 3, #escala temporal
  distribution = "log-Logistic", # distribución logística para estandarizar valores.
  na.rm = TRUE
)

## [1] "Calculating the Standardized Precipitation Evapotranspiration Index (SPEI) at a time scale of 3. Using kernel type 'rectangular', with 0 shift. Fitting the data to a log-Logistic distribution. Using the ub-pwm parameter fitting method. Missing values (`NA`) will not be considered in the calculation. Using the whole time series as reference period. Input type is tsvector. Time series spanning ene 1995 to dic 2024, with frequency = 12."

4.- Se agregó el calculo a la serie original

#agregamos resultado al dataset original
data_indicadores$SPEI3 <- as.numeric(spei3$fitted)

5.- Eliminamos valores NA propios de la agregación trimestral de los primeros 2 meses de 1995.

#eliminamos valores NA para continuar con RF
data_indicadores <- data_indicadores |>
  filter(!is.na(SPEI3))

El índice SPEI3 presentó una distribución aproximadamente simétrica, centrada en cero (media = 0.01; mediana = −0.01), con valores extremos dentro de rangos esperados (−2.35 a 2.34), lo que confirma la correcta implementación de la estandarización.

summary(data_indicadores$SPEI3)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -2.350621 -0.763276 -0.008473  0.008210  0.736939  2.341977

El balance hídrico mensual D presentó valores mayoritariamente negativos, con una mediana de −49 mm, reflejando un régimen climático deficitario característico del área de estudio. Los valores positivos se concentraron en meses húmedos, confirmando la coherencia física del cálculo del balance hídrico.

summary(data_indicadores$D)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -210.850 -133.654  -49.329  -60.978   -1.628  149.352

El boxplot del SPEI-3 muestra una distribución centrada en torno a cero, lo que indica predominio de condiciones climáticas normales durante el período analizado. El rango intercuartílico es equilibrado, reflejando una variabilidad similar entre condiciones secas y húmedas moderadas, mientras que los valores extremos alcanzan umbrales cercanos a ±2, evidenciando la ocurrencia ocasional de eventos de sequía severa y humedad intensa. En conjunto, la figura confirma una correcta estandarización del índice y una representación coherente de la variabilidad climática a corto plazo.

boxplot(
  data_indicadores$SPEI3,
  main = "Boxplot SPEI-3"
)

La anomalía del balance hídrico D_anom presentó una media cercana a cero y una mediana ligeramente negativa, lo que indica una correcta eliminación de la estacionalidad y una leve predominancia de condiciones deficitarias. Los valores extremos positivos reflejan episodios húmedos intensos, consistentes con la variabilidad climática mensual del área de estudio.

summary(data_indicadores$D_anom)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -65.78217 -23.47730  -4.96970   0.04654  15.65318 140.02022

La distribución de la anomalía del balance hídrico (\(D_{\text{anom}}\)), analizada mediante un diagrama de caja, se centra en una mediana cercana a cero, lo que indica que el valor típico de la anomalía es la normalidad. El 50% central de los datos se distribuye simétricamente en un rango moderado.

El análisis de la anomalía del balance hídrico evidencia un predominio de eventos deficitarios, los cuales representan el 58.4 % del total de registros, frente al 41.6 % correspondiente a eventos de superávit. Esta distribución confirma que el régimen climático del área de estudio es mayoritariamente deficitario, con episodios húmedos menos frecuentes y de carácter puntual.

data_indicadores |>
  mutate(tipo_evento = case_when(
    D_anom > 0 ~ "Superávit",
    D_anom < 0 ~ "Déficit",
    TRUE ~ "Neutro"
  )) |>
  count(tipo_evento) |>
  mutate(
    porcentaje = round(100 * n / sum(n), 1)
  )

## # A tibble: 2 × 3
##   tipo_evento     n porcentaje
##   <chr>       <int>      <dbl>
## 1 Déficit       209       58.4
## 2 Superávit     149       41.6

Las mayores anomalías positivas del balance hídrico se concentraron en eventos puntuales de alta magnitud, destacando septiembre de 2023 como el mes con el mayor superávit hídrico registrado. Otros episodios relevantes ocurrieron en marzo de 2022 y en los meses de noviembre de 1997 y diciembre de 1996, evidenciando la ocurrencia ocasional de eventos húmedos extremos dentro de un régimen climático predominantemente deficitario.

data_indicadores |>
  arrange(desc(D_anom)) |>
  head(n = 5) |> 
  select(year, month, D_anom)

## # A tibble: 5 × 3
##    year month D_anom
##   <dbl> <dbl>  <dbl>
## 1  2023     9   140.
## 2  1996    12   127.
## 3  2022     3   114.
## 4  1997    11   112.
## 5  2021     9   104.

Las anomalías negativas más extremas del balance hídrico se registraron en abril de 2023, cinco meses antes del superávit hídrico observado en septiembre de 2023. Otros episodios severos ocurrieron en octubre de 1995, marzo de 1997, mayo de 2015 y abril de 1994, reflejando periodos de sequía intensa y persistente dentro de la serie temporal analizada, con magnitudes claramente inferiores a las anomalías positivas extremas pero con mayor frecuencia de ocurrencia.

data_indicadores |>
  arrange((D_anom)) |>
  head(n = 5) |> 
  select(year, month, D_anom)

## # A tibble: 5 × 3
##    year month D_anom
##   <dbl> <dbl>  <dbl>
## 1  2023     4  -65.8
## 2  1995    10  -63.5
## 3  1997     3  -63.1
## 4  2015     5  -53.6
## 5  1995     4  -51.3

No obstante, se observa una clara asimetría en los valores extremos, con la presencia exclusiva de outliers de alta magnitud en el extremo positivo (superávit hídrico), superando el percentil 99 y alcanzando valores máximos de aproximadamente 140. Esta distribución sugiere que los eventos de gran anomalía positiva son raros pero significativos, mientras que no se registraron déficits extremos comparables.

boxplot(
  data_indicadores$D_anom,
  main = "Boxplot D_anom"
)

En conjunto, los indicadores climáticos calculados evidencian un régimen hídrico caracterizado por un déficit estacional dominante, reflejado en los valores mayoritariamente negativos del balance hídrico mensual (D), coherente con las condiciones climáticas del área de estudio. La transformación a anomalías del balance hídrico (D_anom) permitió eliminar eficazmente la estacionalidad, obteniendo una variable centrada en cero que resalta episodios anómalos de déficit y superávit hídrico, siendo estos últimos menos frecuentes pero de mayor magnitud.

Por su parte, el índice SPEI-3 mostró una distribución estandarizada y simétrica, adecuada para la identificación de eventos de sequía a corto plazo. En conjunto, estos resultados confirman la coherencia física y estadística de los indicadores empleados y justifican el uso de D_anom como variable objetivo para la modelación predictiva, mientras que D y SPEI-3 aportan información complementaria sobre la dinámica hidrológica y la persistencia temporal de las condiciones de sequía.

3.4. Modelado con Random Forest y Tidymodels

3.4.1. División de datos (Time Series Split)

Siguiendo las recomendaciones metodológicas de Harsányi (2025), se aplicó un esquema de validación cruzada junto con una partición explícita de los datos en conjuntos de entrenamiento y prueba. En este estudio se adoptó una división del 70% de los datos para entrenamiento y el 30% para prueba, considerando que el conjunto de datos, agregado a escala mensual, presenta una extensión temporal limitada. Este enfoque permite evaluar la capacidad de generalización del modelo y reducir el riesgo de sobreajuste en series temporales relativamente cortas.

set.seed(123)  # Para reproducibilidad
data_split <- initial_split(data_indicadores, prop = 0.7, strata = D_anom) #dividimos el dataset para entrenar y testear
train_data <- training(data_split) #datos para training
test_data <- testing(data_split) # datos para testing

#cantidad de observaciones
cat("Observaciones Entrenamiento:", nrow(train_data), "\n")

## Observaciones Entrenamiento: 248

cat("Observaciones Prueba:", nrow(test_data))

## Observaciones Prueba: 110

La optimización del modelo se realizó mediante validación cruzada k-fold con k = 10, aplicada exclusivamente al conjunto de entrenamiento, con el fin de seleccionar los hiperparámetros óptimos y evaluar la estabilidad del modelo daua una limitada extensión temporal.

#dividiso en 5-folds para entrenar
p_folds <- vfold_cv(train_data, v=10)

3.4.2. Creación de la Receta

Para la construcción del flujo de modelación se definió una receta que especifica la relación entre la variable objetivo y las variables predictoras (se utilizaron todas las vairables del dataset como predictoras). En esta etapa, se incorporó la normalización de las variables predictoras con el objetivo de estandarizar su representación dentro del framework tidymodels, garantizando un preprocesamiento consistente y reproducible previo al entrenamiento del modelo.

# Crear la receta de preprocesamiento
rf_recipe <- recipe(D_anom ~ ., data = train_data) |> 
  step_normalize(all_predictors())  # Estandarizar las variables predictoras

3.4.3. Especificación del Modelo.

En esta etapa se especificó el modelo Random Forest con tuning, utilizando el motor ranger por ser rápido y compatible con tidymodels y la importancia de tipo permutation por ser más robusta y fiable para estudios científicos. Los árboles elegidos 300 fueron recomendados en el estudio de Harsányi (2025).

#creamos el modelo
rf_tune <- rand_forest(
  mtry = tune(),
  trees = 300, # menos árboles menor sobreajuste. 300 recomienda el estudio.
  min_n = tune()) |> 
  set_mode("regression") |> 
  set_engine("ranger", importance = "permutation") #mayor rigor científico

3.4.4. Workflow.

Inicializo el workflow para trabajar de manera ordenada, luego visualizo los pasos.

#creamos el workflow
tune_wf <- workflow() |> 
  add_recipe(rf_recipe) %>%
  add_model(rf_tune)

#revisamos el resultado
tune_wf

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: rand_forest()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 1 Recipe Step
## 
## • step_normalize()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Random Forest Model Specification (regression)
## 
## Main Arguments:
##   mtry = tune()
##   trees = 300
##   min_n = tune()
## 
## Engine-Specific Arguments:
##   importance = permutation
## 
## Computational engine: ranger

3.4.5. Entrenamiento del modelo con tuning

El desempeño del modelo se evaluó utilizando un conjunto de métricas complementarias, que incluyen el coeficiente de determinación (R²), el error absoluto medio (MAE) y la raíz del error cuadrático medio (RMSE).

Estas métricas permiten evaluar de forma conjunta la capacidad explicativa del modelo, la magnitud promedio del error y la penalización de errores de mayor magnitud en la predicción de la anomalía del balance hídrico.

#métricas para validar
metrics <- metric_set(rsq, mae, rmse)

Para garantizar la reproducibilidad de los resultados, se fijó una semilla aleatoria antes del proceso de entrenamiento. La optimización de los hiperparámetros del modelo Random Forest se realizó mediante búsqueda en rejilla (grid search), evaluando once (11) combinaciones distintas de hiperparámetros. Este proceso se ejecutó utilizando validación cruzada sobre el conjunto de entrenamiento y procesamiento en paralelo, con el fin de mejorar la eficiencia computacional y seleccionar la configuración que maximiza el desempeño del modelo según las métricas definidas.

#entrenamiento con los 10 folds
set.seed(123) # para reproducibilidad
doParallel::registerDoParallel() ##activa el procesamiento paralelo con tune_grid
ranger_tune <-tune_grid(tune_wf,
    resamples = p_folds, #tomo la variable del training folds.
    metrics=metrics,
    grid = 11 #evalua 11 combinaciones distintas de hiperparámetros. 
  )

El mejor modelo según el criterio del Error Absoluto Medio - MAE fue la configuración pre0_mod09_post0. Recojemos todas las métricas asociadas a esa configuración para evaluar el desempeño.

#mejor modelo
show_best(ranger_tune, metric = "mae")

## # A tibble: 5 × 8
##    mtry min_n .metric .estimator  mean     n std_err .config         
##   <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>           
## 1     9     2 mae     standard    7.45    10   0.432 pre0_mod09_post0
## 2    10    13 mae     standard    8.28    10   0.438 pre0_mod10_post0
## 3     6     9 mae     standard    8.36    10   0.454 pre0_mod06_post0
## 4     3     5 mae     standard    9.06    10   0.559 pre0_mod03_post0
## 5     7    24 mae     standard    9.50    10   0.548 pre0_mod07_post0

El modelo Random Forest optimizado mediante validación cruzada mostró un desempeño elevado en la predicción de la anomalía del balance hídrico mensual, alcanzando un coeficiente de determinación promedio de 0.91. El error absoluto medio fue de 7.45 mm y el RMSE de 10.34 mm, lo que indica una buena precisión en relación con el rango de variabilidad de la variable objetivo.

El modelo Random Forest utilizó 9 variables predictoras aleatoriamente mtry = 9 en cada división de los árboles manteniendo árboles profundos min_n = 2. En principio, el modelo puede aumentar el riesgo de sobreajuste con dataset pequeños como la agregación mensual que se tiene. Sin embargo, considerando el std_err que es bajo 0.013 para rsq, se podría asumir un desempeño estable. Adicionalmente, el bagging, el muestreo aleatorio y la agregación de árboles reduce significativamente el sobreajuste global.

#métricas completas.
ranger_tune |> 
  collect_metrics() |> 
  filter(.config == "pre0_mod09_post0")

## # A tibble: 3 × 8
##    mtry min_n .metric .estimator   mean     n std_err .config         
##   <int> <int> <chr>   <chr>       <dbl> <int>   <dbl> <chr>           
## 1     9     2 mae     standard    7.45     10  0.432  pre0_mod09_post0
## 2     9     2 rmse    standard   10.3      10  0.652  pre0_mod09_post0
## 3     9     2 rsq     standard    0.909    10  0.0129 pre0_mod09_post0

Guardamos la configuración dentro de un nuevo objeto para entrenar finalmente el modelo con los datos de test y la mejor configuración.

#guardamos mejor modelo
best_mae <- select_best(ranger_tune)
best_mae

## # A tibble: 1 × 3
##    mtry min_n .config         
##   <int> <int> <chr>           
## 1     9     2 pre0_mod09_post0

Finalizamos el workflow con el mejor modelo seleccionado.

final_rf <- finalize_workflow(
  tune_wf,
  best_mae
)

3.4.6. Predicción en TEST set

Una vez evaluado el desempeño del modelo en entrenamiento, procedemos a evaluar las métricas en testing para validar el desempeño con test dataset.

#metricas para evaluar en test
metrics <- metric_set(rsq, mae, rmse)

Entrenamos el modelo con el dataset de prueba usando last_fit y se pudo determinar que el desempeño del modelo se mantuvo estable al evaluarse sobre el conjunto de prueba independiente, alcanzando un R² de 0.89, un MAE de 7.82 mm y un RMSE de 10.80 mm. La similitud entre las métricas obtenidas en la validación cruzada y en el conjunto de prueba sugiere una adecuada capacidad de generalización y ausencia de sobreajuste.

set.seed(123) # para reproducibilidad
final_res <- last_fit(final_rf, data_split, metrics = metrics)
collect_metrics(final_res)

## # A tibble: 3 × 4
##   .metric .estimator .estimate .config        
##   <chr>   <chr>          <dbl> <chr>          
## 1 rsq     standard       0.891 pre0_mod0_post0
## 2 mae     standard       7.82  pre0_mod0_post0
## 3 rmse    standard      10.8   pre0_mod0_post0

Una vez entrenado el modelo final, se extrajeron las predicciones correspondientes al conjunto de prueba mediante la función collect_predictions(). Este procedimiento permitió obtener los valores predichos de la anomalía del balance hídrico (\(D_{\text{anom}}\)) y compararlos directamente con los valores observados, facilitando la evaluación del desempeño del modelo y el análisis de su capacidad de generalización.

#valores predichos y observados
predictions <- final_res |> 
  collect_predictions()

El modelo de predicción para la Anomalía del Balance Hídrico Climático (\(D_{\text{anom}}\)) demuestra ser robusto y preciso en la mayoría de los rangos, especialmente para déficits y anomalías cercanas a cero. La fuerte agrupación de puntos sobre la línea 1:1 es un indicador de un alto poder predictivo. Sin embargo, existe una limitación o incertidumbre mayor en la predicción de los eventos de superávit hídrico más extremos, donde el modelo tiende a tener una precisión ligeramente menor.

#gráfico de valores observados vs valores predichos
ggplot(predictions, aes(x = D_anom, y = .pred)) +
  geom_point(color = "blue", alpha = 0.6) +  # Puntos azules para las observaciones
  geom_abline(slope = 1, intercept = 0, color = "red") +  # Línea de igualdad
  labs(
    title = "Valores predichos VS Valores observados - D_anom",
    subtitle = "Anomalía del Balance Hídrico Climático",
    x = "Valores observados",
    y = "Valores predichos"
  ) +
  theme_minimal(base_family = "Tahoma") +
  theme(
    plot.title = element_text(size = 15, hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(size = 12, hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text.x = element_text(size = 10),
    axis.text.y = element_text(size = 10, face = "bold"))

Se calcularon los residuos del modelo que permitieron evaluar la adecuación del modelo y detectar posibles sesgos sistemáticos.

#residuos
predictions <- predictions |> 
  mutate(residuo = .pred - D_anom)

El diágnostico de residuos revela ciertos problemas de sesgos no lineales con ligeras sobreestimaciones con los déficits y subestimaciones con los superávits. Esto se traduce a que el modelo es menos fiable para predecir anomalías positivas extremas. Esto se pudo evidenciar con el boxplot de D_anom que contenía outliers positivos extremos.

#grafico valores observados vs residuos
ggplot(predictions, aes(x = D_anom, y = residuo)) +
  geom_point(alpha = 0.6, color = "blue") +
  geom_hline(yintercept = 0, color = "red") +
  geom_smooth(method = "loess", se = FALSE, color = "black") +
  labs(
    x = "Valores observado",
    y = "Residuo (Predicho − Observado)",
    title = "Residuos del modelo vs D_anom observado",
    subtitle = "Anomalía del Balance Hídrico Climático") +
  theme_minimal(base_family = "Tahoma") +
  theme(
    plot.title = element_text(size = 15, hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(size = 12, hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text.x = element_text(size = 10),
    axis.text.y = element_text(size = 10, face = "bold"))

3.4.7. Importancia de variables

La principal variable relevante para el modelo creado fue la precipitación seguida del promedio del balance hídrico climático, el balance hídrico climático, temperatura media. En menor medida el PET, SPEI3 y radiación.

#para obtener variables importantes
library(vip)

# Ajustamos modelo final
final_fit <- final_rf |> 
  fit(data = train_data)
# graficamos importania final
vip(pull_workflow_fit(final_fit)$fit, geom = "point") +
  labs(title = "Importancia de Variables (Random Forest)")

4. Conclusión

• Los resultados obtenidos permiten responder afirmativamente a la pregunta de investigación, evidenciando que es posible modelar anomalías del balance hídrico mensual mediante algoritmos de aprendizaje automático y variables climáticas.
• El modelo Random Forest cumplió el objetivo del estudio al presentar un alto desempeño predictivo y una adecuada capacidad de generalización, especialmente para condiciones normales y de déficit moderado. Aunque el modelo muestra limitaciones en la predicción de eventos extremos de superávit hídrico, su desempeño global confirma su utilidad como herramienta de apoyo para el análisis y monitoreo de la sequía.

5. Recomendación

Se podría eliminar los outliers (superhábit) identificados con el boxplot de D_anom para entrenar nuevamente el modelo y evaluar su desempeño y ver si mejoran los residuos y su distribución cercano a 0.