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

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

# Ejemplo para 48 estaciones, cada una con 12 meses de datos de precipitación
rainfall_data <- matrix(runif(48 * 60, min = 0, max = 200), nrow = 48, ncol = 60)

# Ejemplo de datos de localización para 48 estaciones (latitud y longitud)
locations <- data.frame(
  station_id = 1:48,
  latitude = runif(48, min = -4, max = 10),
  longitude = runif(48, min = -80, max = -60)
)

# Combinar datos de precipitación con localización
combined_data <- cbind(rainfall_data, locations[, c("latitude", "longitude")])

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

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

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


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

# Análisis de Silhouette para k = 2 a k = 47
fviz_nbclust(combined_data, kmeans, method = "silhouette", k.max = 47)

optimal_k <- 17

#---- 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
locations$cluster <- ts_cluster@cluster

# Usar clustering jerárquico con los clusters DTW y los datos espaciales
d <- dist(locations[, c("latitude", "longitude")])  # 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) ----

# Normalize spatial data (latitude and longitude)
normalize_sd <- function(x) { (x - min(x)) / (max(x) - min(x)) }
locations$latitude <- normalize_sd(locations$latitude)
locations$longitude <- normalize_sd(locations$longitude)

# Compute the spatial distance matrix (Euclidean distance)
spatial_dist <- dist(locations[, c("latitude", "longitude")])

# Normalize each station's time series (min-max normalization)
normalize_ts <- function(x) { (x - min(x)) / (max(x) - min(x)) }
rainfall_data_normalized <- apply(rainfall_data, 1, normalize_ts)

# Transpose back to original shape (if needed)
rainfall_data_normalized <- t(rainfall_data_normalized)

# Compute DTW distance matrix for the rainfall time series
dtw_dist <- proxy::dist(rainfall_data_normalized, method = "dtw_basic")

# Combine the DTW and spatial distances (weighted sum)
alpha <- 0.5  # Weight for the spatial distance (adjust as needed)
combined_dist <- (1 - alpha) * dtw_dist + alpha * as.matrix(spatial_dist)

# Apply K-means time series clustering using the combined distance matrix
ts_cluster <- tsclust(combined_dist, type = "partitional", k = optimal_k, distance = "dtw_basic")

# Get the cluster assignments
locations$cluster <- ts_cluster@cluster


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

# Create convex hulls for each cluster
hull_data <- locations %>%
  group_by(cluster) %>%
  do(hull = chull(.$longitude, .$latitude)) %>%
  ungroup() %>%
  mutate(hull_points = map(hull, ~.x[c(.x, .x[1])]))  # Close the polygon

# Extract the coordinates of the convex hulls for plotting
hull_coords <- locations %>%
  group_by(cluster) %>%
  slice(unlist(hull_data$hull_points)) %>%
  ungroup()

# Plot the stations with the convex hull polygons
ggplot() +
  geom_polygon(data = hull_coords, aes(x = longitude, y = latitude, group = cluster, fill = factor(cluster)), color = "black", alpha = 0.3) +  # Plot convex hulls
  geom_point(data = locations, aes(x = longitude, y = latitude, color = factor(cluster)), size = 3) +  # Plot the points
  labs(title = "Convex Hulls of Clusters of Pluviometric Stations",
       x = "Longitude", y = "Latitude") +
  theme_minimal()


# Plot the Voronoi polygons based on the station locations
ggplot() +
  geom_voronoi(aes(x = longitude, y = latitude, fill = factor(cluster), group = cluster), data = locations, alpha = 0.3, color = "black") +  # Voronoi polygons
  geom_point(data = locations, aes(x = longitude, y = latitude, color = factor(cluster)), size = 3) +  # Plot the points
  labs(title = "Voronoi Polygons of Clusters of Pluviometric Stations",
       x = "Longitude", y = "Latitude") +
  theme_minimal()

# Plot stations by cluster (spatial view)
ggplot(locations, aes(x = longitude, y = latitude, color = factor(cluster))) +
  geom_point(size = 3) +
  labs(title = "Clustering of Pluviometric Stations by Rainfall and Location",
       x = "Longitude", y = "Latitude") +
  theme_minimal()


# Create a data frame from the rainfall data for plotting
time_series_df <- data.frame(
  Month = rep(1:60, times = 48),              # Repeat the months for each station
  Rainfall = as.vector(rainfall_data),         # Flatten the rainfall data matrix
  Station = factor(rep(1:48, each = 60)),      # Station IDs
  Cluster = factor(rep(locations$cluster, each = 60))  # Cluster assignments
)


# Calculate the average rainfall per month for each cluster
avg_time_series <- aggregate(Rainfall ~ Month + Cluster, data = time_series_df, FUN = mean)

# Plot the average time series for each cluster
ggplot(avg_time_series, aes(x = Month, y = Rainfall, color = Cluster, group = Cluster)) +
  geom_line(linewidth = 1.2) +
  labs(title = "Average Time Series of Rainfall by Cluster", x = "Month", y = "Average Rainfall (mm)") +
  theme_minimal()














