# Programa para determinar clusters de series de tiempo

#---- 0. LIBRERÍAS ----
library(factoextra)
library(dtwclust)
library(cluster)
library(ggplot2)
library(ggvoronoi)
library(dplyr)
library(purrr)
library(readr)
library(readxl)
library(tibble)
library(writexl)

rm(list = ls())
gc()
options(scipen = 999)

#---- 1. DATOS ----

# Cargue de las series de tiempo completadas mediante imputación con IDW
rainfall_data <- read_csv("D:/guvelez/Personales/Otros/Scripts/series_completas.csv", 
                          col_types = cols(...1 = col_skip(), 
                                           FECHALECTURA = col_skip(), 
                                           day = col_skip()))
# Transposición de la matriz de datos
rainfall_data <- t(rainfall_data)

# Creación del dataframe para hallar el número óptimo de clusters
rainfall_data_c <- as.data.frame(rainfall_data)
rainfall_data_c <- rownames_to_column(rainfall_data_c, var = "station_id")

# Cargue de las coordenadas de las estaciones pluviométricas
stations <- read_xlsx("D:/guvelez/Personales/Otros/Scripts/Coordenadas_1.xlsx")
names(stations)[1] <- "station_id"

# Combinar datos de precipitación con localización
combined_data <- merge(rainfall_data_c, stations[, c("station_id", "NORTE", "ESTE")], 
                       by = "station_id")
combined_data$station_id <- NULL
combined_data <- as.matrix(combined_data)


#---- 2. DETERMINACIÓN DEL NÚMERO ÓPTIMO DE CLUSTERS CON EL MÉTODO DEL CODO ----

# Calcular el total de la suma de cuadrados en cada cluster para diferentes k
wss <- function(k) {
  kmeans(combined_data, centers = k, nstart = 55)$tot.withinss
}

# Iterar la función sobre valores de 1 a 48
k_values <- 1:55
wss_values <- sapply(k_values, wss)


# Trazar la gráfica del método del codo
fviz_nbclust(combined_data, kmeans, method = "wss", k.max = 55)

# Hallar el valor del "codo" o número óptimo de clusters
# Hallar la segunda derivada del gráfico del codo
wss_diff <- diff(wss_values, differences = 2)

# El número óptimo de clusters corresponde al valor en el que se tiene un máximo de la segunda derivada
elbow_k <- which.max(abs(wss_diff)) + 1 
optimal_k <- elbow_k

#optimal_k <- 3

#---- 3. FORMACIÓN DE CLUSTERS (MÉTODO 1) ----

# Clusterizar los datos de precipitación usando DTW
ts_cluster <- tsclust(rainfall_data, type = "partitional", k = optimal_k, distance = "dtw_basic")

# Agregar los clusters DTW a los datos de localización
stations$cluster <- ts_cluster@cluster

# Usar clustering jerárquico con los clusters DTW y los datos espaciales
d <- dist(stations[, c("NORTE", "ESTE")])  # Distancia basada en la información espacial
hc <- hclust(d)

# "Cortar el árbol" para obtener los clusters
locations$spatial_cluster <- cutree(hc, k = optimal_k)

#---- 4. FORMACIÓN DE CLUSTERS (MÉTODO 2) ----

set.seed(123)
kmeans_result <- kmeans(cbind(rainfall_data, locations[, c("latitude", "longitude")]), centers = optimal_k)
locations$cluster_2 <- kmeans_result$cluster  # Assign clusters to the stations


#---- 5. FORMACIÓN DE CLUSTERS (MÉTODO 3: COMBINAR INFORMACIÓN ESPACIAL Y DE SERIES DE TIEMPO) ----

# Normalización de los datos espaciales (coordenadas este y norte)
normalize_sd <- function(x) { (x - min(x)) / (max(x) - min(x)) }
stations$NORTE <- normalize_sd(stations$NORTE)
stations$ESTE <- normalize_sd(stations$ESTE)

# Cálculo de la matriz de distancias espaciales (distancias euclidianas)
spatial_dist <- dist(stations[, c("NORTE", "ESTE")])

# Normalización de los datos de cada estación 
normalize_ts <- function(x) { (x - min(x)) / (max(x) - min(x)) }
rainfall_data_normalized <- apply(rainfall_data, 1, normalize_ts)

# Transposición de la matriz de datos de series normalizada
rainfall_data_normalized <- t(rainfall_data_normalized)

# Cálculo de la matriz de distancias DTW para las series de precipitación
dtw_dist <- proxy::dist(rainfall_data_normalized, method = "dtw_basic")

# Combinar las distancias espaciales y DTW mediante una suma ponderada
alpha <- 0.5  # Peso para la distancia espacial
combined_dist <- (1 - alpha) * dtw_dist + alpha * as.matrix(spatial_dist)

# Aplicar clusterización k-means a las series de tiempo, usando el k óptimo del 
# método del codo
ts_cluster <- tsclust(combined_dist, type = "partitional", k = optimal_k, distance = "dtw_basic", centroid = "pam")

# Obtener la asignación de cluster para cada estación
stations$cluster <- ts_cluster@cluster


#---- 6. VISUALIZAR LOS CLUSTERS ----

# Graficar las estaciones por cluster
ggplot(stations, aes(x = ESTE, y = NORTE, color = factor(cluster))) +
  geom_point(size = 3) +
  labs(title = "Clusterización de las Estaciones Pluviométricas por Precipitación y Localización",
       x = "Este", y = "Norte") +
  theme_minimal()

# Plot the Voronoi polygons based on the station locations
ggplot() +
  geom_voronoi(aes(x = ESTE, y = NORTE, fill = factor(cluster), group = cluster), data = stations, alpha = 0.3, color = "black") +  
  geom_point(data = stations, aes(x = ESTE, y = NORTE, color = factor(cluster)), size = 3) +  
  labs(title = "Clusterización de las Estaciones Pluviométricas por Precipitación y Localización",
       x = "ESTE", y = "NORTE") +
  theme_minimal()

#---- 7. VISUALIZAR LAS SERIES DE TIEMPO ----


# Dataframe con los datos de precipitación
time_series_df <- data.frame(
  Day = rep(1:4017, times = 56),              # Repite los días por cada estación
  Rainfall = as.vector(rainfall_data),         # "Aplana" la matriz de datos 
  Station = factor(rep(stations$station_id, each = 4017)),      # IDs de las estaciones
  Cluster = factor(rep(stations$cluster, each = 4017))  # Asignación de cluster
)


# Cálculo de la serie de tiempo promedio para cada cluster
avg_time_series <- aggregate(Rainfall ~ Day + Cluster, data = time_series_df, FUN = mean)

# Gráfica de la serie de tiempo promedio para cada cluster
ggplot(avg_time_series, aes(x = Day, y = Rainfall, color = Cluster, group = Cluster)) +
  geom_line(linewidth = 1.2) +
  labs(title = "Average Time Series of Rainfall by Cluster", x = "Day", y = "Average Rainfall (mm)") +
  theme_minimal()

write_csv(stations, "estaciones_cluster_diario.csv")

# Exportar las series de tiempo promedio para cada cluster
avg_rainfall_series_wide <- avg_time_series %>%
  pivot_wider(names_from = Cluster, values_from = Rainfall)

avg_rainfall_series_wide$FECHALECTURA <- seq(as.POSIXct("2013-01-01", tz = "UTC"), as.POSIXct("2023-12-31", tz = "UTC"), by = "day")
avg_rainfall_series_wide$Day <- NULL
names(avg_rainfall_series_wide) <- c("PRECIPITACION_CL1", "PRECIPITACION_CL2", "FECHALECTURA")

avg_rainfall_series_wide <- avg_rainfall_series_wide %>%
  select(FECHALECTURA, everything())

for(i in 1:2){
  serie <- avg_rainfall_series_wide[,c(1,i+1)]
  ruta <- paste0("D:/guvelez/Personales/Otros/Scripts/DATOS_PROMEDIO_CLUSTERS/CLUSTER_", i, ".xlsx")
  write_xlsx(serie, ruta)
}
