# Load the necessary library
library(terra)

# Load the raster files
raster1 <- rast("c:/Users/oterritorial/OneDrive - Universidad Nacional de Chimborazo/_00UNACH/Clases2024_02S/06_Ordenamiento Territorial/PAE/UH Relieve/volcan/Dem_Chimborazo.tif")
raster2 <- rast("c:/Users/oterritorial/OneDrive - Universidad Nacional de Chimborazo/_00UNACH/Clases2024_02S/06_Ordenamiento Territorial/PAE/UH Relieve/volcan/roughness.tif")
raster3 <- rast("c:/Users/oterritorial/OneDrive - Universidad Nacional de Chimborazo/_00UNACH/Clases2024_02S/06_Ordenamiento Territorial/PAE/UH Relieve/volcan/slope.tif")

# Stack the rasters for easier handling (optional)
raster_stack <- c(raster1, raster2, raster3)

# Mask out any NA values across all rasters
raster_stack <- mask(raster_stack, mask = any(is.na(raster_stack)))

# Extract raster values as a data frame
raster_values <- as.data.frame(values(raster_stack))

# Calculate pairwise correlation coefficients
correlation_matrix <- cor(raster_values, use = "complete.obs")

# Display the correlation matrix
print(correlation_matrix)

# Remove rows with NA or Inf values
valid_values <- complete.cases(raster_values)
raster_values_clean <- raster_values[valid_values, ]

# Perform PCA
pca_result <- prcomp(raster_values_clean, center = TRUE, scale. = TRUE)

# View PCA summary
summary(pca_result)

# Extract the PCA loadings
loadings <- pca_result$rotation
print(loadings)


# Project PCA components back onto rasters (optional)
pca_rasters <- list()
valid_index <- which(valid_values) # Indices of valid cells
for (i in seq_along(pca_result$x[1, ])) {
  # Create raster for each PCA component
  pca_layer <- rast(raster_stack[[1]])
  pca_values <- rep(NA, ncell(pca_layer)) # Initialize with NA
  pca_values[valid_index] <- pca_result$x[, i] # Assign PCA values
  values(pca_layer) <- pca_values
  pca_rasters[[i]] <- pca_layer
  names(pca_rasters[[i]]) <- paste0("PCA_", i)
}
# Combine PCA rasters into a stack
pca_stack <- rast(pca_rasters)

# Save PCA rasters (optional)
writeRaster(pca_stack, "c:/Users/oterritorial/OneDrive - Universidad Nacional de Chimborazo/_00UNACH/Clases2024_02S/06_Ordenamiento Territorial/PAE/UH Relieve/volcan/pca_stack_cleaned.tif", overwrite = TRUE)

# Plot the first PCA component
plot(pca_stack[[1]], main = "PCA Component 3")

############### Análisis de Clusters k-means ######

# Convertir el raster stack a una matriz (valores en filas y bandas en columnas)
pca_values <- as.data.frame(values(pca_stack))

# Filtrar valores válidos (sin NA)
valid_pixels <- complete.cases(pca_values)
pca_values_clean <- pca_values[valid_pixels, ]

# Aplicar k-means clustering (5 clusters)
set.seed(123) # Para reproducibilidad
kmeans_result <- kmeans(pca_values_clean, centers = 5, nstart = 25)

# Crear un raster con los resultados del clustering
cluster_raster <- rast(pca_stack[[1]]) # Crear un raster vacío basado en el primero
cluster_values <- rep(NA, ncell(cluster_raster)) # Inicializar con NA
cluster_values[which(valid_pixels)] <- kmeans_result$cluster # Asignar resultados del clustering
values(cluster_raster) <- cluster_values

# Renombrar el raster
names(cluster_raster) <- "Cluster"

# Visualizar el raster clasificado
plot(cluster_raster, main = "Clasificación en 5 Clusters (K-means)")

# Guardar el resultado en un archivo (opcional)
writeRaster(cluster_raster, "c:/Users/oterritorial/OneDrive - Universidad Nacional de Chimborazo/_00UNACH/Clases2024_02S/06_Ordenamiento Territorial/PAE/UH Relieve/volcan/pca_clusters.tif", overwrite = TRUE)