semana4

# ----------------------------------------------------------
# PASO 0: SELECCIÓN DE AÑO PARA ANÁLISIS
# ----------------------------------------------------------
año_seleccionado <- 2006

# ----------------------------------------------------------
# PASO 1: INTEGRACIÓN DE DATOS DESDE ARCHIVOS TXT
# ----------------------------------------------------------

# Lista de estaciones meteorológicas
estaciones <- c("Puno", "Isla Taquile", "Capachica", "Azangaro", "Arapa",
                "Progreso", "Muñani", "Crucero", "Juli", "Desaguadero",
                "Pizacoma", "Ilave", "Capazo", "Mazo Cruz", "Huancane","Cojata",
                "Lampa","Putina","Huaraya Moho","Isla Suana","Tahuaco Yunguyo",
                "Cuyo Cuyo","Cabanillas","Ananea","Pucara","Pampahuta","Chuquibambilla",
                "Santa Rosa","Ayaviri","Crucero Alto")

# Definir nombres de meses en español para conversión posterior
meses_español <- c("Ene", "Feb", "Mar", "Abr", "May", "Jun", 
                  "Jul", "Ago", "Sep", "Oct", "Nov", "Dic")

# Cargar librerías necesarias

library(dplyr)    # Para manipulación de datos
library(tidyr)    # Para reformatear datos
library(purrr)    # Para funciones de iteración
library(fda)


# Leer y combinar todos los archivos

datos_completos <- bind_rows(lapply(estaciones, function(est) {
  # Leer cada archivo de texto
  datos_estacion <- read.table(
    file = paste0(est, ".txt"),  # Construye nombre del archivo
    header = FALSE,              # Los archivos no tienen encabezado
    col.names = c("AÑO", "MES", "DIA", "PRECIPITACION", "TEMP_MAX", "TEMP_MIN")  # Asignar nombres de columnas
  )
  
  # Limpieza y transformación de datos
  datos_estacion %>%
    mutate(
      Estacion = est,  # Añadir columna con nombre de estación
      # Reemplazar valores faltantes codificados como -99.9
      PRECIPITACION = ifelse(PRECIPITACION == -99.9, NA, PRECIPITACION),
      TEMP_MAX = ifelse(TEMP_MAX == -99.9, NA, TEMP_MAX),
      TEMP_MIN = ifelse(TEMP_MIN == -99.9, NA, TEMP_MIN),
      # Manejar posibles valores infinitos
      across(c(PRECIPITACION, TEMP_MAX, TEMP_MIN), 
             ~ifelse(is.infinite(.), NA, .))
    )
})) %>%
  filter(AÑO == año_seleccionado)  # FILTRADO POR AÑO SELECCIONADO

# ----------------------------------------------------------
# PASO 2: PREPARACIÓN DE DATOS PARA ANÁLISIS FUNCIONAL
# ----------------------------------------------------------


datos_agregados <- datos_completos %>%
  mutate(
    NOMBRE_MES = meses_español[MES],  # Columna con nombres en español
    MES_NUM = MES                      # Conservar número original (1-12)
  ) %>%
  group_by(Estacion, MES_NUM, NOMBRE_MES) %>%  # Agrupar por estación y mes
  summarise(
    # Calcular promedios mensuales para temperaturas
    TEMP_MAX = mean(TEMP_MAX, na.rm = TRUE),
    TEMP_MIN = mean(TEMP_MIN, na.rm = TRUE),
    # Calcular precipitación acumulada mensual
    PRECIPITACION = sum(PRECIPITACION, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  rename(MES = MES_NUM)  # Renombrar para claridad


# Reformatear datos a formato ancho para análisis funcional

datos_ancho <- datos_agregados %>%
  pivot_wider(
    names_from = Estacion,          # Crear columnas por estación
    values_from = c(TEMP_MAX, TEMP_MIN, PRECIPITACION)  # Variables a expandir
  )

# ----------------------------------------------------------
# PASO 3: CREAR MATRICES PARA ANÁLISIS FUNCIONAL (VERSIÓN SEGURA)
# ----------------------------------------------------------

# Función para filtrar columnas con datos completos (sin NAs)
filtrar_datos_completos <- function(matriz) {
  matriz[, colSums(is.na(matriz)) == 0, drop = FALSE]  # Columnas sin NAs
}

# 3.1 Matriz para Temperatura Máxima 

temp_max_matrix <- datos_ancho %>%
  dplyr::select(MES, dplyr::starts_with("TEMP_MAX_")) %>% 
  arrange(MES) %>%
  tibble::column_to_rownames("MES") %>%
  as.matrix() %>%
  filtrar_datos_completos()

# 3.2 Matriz para Temperatura Mínima 

temp_min_matrix <- datos_ancho %>%
  dplyr::select(MES, dplyr::starts_with("TEMP_MIN_")) %>%
  arrange(MES) %>%
  tibble::column_to_rownames("MES") %>%
  as.matrix() %>%
  filtrar_datos_completos()

# 3.3 Matriz para Precipitación 

precip_matrix <- datos_ancho %>%
  dplyr::select(MES, dplyr::starts_with("PRECIPITACION_")) %>%
  arrange(MES) %>%
  tibble::column_to_rownames("MES") %>%
  as.matrix() %>%
  filtrar_datos_completos()

# ----------------------------------------------------------
# PASO 4: ANÁLISIS FUNCIONAL PARA TEMPERATURA MÁXIMA
# ----------------------------------------------------------

# 4.1 Configuración base Fourier (ajustada a 12 meses)
base_fourier_temp <- create.fourier.basis(rangeval = c(1, 12), nbasis = 5)

# 4.2 Conversión a objeto funcional
fd_temp_max <- Data2fd(argvals = 1:12, y = temp_max_matrix, basisobj = base_fourier_temp)

# 4.3 Visualizacion
par(mar = c(5, 4, 4, 7), xpd = TRUE, mgp = c(2, 0.5, 0))
plot(fd_temp_max, 
     main = "Temperatura Máxima por Estación",
     xlab = "Mes", 
     ylab = "°C",
     axes = FALSE)

## [1] "done"

# Ejes personalizados
axis(1, at = 1:12, labels = meses_español, las = 2, cex.axis = 0.8)
axis(2, las=1)
box()
grid(nx = NA, ny = NULL, col = "gray" , lty=3)
legend("topright", inset = c(-0.25, 0), legend = estaciones,
       col = rainbow(length(estaciones)), lty = 1, pch = 19, seg.len = 0.8,
       cex = 0.6, xpd = TRUE, bty = "n",  y.intersp = 0.8)

# 4.4 Análisis de tendencia

media_temp_max <- mean.fd(fd_temp_max)
plot(media_temp_max, lwd = 2, col = "red", main = "Temperatura Máxima Media")

## [1] "done"

# ----------------------------------------------------------
# PASO 5: ANÁLISIS FUNCIONAL PARA TEMPERATURA MÍNIMA
# ----------------------------------------------------------


# 5.1 Usamos la misma base Fourier
fd_temp_min <- Data2fd(argvals = 1:12, y = temp_min_matrix, basisobj = base_fourier_temp)

# 5.2 Visualización específica
plot(fd_temp_min, main = "Temperatura Mínima por Estación",
     xlab = "Mes", ylab = "°C")

## [1] "done"

# 5.3 Análisis de tendencia
media_temp_min <- mean.fd(fd_temp_min)
plot(media_temp_min, lwd = 2, col = "blue", main = "Temperatura Mínima Media")

## [1] "done"

# -----------------------------------------------
# PASO 6: ANÁLISIS FUNCIONAL PARA PRECIPITACIÓN 
# -----------------------------------------------


# 6.1 Base diferente para precipitación
base_fourier_precip <- create.fourier.basis(rangeval = c(1, 12), nbasis = 5)

# 6.2 Conversión a objeto funcional
fd_precip <- Data2fd(argvals = 1:12, y = precip_matrix, basisobj = base_fourier_precip)

# 6.3 Visualización
plot(fd_precip, main = "Precipitación Mensual", xlab = "Mes", ylab = "mm")

## [1] "done"