Este cuaderno es una guía para facilitar el desarrollo de la actividades en clase de los cursos: Hidroclimatológia (pregrado y posgrado) y Modelación de procesos hidrológicos de la FEAR-PUJ, y fue consolidado con la ayuda de Diego Cortés (estudiante de Maestría en Recursos Hídricos en la Unal).

El presente cuaderno muestra una introducción al procesamiento de series de datos hidrometeorológicos mediante el uso de R. Se muestra una comparación de datos de precipitación y caudal con respecto a los índices ENOS.

Antes de iniciar, es recomendable iniciar limpiando el espacio de trabajo y llamando de una vez todos las librerías que se vayan a necesitar. Si no se tiene instalada alguna de las librerías es necesario hacerlo con el comando “install.packages”. Así mismo, se selecciona el directorio de trabajo donde se tiene el archivo descargado del DHIME.

rm(list = ls())

require(xts); # Series Temporales
require(lubridate); # Utilidades para manejo de fechas
require(dplyr); # To work with data frame like objects
require(rsoi); #Import Various Northern and Southern Hemisphere Climate Indices https://cran.r-project.org/web/packages/rsoi/index.html
require(ggplot2); # Plot
require(ggpubr) # Plot

Una vez instalados los paquetes requeridos, por favor indique con el comando “setwd” la ubicación de la CARPETA donde tiene los archivos con los que vamos a trabajar.Por favor recuerde que para correr este comando debe usar DOBLE SLASH. Es recomendable que cree una carpeta de trabajo para este curso y que esta no este guardada en la nube (OneDrive) para evitar problemas de sincronización por falta de internet durante las clases.

# Carpeta de trabajo
setwd("D://Javeriana//Scripts - R")

Iniciamos cargando los datos que se analizaran. Estos deben corresponder a precipitación y caudal. En el ejemplo que se desarrollará a continuación usted tiene la posibilidad de realizar los análisis y gráficas de comparación entre las variables y ONI para varias estaciones en la misma variable. Si usted desea realizar el análisis para UNA ÚNICA ESTACIÓN por favor modifique el código en las funciones for y podría considerar usar la función read.csv en lugar de list.files.

fls <- list.files("Registros Filtro 1/", ".csv$", full.names = T) # Lectura y Carga de Archivos con Registros
vrs <- c("PTPM", "Q_M") # Variables a analizar

Una vez cargados los archivos, se asignan los nombres de las variables a analizar a los archivos cargados, con el fin de asegurarnos que todas las variables auqnue correspondan a diferentes estaciones tengan el mismo nombre. Adicionalmente, descargamos los registros del ONI para iniciar el proceso de comparación.

fls <- fls[ sapply(vrs, function(x) grepl(x, fls)) %>% apply( 1, any) ] # se cambian los nombres por si los archivos no se llaman igual
oni <- download_oni() # Indice Oceanic Nino Index
oni <- as.xts(oni[,-3],
        order.by = as.POSIXct(oni$Date %>% as.character, tz = "GMT")) #en este caso usamos como zona horaria GMT

Ahora usaremos la función for para realizar el análisis a varias estaciones las cuales estan contenidas en el archivo incial. Es decir, la variable fls contiene los datos de precipitación y caudal para varias estaciones.Adicionalmente, cálculamos los datos atípicos de cada variable con el fin de hacer más fácil el proceso de identificación de correspondencia entre el ONI y el valor de la variable. Finalmente, se generan las gráficas de comparación y un archivo csv por cada estación que arca el valor atípico y el valor registrados por el ONI en esa msima fecha. Para que este proceso sea ágil cuando se maneja un gran volumen de información hidroclimatológica es importante crear lso espacios de almacenamiento antes de ejecutar el código. Para este ejericico se han creado en la dirección de trabajo actual las carpetas OutliersPlot/ y Outliers/.

for (i in fls){ # Para todos los archivos que contenga la variable fls
  var <- sub(".+/", "", i); var <- sub(".csv", "", var)
  data <- read.csv(i) # Datos
  f <- ifelse(grepl(":", data[1,1]), "%Y-%m-%d %T", "%Y-%m-%d")
  data <- as.xts(data[,-1], order.by = as.POSIXct(data[,1], format = f, tz = "GMT"))
  
  coef <- apply(as.matrix(data), 2, function(x) moments::kurtosis(x, na.rm = T)) %>%
    median
  coef <- 1.5 * 0.35*coef; coef <- ifelse(coef>5, 5, coef)
  
  dys <- yday(index(data))
  dys[dys == 366] <- 365
  
  for (j in 1:ncol(data)){ # Para cada estación - se lee por código
    cod <- colnames(data)[j] %>% substring(2)
    dj <- data[,j]
    outliers <- data; outliers[] <- NA
    for (k in 1:365){ # Para cada dma del aqo
      ind <- seq(k-15, k+15)
      ind[ind < 1] <- ind[ind < 1] + 365
      ind[ind > 365] <- ind[ind > 365] - 365
      
      dk <- dj[dys %in% ind] %>% as.numeric
      otl <- boxplot.stats(dk, coef)$out # Identificación de atípicos mediante boxplot
      
      if (any(dj[dys == k] %>% as.numeric %in% otl)){
        outliers[index(dj)[dys == k][dj[dys == k] %in% otl],j] <- T
      }
    }
    
     # Se genera como salida las gráficas de la variable hidroclimatológica y el ONI
    oni.plot <- oni[cut(index(dj), "month") %>% unique %>% as.POSIXct(tz = "GMT"),]
    oni.plot <- data.frame(time = index(oni.plot), value = oni.plot$ONI %>% as.numeric)
    otl <- dj[outliers[,j]==1]
    otl <- data.frame(time = index(otl), value = otl %>% as.numeric)
    
    dj <- data.frame(time = index(dj), value = dj %>% as.numeric)
    p <- ggplot(dj %>% as.data.frame, aes(time, value)) + geom_line(color = "gray60") +
      geom_point(data = otl, aes(time, value), color = "indianred4", size = 0.8) +
      labs(title = sprintf("%s: %s", var, cod), x = "Fechas", y = "")
      
    o <- ggplot(oni.plot, aes(time, value)) + geom_col() +
      geom_hline(yintercept = c(-0.5, 0.5)) + 
      labs(y = "ONI")
    
    tp <- ggarrange(p, o, nrow = 2, ncol = 1, heights = c(1, 0.8))
    ggsave(sprintf("OutliersPlot/%s_%s.png", var, cod))
    if (nrow(otl)>0){
      otl[,"ONI"] <- oni[cut(otl[,1], "month") %>% as.POSIXct(tz = "GMT"),"ONI"] %>% as.numeric
      write.csv(otl, sprintf("Outliers/%s_%s.csv", var, cod))
    }
  }
}