---
title: "Regionalización Hidroclimática de Cuencas en el Perú"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Regionalización Hidroclimática de Cuencas en el Perú}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 6,
  out.width = "100%"
)
```

La **regionalización hidrológica** es un procedimiento fundamental en la hidrología científica y aplicada. Su objetivo principal es clasificar cuencas hidrográficas en regiones homogéneas que comparten características físicas y climáticas similares. Esta técnica es de gran utilidad para:

1. **Estimación en cuencas no instrumentadas**: Transferir parámetros de modelos o firmas hidrológicas desde cuencas con datos (donantes) hacia cuencas sin mediciones de caudal (receptoras).
2. **Comprensión del régimen de caudales**: Identificar los principales factores físicos (pendiente, elevación) y climáticos (precipitación, aridez) que gobiernan la respuesta hidrológica regional.

Este tutorial muestra cómo utilizar el paquete `rcamelspe` para cargar atributos climáticos y topográficos, clasificar las 136 cuencas del Perú utilizando un algoritmo de agrupamiento (*K-means*), y mapear las regiones resultantes.

---

## 1. Carga de Librerías y Configuración

Primero, cargamos el paquete `rcamelspe` junto con `ggplot2` y `sf` para el manejo espacial y la visualización.

```{r setup, message=FALSE, warning=FALSE}
if (!requireNamespace("rcamelspe", quietly = TRUE)) {
  pkg_dir <- if (file.exists("../DESCRIPTION")) ".." else if (file.exists("DESCRIPTION")) "." else NULL
  if (!is.null(pkg_dir)) {
    r_files <- list.files(file.path(pkg_dir, "R"), full.names = TRUE, pattern = "\\.[rR]$")
    invisible(lapply(r_files, sys.source, envir = globalenv()))
  }
}
if (requireNamespace("rcamelspe", quietly = TRUE)) {
  library(rcamelspe)
}
library(ggplot2)
library(sf)

old_camels_path <- getOption("rcamelspe.path")

# Configurar la ruta al dataset CAMELS-PE si está disponible localmente
# Buscamos rutas relativas comunes en el espacio de trabajo
path <- "raw-data/CAMELS-PE"
if (dir.exists(path) && exists("set_camels_pe_path", mode = "function")) {
  set_camels_pe_path(path)
} else if (dir.exists("../raw-data/CAMELS-PE") && exists("set_camels_pe_path", mode = "function")) {
  set_camels_pe_path("../raw-data/CAMELS-PE")
} else if (dir.exists("../../raw-data/CAMELS-PE") && exists("set_camels_pe_path", mode = "function")) {
  set_camels_pe_path("../../raw-data/CAMELS-PE")
}

# Verificar si el dataset está disponible antes de ejecutar
has_data <- exists("get_camels_pe_path", mode = "function") &&
  !is.null(get_camels_pe_path()) &&
  dir.exists(get_camels_pe_path())
```

---

## 2. Preparación y Carga de Datos

Cargamos de forma eficiente los atributos físicos de interés y los límites geoespaciales de las cuencas. Utilizaremos variables climáticas y topográficas clave:

* **Física/Topografía**: Área de la cuenca (`area`), elevación media (`elev_mean`) y pendiente media (`slope_mean`).
* **Clima**: Precipitación media diaria (`p_mean`) y el índice de aridez (`aridity`).

```{r load_data, message=FALSE}
# Verificar si el dataset está disponible antes de ejecutar
has_data <- exists("get_camels_pe_path", mode = "function") &&
  !is.null(get_camels_pe_path()) &&
  dir.exists(get_camels_pe_path())

if (has_data) {
  # Cargar atributos topográficos y climáticos mezclados por gauge_id
  attrs <- load_pe_attributes(attributes = c("topographic", "climatic"))
  
  # Cargar geometrías vectoriales de las cuencas (polígonos GPKG)
  catchments <- load_pe_geospatial(type = "catchments")
  
  # Seleccionar las variables numéricas para el análisis
  clustering_vars <- c("area", "elev_mean", "slope_mean", "p_mean", "aridity")
  
  # Filtrar registros completos
  data_clustering <- attrs[, c("gauge_id", clustering_vars)]
  data_clustering <- na.omit(data_clustering)
  
  # Vista rápida de los atributos cargados
  head(data_clustering)
}
```

---

## 3. Preprocesamiento y Agrupamiento (Clustering K-means)

Dado que las variables físicas y climáticas tienen unidades y rangos drásticamente diferentes (por ejemplo, el área en $km^2$ vs el índice de aridez que es adimensional), debemos **escalar** (estandarizar) los datos para que tengan media igual a 0 y desviación estándar de 1.

Posteriormente, aplicamos el algoritmo clásico de particionamiento **K-means** seleccionando $K = 3$ grupos para identificar tres macrorregiones hidroclimáticas representativas en el Perú.

```{r clustering}
if (has_data) {
  # Estandarizar atributos
  scaled_data <- scale(data_clustering[, clustering_vars])
  
  # Ejecutar algoritmo K-means si hay suficientes cuencas, o partición adaptativa
  set.seed(123) # Para reproducibilidad
  if (nrow(data_clustering) >= 3) {
    k_centers <- min(3, nrow(data_clustering))
    kmeans_res <- kmeans(scaled_data, centers = k_centers, nstart = 25)
    data_clustering$region <- as.factor(kmeans_res$cluster)
  } else {
    data_clustering$region <- as.factor(seq_len(nrow(data_clustering)))
  }
  
  # Contar cuántas cuencas pertenecen a cada región
  table(data_clustering$region)
}
```

---

## 4. Unión Geoespacial de Resultados

Unimos nuestros resultados de regionalización con los límites geográficos de las cuencas. Aprovechamos la función de alto rendimiento `collapse::join()` para acelerar el proceso y nos aseguramos de preservar la clase de objeto geoespacial `sf`.

```{r geospatial_join}
if (has_data) {
  # Realizar un Left Join rápido entre el objeto espacial y los resultados de clustering
  catchments_regions <- collapse::join(
    catchments, 
    data_clustering[, c("gauge_id", "region")], 
    on = "gauge_id", 
    how = "left"
  )
  
  # Forzar a sf si collapse limpió los atributos espaciales
  if (!inherits(catchments_regions, "sf")) {
    catchments_regions <- sf::st_as_sf(catchments_regions)
  }
}
```

---

## 5. Visualización de las Regiones Hidroclimáticas

Finalmente, graficamos el mapa temático de las cuencas hidrográficas coloreándolas de acuerdo a la región a la que fueron asignadas por el algoritmo.

```{r plot_map, fig.align='center', eval = isTRUE(has_data)}

  ggplot(catchments_regions) +
    geom_sf(aes(fill = region), color = "grey40", linewidth = 0.15) +
    scale_fill_brewer(palette = "Set2", name = "Región Hidroclimática\n(Grupo)") +
    labs(
      title = "Regionalización Hidroclimática de Cuencas en el Perú",
      subtitle = "Clasificación basada en K-means (K = 3) sobre atributos físicos y climáticos",
      caption = "Datos físicos y climáticos extraídos de CAMELS-PE",
      x = NULL,
      y = NULL
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(face = "bold", size = 12),
      plot.subtitle = element_text(color = "grey30", size = 9),
      legend.position = "right",
      panel.grid.major = element_line(color = "grey90", linewidth = 0.2)
    )

```

---

## 6. Interpretación de las Regiones Hidroclimáticas

Para entender qué representa cada clúster, podemos resumir el promedio de los atributos físicos y climáticos para cada grupo identificado:

```{r interpretation}
if (has_data) {
  # Calcular promedios grupales utilizando collapse::collap
  resumen_regiones <- collapse::collap(
    data_clustering[, c("region", clustering_vars)],
    by = ~region,
    FUN = mean
  )
  
  print(resumen_regiones)
}
```

* **Región 1**: Cuencas con menor precipitación media y alta aridez (típicamente asociadas a la vertiente árida del Pacífico).
* **Región 2**: Cuencas de gran altitud y pendientes pronunciadas (zonas andinas de cabecera de cuenca).
* **Región 3**: Cuencas con altas precipitaciones y menor aridez (regiones húmedas de la vertiente del Atlántico/Amazonía).

```{r cleanup, include=FALSE}
options(rcamelspe.path = old_camels_path)
```