Curso-Taller: Análisis Hidrológico de Cuencas con R
Clase 3: Modelación hidrológica conceptual con GR4J (parte 3)
Harold Llauca
📦 Preparación del entorno de trabajo
Se configura el entorno en R, definiendo el directorio de trabajo y cargando los paquetes necesarios para asegurar un flujo de trabajo ordenado y reproducible.
Instalación y carga de paquetes
# Establecer mirror de CRAN
options(repos = c(CRAN = "https://cloud.r-project.org"))
# Lista de paquetes necesarios
packages <- c(
"openxlsx", # lectura y escritura de archivos Excel (.xlsx)
"zoo", # manejo de series temporales
"dplyr", # manipulación de datos (data frames)
"terra", # manejo de datos geoespaciales (raster y vector)
"airGR", # implementación de modelos hidrológicos tipo GR (GR4J, GR5J, etc.)
"rtop", # optimización SCE-UA y herramientas de regionalización hidrológica
"hydroGOF", # métricas de evaluación hidrológica (NSE, KGE, RMSE, PBIAS, etc.)
"hydroTSM", # análisis, agregación y visualización de series temporales hidrológicas
"airGRiwrm", # modelación hidrológica integrada con modelos GR en redes de cuencas
"tidyr", # transformación de datos entre formato ancho y largo
"ggplot2" # elaboración de gráficos
)
# Instalar paquetes si no están instalados
for (pkg in packages) {
if (!requireNamespace(pkg, quietly = TRUE)) {
install.packages(pkg)
}
}
# Cargar paquetes
invisible(lapply(packages, function(pkg) {
suppressPackageStartupMessages(library(pkg, character.only = TRUE))
}))## Warning: package 'ggplot2' was built under R version 4.5.2☁️ Cargar salidas del modelo calibrado-validado
El objeto Outputs generado por el modelo GR4J contiene información tanto de las variables de entrada como de los procesos internos y las salidas del modelo. A continuación, se describen sus principales componentes:
DatesR corresponde a las fechas de la simulación. Precip representa la precipitación total (mm/d) y PotEvap la evapotranspiración potencial (mm/d).
Prod indica el almacenamiento en el reservorio de producción (mm). Pn es la precipitación neta que ingresa al sistema luego de considerar la evapotranspiración potencial. Ps corresponde a la fracción de Pn que alimenta el reservorio de producción. AE es la evapotranspiración real y Perc representa la percolación desde el reservorio de producción hacia el sistema de enrutamiento.
PR es la precipitación efectiva que alimenta los hidrogramas unitarios del modelo.
Q9 y Q1 corresponden a los flujos generados por los hidrogramas unitarios UH1 y UH2, asociados a la respuesta lenta (≈90%) y rápida (≈10%) del sistema, respectivamente.
Rout representa el almacenamiento en el reservorio de enrutamiento, mientras que Exch es el intercambio de agua con el acuífero, que puede ser positivo o negativo. AExch1, AExch2 y AExch corresponden a los intercambios acumulados en las diferentes componentes del modelo.
QR es el caudal proveniente del reservorio de enrutamiento (flujo base), y QD es el caudal directo o escorrentía rápida. La suma de ambos da como resultado Qsim, que es el caudal total simulado por el modelo (mm/d).
RunOptions contiene la configuración utilizada para ejecutar el modelo, mientras que StateEnd guarda el estado final de los reservorios, útil como condición inicial para simulaciones posteriores.
Finalmente, Qobs corresponde al caudal observado, IdCal e IdVal son los índices de los periodos de calibración y validación, respectivamente, WarmUp indica el periodo de calentamiento, Param contiene los parámetros calibrados del modelo (X1, X2, X3, X4), y Outputs$area corresponde al área de la cuenca (km²).
#-----------------------------------------------------------
# Cargar salidas y extraer variables del modelo
#-----------------------------------------------------------
# Leer archivo RDS con las salidas del modelo
Outputs <- readRDS('./Modules/Clase_3/outputs/model_outputs_lumped.rds')
# Ver nombres de los elementos disponibles
names(Outputs)## [1] "DatesR" "PotEvap" "Precip" "Prod" "Pn"
## [6] "Ps" "AE" "Perc" "PR" "Q9"
## [11] "Q1" "Rout" "Exch" "AExch1" "AExch2"
## [16] "AExch" "QR" "QD" "Qsim" "RunOptions"
## [21] "StateEnd" "Qobs" "IdCal" "IdVal" "WarmUp"
## [26] "Param" "Area"# Extraer variables principales
area <- Outputs$Area # área de la cuenca [km2]
Qsim <- Outputs$Qsim # caudal simulado [mm/d]
Qobs <- Outputs$Qobs # caudal observado [mm/d]
Ind_cal <- Outputs$IdCal # índices del periodo de calibración
Ind_val <- Outputs$IdVal # índices del periodo de validación
wup <- Outputs$WarmUp # periodo de warm-up⚙️Desempeño del modelo
En esta sección se evalúa visualmente el desempeño del modelo GR4J utilizando las funciones gráficas de airGR. Estas figuras permiten analizar la coherencia de las variables simuladas, comparar los caudales observados y simulados, e identificar posibles sesgos en la estacionalidad, la distribución de caudales y la relación entre valores observados y simulados.
Primero se visualizan las principales variables hidrológicas simuladas por el modelo, incluyendo la precipitación, la evapotranspiración real y los flujos generados. Esta figura permite revisar si la respuesta del modelo es coherente con los eventos de lluvia.
#-----------------------------------------------------------
# Visualizar precipitación, evapotranspiración real y caudales
#-----------------------------------------------------------
# Mostrar precipitación, evapotranspiración real y componentes de caudal
plot(
Outputs,
which = c("Precip", "ActuEvap", "Flows")
)
Luego se compara el caudal observado con el caudal simulado. Además, se muestra el error de simulación, lo cual permite identificar periodos donde el modelo sobreestima o subestima los caudales observados.
#-----------------------------------------------------------
# Comparar caudal observado, caudal simulado y error
#-----------------------------------------------------------
# Mostrar hidrograma observado/simulado y error de simulación
plot(
Outputs,
Qobs = Outputs$Qobs,
which = c("Flows", "Error")
)## Warning in plot.OutputsModel(Outputs, Qobs = Outputs$Qobs, which = c("Flows", :
## zeroes detected in 'Qobs': some plots in the log space will not be created
## using all time-steps
Finalmente, se evalúa el desempeño del modelo desde una perspectiva estadística. El régimen medio permite comparar la estacionalidad observada y simulada; la curva de duración muestra si el modelo reproduce adecuadamente la distribución de caudales; y el gráfico observado-simulado permite analizar la dispersión y el sesgo general de la simulación.
#-----------------------------------------------------------
# Diagnóstico estadístico del desempeño del modelo
#-----------------------------------------------------------
# Mostrar régimen medio, curva de duración y gráfico observado-simulado
plot(
Outputs,
Qobs = Outputs$Qobs,
which = c("Regime", "CumFreq", "CorQQ"),
LayoutMat = matrix(c(1, 2, 3, 3), ncol = 2),
LayoutWidths = c(1.5, 1),
LayoutHeights = c(1, 1, 1)
)## Warning in plot.OutputsModel(Outputs, Qobs = Outputs$Qobs, which = c("Regime",
## : zeroes detected in 'Qobs': some plots in the log space will not be created
## using all time-steps
En esta sección se calculan estadísticas de desempeño del modelo GR4J para los periodos de calibración, validación y el periodo total. Para ello, se comparan los caudales simulados y observados, previamente convertidos a unidades de m³/s y excluyendo el periodo de calentamiento (warm-up). Finalmente, los resultados se exportan para su análisis y reporte.
#-----------------------------------------------------------
# Cálculo de métricas de desempeño del modelo
#-----------------------------------------------------------
stats <- list()
#---------------------------
# Calibración
#---------------------------
x <- Qsim[Ind_cal] * area / 86.4 # mm/d → m3/s
y <- Qobs[Ind_cal] * area / 86.4
x <- x[-wup] # remover warm-up
y <- y[-wup]
stats$calibracion <- gof(x, y)
#---------------------------
# Validación
#---------------------------
x <- Qsim[Ind_val] * area / 86.4
y <- Qobs[Ind_val] * area / 86.4
stats$validacion <- gof(x, y)## Warning: 'rNSE' can not be computed: some elements in 'obs' are zero !## Warning: 'wNSE' can not be computed: some elements in 'obs' are zero !## Warning: 'rd' can not be computed: some elements in 'obs' are zero !#---------------------------
# Periodo total
#---------------------------
x <- Qsim * area / 86.4
y <- Qobs * area / 86.4
x <- x[-wup]
y <- y[-wup]
stats$total <- gof(x, y)## Warning: 'rNSE' can not be computed: some elements in 'obs' are zero !## Warning: 'wNSE' can not be computed: some elements in 'obs' are zero !## Warning: 'rd' can not be computed: some elements in 'obs' are zero !#---------------------------
# Resumen
#---------------------------
stats <- do.call(cbind, stats)
colnames(stats) <- c("Calibración", "Validación", "Total")
print(stats)## Calibración Validación Total
## ME 2.99 4.75 3.55
## MAE 15.04 15.77 15.27
## MSE 1164.25 1110.69 1147.11
## RMSE 34.12 33.33 33.87
## ubRMSE 33.99 32.99 33.68
## NRMSE % 66.80 80.00 70.20
## PBIAS % 11.90 18.70 14.10
## RSR 0.67 0.80 0.70
## rSD 0.89 1.17 0.97
## NSE 0.55 0.36 0.51
## mNSE 0.49 0.42 0.47
## rNSE -21.52 -Inf -Inf
## wNSE 0.57 0.58 0.57
## wsNSE 0.36 0.30 0.34
## d 0.86 0.85 0.86
## dr 0.75 0.71 0.73
## md 0.73 0.71 0.73
## rd -6.12 -Inf -Inf
## cp -1.77 -3.13 -2.09
## r 0.76 0.75 0.75
## R2 0.55 0.36 0.51
## bR2 0.42 0.35 0.41
## VE 0.40 0.38 0.39
## KGE 0.71 0.64 0.71
## KGElf 0.13 NaN NaN
## KGEnp 0.70 0.64 0.68
## KGEkm 0.72 0.65 0.70#-----------------------------------------------------------
# Exportar resultados
#-----------------------------------------------------------
openxlsx::write.xlsx(
data.frame(Index = rownames(stats), stats),
"./Modules/Clase_3/outputs/model_gof_lumped.xlsx"
)Se presenta un gráfico resumen de las principales métricas de desempeño del modelo para los periodos de calibración, validación y total. Esto permite comparar rápidamente la eficiencia del modelo, el error medio y el sesgo en cada periodo evaluado.
#-----------------------------------------------------------
# Gráfico resumen de métricas de desempeño
#-----------------------------------------------------------
stats_plot <- stats %>%
as.data.frame() %>%
tibble::rownames_to_column("Metrica") %>%
dplyr::filter(Metrica %in% c("NSE", "KGE", "RMSE", "PBIAS %")) %>%
tidyr::pivot_longer(
cols = -Metrica,
names_to = "Periodo",
values_to = "Valor"
)
# Ordenar métricas
stats_plot$Metrica <- factor(
stats_plot$Metrica,
levels = c("KGE", "NSE", "RMSE", "PBIAS %")
)
# Ordenar periodos
stats_plot$Periodo <- factor(
stats_plot$Periodo,
levels = c("Calibración", "Validación", "Total")
)
# Crear gráfico
ggplot(stats_plot, aes(x = Periodo, y = Valor, fill = Periodo)) +
geom_col(width = 0.65, color = "black", linewidth = 0.2) +
facet_wrap(~ Metrica, scales = "free_y", ncol = 2) +
labs(
x = NULL,
y = NULL,
title = "Resumen de métricas de desempeño del modelo GR4J"
) +
theme_bw(base_size = 10) +
theme(
legend.position = "none",
plot.title = element_text(hjust = 0.5, face = "bold"),
strip.text = element_text(face = "bold"),
panel.grid.minor = element_blank()
)