###### CARGUE DE SERIES DE PRECIPITACIÓN E IMPUTACIÓN DE VALORES FALTANTES CON IDW

#---- 1. LIBRERÍAS ----
library(data.table)
library(dplyr)
library(gstat)
library(tidyr)
library(sp)
library(stringr)
library(writexl)
library(readxl)

#---- 2. LIMPIEZA DE MEMORIA ----
rm(list = ls())
gc()
options(scipen = 999)

#---- 3. CARGUE DE SERIES DE TIEMPO Y AJUSTE DE FECHAS FALTANTES ----

# Estaciones con coordenadas geográficas y planas
stations <- read_xlsx("D:/guvelez/Personales/Otros/Scripts/Coordenadas_1.xlsx")
names(stations)[1] <- "station_id"

# Ruta con los excel de las series de tiempo de precipitación
ruta <- "D:/guvelez/Personales/Otros/Scripts/DATOS"
all_files <- list.files(path = ruta, pattern = "\\.xlsx$", full.names = TRUE)

# Dataframe con todas las fechas por día, desde el 2013-01-01 hasta el 2023-12-31
all_dates <- data.frame(date = seq(as.POSIXct("2013-01-01", tz = "UTC"), as.POSIXct("2023-12-31", tz = "UTC"), by = "day"))
#all_dates <- data.frame(date = seq(as.Date("2013-01-01"), as.Date("2023-12-31"), by = "day"))

# Lectura de bases de datos de precipitación, completando las fechas faltantes con NA
rainfall_series <- data.frame()
i <- 1
for(file in all_files){
  df <- read_xlsx(file)
  df1 <- df
  df$FECHALECTURA <- as.Date(df$FECHALECTURA, format = "%d/%m/%Y")
  # Eliminar duplicados por fecha
  df <- df[!duplicated(df$FECHALECTURA) & !(is.na(df$FECHALECTURA) & duplicated(df$FECHALECTURA)), ]
  
  # Convertir a fecha la columna FECHALECTURA
  #df <- convertir_a_fecha(df)
  
  # Completar las fechas faltantes en cada serie
  df<- df %>%
    complete(FECHALECTURA = all_dates$date) %>%  
    arrange(FECHALECTURA)
  
  if(i == 1) {
    rainfall_series <- df
    names(rainfall_series)[2] <- paste0("E",i)
  } else {
    names(df)[2] <- paste0("E",i)
    rainfall_series <- merge(rainfall_series, df, by = "FECHALECTURA")
  }
  
  i <- i + 1
  
}

#---- 4. IMPUTACIÓN DE VALORES FALTANTES POR IDW ----

# Tabla de datos en formato "largo"
rainfall_long <- rainfall_series %>%
  mutate(day = 1:4017) %>%
  tidyr::pivot_longer(cols = starts_with("E"), names_to = "station_id", values_to = "rainfall")

# Vínculo con las coordenadas de las estaciones
rainfall_long <- left_join(rainfall_long, stations, by = "station_id")

# Función que realiza la imputación de datos faltantes por IDW
# Se fija un total de 4 vecinos (criterio "torre" de vecindad)
idw_fill_missing <- function(rainfall_data, max_neighbors = 4) {
  
  # Conversión a objeto espacial
  coordinates(rainfall_data) <- ~ ESTE + NORTE
  
  # Inicializar una lista para guardar los resultados
  filled_data <- rainfall_data
  
  # Iterar sobre cada día para llenar los valores faltantes empleando IDW
  for (day in unique(rainfall_data$day)) {
    day_data <- rainfall_data[rainfall_data$day == day, ]
    
    # Dataframe de valores nulos
    missing_data <- day_data[is.na(day_data$rainfall), ]
    
    if (nrow(missing_data) > 0) {
      # Generación de un modelo IDW empleando los datos disponibles (no nulos) de precipitación
      idw_model <- gstat(formula = rainfall ~ 1, 
                         data = day_data[!is.na(day_data$rainfall), ], 
                         nmax = max_neighbors,
                         set = list(idp = 2)) # Potencia beta (de la ecuación IDW)
      
      # Aplicación de interpolación IDW para calcular los valores nulos
      idw_results <- predict(idw_model, newdata = missing_data)
      
      # Llenado de valores nulos
      filled_data$rainfall[filled_data$day == day & is.na(filled_data$rainfall)] <- idw_results@data$var1.pred
    }
  }
  return(filled_data)
}

# Se recupera la tabla con los valores completos
a <- idw_fill_missing(rainfall_long)
rainfall_long_filled <- a@data[,c(1:4)]

# Se cambia la tabla a formato ancho (con los datos de cada estación por columna)
rainfall_series_comp <- rainfall_long_filled %>%
  pivot_wider(names_from = station_id, values_from = rainfall)

# Se guarda la tabla resultante como CSV
write.csv(rainfall_series_comp, "D:/guvelez/Personales/Otros/Scripts/series_completas.csv")

# Se guardan las series completas como archivos separados de Excel
for(i in 1:56){
  serie <- rainfall_series_comp[,c(1,i+2)]
  names(serie)[2] <- "PRECIPITACION"
  ruta <- str_replace(all_files[i], "DATOS", "DATOS_COMPLETOS")
  write_xlsx(serie, ruta)
}


# PENDIENTE. ¿Cómo determinar el error y la bondad de ajuste de esta imputación?

