---
title: "Curso-Taller: Análisis Hidrológico de Cuencas con R - 2026"
subtitle: "Clase 3: Modelación hidrológica conceptual con GR4J (parte 1)"
author: "Harold Llauca"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_depth: 3
    toc_float:
      collapsed: false
      smooth_scroll: true
    theme: cosmo
    highlight: tango
    
knitr:
  opts_chunk:
    echo: true
    warning: false
    message: false
---

# 📦 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.

```{r setup, include=FALSE}
# Configuración inicial del documento
knitr::opts_chunk$set(echo = TRUE)

# Cambiar esta ruta según la ubicación del proyecto en su equipo.
# Se recomienda usar rutas relativas para mayor portabilidad.

knitr::opts_knit$set(
  root.dir = "C:/Users/hllauca/Desktop/Analisis_Hidrologico_Cuencas_R_2026"
  )
```

## Instalación y carga de librerias

```{r}
# Establecer mirror de CRAN
options(repos = c(CRAN = "https://cloud.r-project.org"))

# Lista de paquetes necesarios
packages <- c(
  "Evapotranspiration", # cálculo de evapotranspiración potencial (Hargreaves-Samani)
  "openxlsx",           # lectura de archivos Excel (.xlsx)
  "zoo",                # manejo de series temporales
  "dplyr",              # manipulación de datos (data frames)
  "exactextractr",      # extracción precisa de valores raster sobre polígonos
  "terra",              # manejo de datos geoespaciales (raster y vector)
  "ggplot2",            # elaboración de gráficos
  "tidyr"               # transformación de datos a formato largo/ancho
)

# 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))
}))
```

# 🗺 Cargar datos geoespaciales de la cuenca

En esta sección se cargan los datos geoespaciales de la cuenca de estudio. Se utiliza un Modelo Digital de Elevación (DEM) y el polígono de la cuenca hidrográfica, ambos en sistema de coordenadas WGS84. Estos insumos serán fundamentales para la preparación de variables fisiográficas requeridas por el modelo hidrológico.

```{r}
# Cargar el Modelo Digital de Elevación (DEM)
dem  <- rast('./Data/SIG/dem.tif')  # Sistema de referencia: WGS84

# Cargar el shapefile de la cuenca
cat  <- vect('./Data/SIG/basin.shp')  # Sistema de referencia: WGS84

# Visualizar cuenca
plot(cat)
```

```{r}
# Calcular el área del polígono de la cuenca en km2
area <- terra::expanse(cat, unit = 'km')

# Mostrar el resultado
print(area)
```

# 📁️ Preparación de datos de entrada

Periodo de modelación hidrológica seleccionada:

```{r}
# Crear vector de tiempo
time <- seq(as.Date('1981-01-01'), as.Date('2020-12-31'), by='days')
```

## Series de precipitación diaria

En esta sección se carga la precipitación diaria del producto PISCO. Esta información corresponde a datos grillados en formato NetCDF, los cuales serán agregados espacialmente sobre el área de la cuenca para obtener una serie temporal representativa de precipitación media areal.

```{r}
# Cargar el archivo NetCDF de precipitación diaria
pisco_pr <- rast('./Data/PISCO/PISCOp_v3.nc')
```

Luego de cargar la precipitación grillada, se calcula la precipitación promedio sobre la cuenca. Para ello se realiza una extracción espacial ponderada por área, lo cual permite obtener un valor medio diario considerando la intersección entre las celdas del raster y el polígono de la cuenca.

```{r}
# Calcular la precipitación media diaria ponderada por área
dpr <- exact_extract(
  pisco_pr,
  as(cat, 'Spatial'),
  fun      = 'mean',
  weights  = 'area',
  progress = FALSE
)

# Transponer el resultado y convertirlo a data frame
dpr <- t(dpr) %>% 
  as.data.frame()

# Convertir la serie de precipitación a objeto zoo
prec <- zoo(dpr, time)

head(prec)
```

## Series de evapotranspiración potencial

En esta sección se cargan las temperaturas máximas y mínimas diarias del producto PISCOt v1.2 (1981-2020). Estas variables serán utilizadas para estimar la evapotranspiración potencial mediante el método de Hargreaves-Samani.

```{r}
# Cargar temperatura máxima diaria
pisco_tx <- rast('./Data/PISCO/PISCOtx_v1.2.nc')

# Cargar temperatura mínima diaria
pisco_tn <- rast('./Data/PISCO/PISCOtn_v1.2.nc')
```

A partir de los datos grillados de temperatura, se extraen los valores medios diarios sobre la cuenca. La extracción se realiza de forma ponderada por área para representar mejor la contribución de cada celda dentro del polígono de la cuenca.

```{r}
#-----------------------------------------------------------
# Extraer temperatura máxima y mínima media sobre la cuenca
#-----------------------------------------------------------

# Extraer temperatura máxima media diaria sobre la cuenca
dtx <- exact_extract(
  pisco_tx,
  as(cat, 'Spatial'),
  fun      = 'mean',
  weights  = 'area',
  progress = FALSE
) %>%
  t() %>%
  as.data.frame()

# Extraer temperatura mínima media diaria sobre la cuenca
dtn <- exact_extract(
  pisco_tn,
  as(cat, 'Spatial'),
  fun      = 'mean',
  weights  = 'area',
  progress = FALSE
) %>%
  t() %>%
  as.data.frame()
```

El método de Hargreaves-Samani requiere información adicional de la cuenca, como la elevación media y la latitud. Estas variables permiten estimar la radiación extraterrestre y ajustar el cálculo de la evapotranspiración potencial a las condiciones geográficas de la zona de estudio.

```{r}
#-----------------------------------------------------------
# Calcular variables fisiográficas para Hargreaves-Samani
#-----------------------------------------------------------

# Calcular la elevación media de la cuenca a partir del DEM
elevation <- exact_extract(
  dem,
  as(cat, 'Spatial'),
  fun      = 'mean',
  weights  = 'area',
  progress = FALSE
) %>%
  t() %>%
  as.data.frame()

# Calcular la latitud del centroide de la cuenca
lat_roi <- crds(centroids(cat))[, 2]

# Función para convertir grados a radianes
deg2rad <- function(deg) {
  (deg * pi) / 180
}

# Convertir latitud de grados a radianes
lat_rad <- deg2rad(lat_roi)
```

Con las series de temperatura y las características geográficas de la cuenca, se organizan los datos y constantes requeridos por la función de Hargreaves-Samani. Este paso permite estructurar las entradas en el formato esperado por la función de cálculo.

```{r}
#-----------------------------------------------------------
# Preparar datos para el cálculo de evapotranspiración potencial
#-----------------------------------------------------------

# Definir constantes requeridas por el método
constants <- list(
  Elev    = elevation,
  lambda  = 2.45,  # Calor latente de vaporización [MJ kg-1]
  lat_rad = lat_rad,
  Gsc     = 0.082  # Constante solar [MJ m-2 min-1]
)

# Organizar datos climáticos diarios
data <- list(
  Date.daily = time,
  J          = as.POSIXlt(time)$yday,
  Tmax       = zoo(dtx, time),
  Tmin       = zoo(dtn, time)
)
```

La evapotranspiración potencial se calcula mediante el método de Hargreaves-Samani, el cual utiliza la temperatura máxima, temperatura mínima y la radiación solar estimada a partir de la ubicación geográfica. El resultado corresponde a una serie diaria de ETP media para la cuenca.

```{r}
#-----------------------------------------------------------
# Calcular evapotranspiración potencial diaria
#-----------------------------------------------------------

# Estimar la evapotranspiración potencial diaria
evap <- ET.HargreavesSamani(
          data            = data,
          constants       = constants,
          ts              = 'daily',
          message         = 'no',
          AdditionalStats = 'no',
          save.csv        = 'no'
        )$ET.Daily

head(evap)
```

## Series de caudales observados

En esta sección se cargan los caudales observados en la estación hidrométrica ubicada a la salida de la cuenca. Esta información será utilizada como serie de referencia para la calibración y evaluación del modelo hidrológico GR4J.

```{r}
#-----------------------------------------------------------
# Cargar caudales observados
#-----------------------------------------------------------

# Leer archivo Excel con datos de caudal observado
raw <- read.xlsx(
  './Data/OBS/data_caudal.xlsx',
  colNames    = FALSE,
  sheet       = 1,
  detectDates = TRUE
)

# Extraer valores de caudal y convertirlos a numérico
obs <- data.frame(
  sapply(raw[-1:-6, -1], as.numeric)
)
```

Luego de cargar los datos observados, se reemplazan los códigos de datos faltantes por valores NA. Dado que el modelo hidrológico trabaja con variables expresadas en unidades de lámina (mm/día), es necesario convertir el caudal observado de m³/s a mm/día. Esta conversión se realiza utilizando el área de la cuenca, permitiendo así trabajar de manera consistente con la precipitación y la evapotranspiración potencial.

```{r}
#-----------------------------------------------------------
# Limpiar y convertir caudal observado
#-----------------------------------------------------------

# Reemplazar código de dato faltante por NA
obs[obs == -999.9] <- NA

# Opcional: remover valores sospechosos si han sido identificados previamente
# obs[obs >= 500] <- NA

# Convertir caudal de m3/s a mm/día
flow <- zoo(obs * 86.4 / area, time)

head(flow)
```

## Exportar datos de entrada

En este paso se integran las variables hidrometeorológicas procesadas en un único conjunto de datos. Este dataset contiene la precipitación, evapotranspiración potencial y caudal observado a escala diaria, y corresponde al formato de entrada requerido para el modelo hidrológico GR4J.

```{r}
#-----------------------------------------------------------
# Construcción del dataset de entrada (formato airGR)
#-----------------------------------------------------------

# Integrar variables en un solo data frame
Inputs <- data.frame(
  DatesR = as.POSIXct(time),  # fechas en formato POSIXct requerido por airGR
  P      = as.numeric(prec),       # precipitación [mm/d]
  E      = as.numeric(evap),       # evapotranspiración potencial [mm/d]
  Q      = as.numeric(flow)        # caudal observado [mm/d]
)

# Visualizar las primeras filas
head(Inputs)
```

```{r}
#-----------------------------------------------------------
# Visualización conjunta de variables hidrometeorológicas
#-----------------------------------------------------------

# Convertir a formato largo
Inputs_long <- Inputs %>%
  tidyr::pivot_longer(
    cols      = all_of(c("P", "E", "Q")),
    names_to  = "variable",
    values_to = "value"
  ) %>%
  dplyr::mutate(
    variable = factor(
      variable,
      levels = c("P", "E", "Q"),
      labels = c(
        "Precipitación [mm/d]",
        "Evapotranspiración potencial [mm/d]",
        "Caudal observado [mm/d]"
      )
    )
  )

# Crear gráfico en tres paneles
ggplot(Inputs_long, aes(x = DatesR, y = value, color = variable)) +
  
  geom_line(linewidth = 0.35) +
  facet_wrap(~ variable, ncol = 1, scales = "free_y") +
  scale_color_manual(values = c(
    "Precipitación [mm/d]" = "dodgerblue",
    "Evapotranspiración potencial [mm/d]" = "darkgreen",
    "Caudal observado [mm/d]" = "black"
  )) +
  labs(
    x = "Tiempo",
    y = NULL,
    title = "Variables hidrometeorológicas de entrada al modelo GR4J"
  ) +
  theme_bw(base_size = 11) +
  theme(
    legend.position = "none",  # quitamos leyenda porque ya está en los paneles
    plot.title = element_text(hjust = 0.5, face = "bold"),
    strip.text = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )
```

```{r}
# Guardar dataset en formato RDS
saveRDS(
  Inputs,
  file = "./Modules/Clase_3/inputs/model_inputs_lumped.rds"
)

# Exportar dataset en formato CSV
write.csv(
  Inputs,
  file = "./Modules/Clase_3/inputs/model_inputs_lumped.csv",
  row.names = FALSE
)
```