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) y Cristian Díaz (estudiante doctoral en la FEAR-PUJ).

El presente cuaderno muestra una introducción al procesamiento de series de datos hidrometeorológicos mediante el uso de R. Se muestra una caracterización estadística inicial de las series, un análisis exploratorio para la identificación de atípicos mediante el uso de diagramas de cajas y bigotes (boxplot) y una agrupación de los datos a escalas mensuales, anuales y multianuales, así como diferentes gráficas a lo largo de todo el proceso. Para el ejercicio se utilizan los datos diarios de precipitación en la estación pluviométrica de Nataga [21050090] operada por el IDEAM y localizada en el municipio del mismo nombre en el departamento del Huila. Esta información fue descargada directamente de la plataforma DHIME http://dhime.ideam.gov.co/atencionciudadano/.

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()) # Limpia el espacio de trabajo

# install.packages("xts")
# install.packages("lubridate")
# install.packages("moments")
# install.packages("ggplot2")
# install.packages("ggpubr")
# install.packages("tidyverse")

library(xts) # Series Temporales
library(lubridate) # Utilidades para manejo de fechas
library(moments) # Momentos estadísiticos
library(ggplot2) # Plot
library(ggpubr) # Plot
library(tidyverse) # Manejo de datos
library(DataExplorer) # Analisis de variables

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

Caracterización Inicial

Iniciamos cargando y leyendo el archivo descargado de DHIME. Para preservar las tildes y caracteres especiales se especifica el tipo de codificación como “UTF-8”.

De toda la información contenida en el archivo centraremos nuestro interés en las columnas de fecha y valor. Con estas dos columnas se crea una variable tipo “xts”, la cual se utiliza para el manejo de series temporales. El DHIME maneja las fechas en formato estandar AAAA-MM-DD, si se tuviera otro formato es necesario especificarlo; si se deseara utilizar las fechas junto con la hora de adquisición se debería utilizar variables tipo “POSIXct”.

En algunos casos, es necesario indicar también el separador (“,” ” ” “;”) que usa el archivo para consolidar la matriz de manera correcta. Para esto use el comando “sep= ”,””. Si el separador de datos es “;” por favor cambielo en la primer línea de este código.

# Lee el archivo del DHIME
datos <- read.csv("nataga.csv", encoding = "UTF-8") 
# Conocer los nombres de las variables de la base de datos 
names(datos) 
##  [1] "CodigoEstacion"   "NombreEstacion"   "Latitud"          "Longitud"        
##  [5] "Altitud"          "Categoria"        "Entidad"          "AreaOperativa"   
##  [9] "Departamento"     "Municipio"        "FechaInstalacion" "FechaSuspension" 
## [13] "IdParametro"      "Etiqueta"         "DescripcionSerie" "Frecuencia"      
## [17] "Fecha"            "Valor"            "Grado"            "Calificador"     
## [21] "NivelAprobacion"
# Organiza los valores como serie temporal
datos <- as.xts(datos$Valor, order.by = as.Date(datos$Fecha))

El IDEAM no incluye celdas cuando en ese día no se registró información, por lo que el número de días total de la serie y el número de registros varían.

cat("Número de días con registro: ", length(datos))
## Número de días con registro:  10944
cat("Número de días totales: ", max(index(datos))-min(index(datos))+1)
## Número de días totales:  10958

Para esto solucionar esto se utiliza un vector xts vacío pero cargado con la totalidad de los días que hay en la serie. La función cbind cuando se utiliza en variables xts hace la combinación de columnas emparejando las fechas.

# Fechas completas
fechas <- seq(min(index(datos)), max(index(datos)), by = "day") 
# Serie temporal vacia
empty.vector <- xts(order.by = fechas) 
# Serie con fechas completas
datos <- cbind(datos, empty.vector) 

cat("\fNúmero de días con registro: ", length(datos))
## Número de días con registro:  10958
cat("\nNúmero de días totales: ", max(index(datos))-min(index(datos))+1)
## 
## Número de días totales:  10958

Posteriormente, se hace la caracterización estadística inicial. Para el presente ejercicio se seleccionaron las estadísticas nombradas en la variable “nombres”. Al concluir la eliminación de valores atípicos y si se realizara llenado de datos faltantes estas estadísticas deben ser calculadas de nuevo.

# Estadísticas a calcular
nombres <- c("Media", "Desv. Estandar", "Asimetria", "Curtosis", 
"Máximo", "Mínimo", "Mediana",
"Coef. Variacíón", "Número Total de Datos",
"Número de Datos Faltantes", "Porcentaje de Faltantes")

# Calculo de las estadísticas
# Vector vacio para caracterización estadística
carc <- vector(length = length(nombres)) 
carc[1] <- mean(datos, na.rm = T, digits=2)
carc[2] <- sd(datos, na.rm = T)
carc[3] <- skewness(datos, na.rm = T)
carc[4] <- kurtosis(datos, na.rm = T)
carc[5] <- max(datos, na.rm = T)
carc[6] <- min(datos, na.rm = T)
carc[7] <- median(datos, na.rm = T)
carc[8] <- carc[2] / carc[1]
carc[9] <- sum(!is.na(datos))
carc[10] <- sum(is.na(datos))
carc[11] <- carc[10] / length(datos)

carc <- data.frame("Estadística" = nombres, Valor = carc)
carc
##                  Estadística        Valor
## 1                      Media 5.454130e+00
## 2             Desv. Estandar 1.097241e+01
## 3                  Asimetria 3.510595e+00
## 4                   Curtosis 2.024767e+01
## 5                     Máximo 1.320000e+02
## 6                     Mínimo 0.000000e+00
## 7                    Mediana 0.000000e+00
## 8            Coef. Variacíón 2.011762e+00
## 9      Número Total de Datos 1.094400e+04
## 10 Número de Datos Faltantes 1.400000e+01
## 11   Porcentaje de Faltantes 1.277605e-03

Recuerde que cuando trabaja con datos de Precipitación, las estadísticas corresponden a los valores temporales ingresados. En este ejemplo son datos con escala temporal diaria.

Ejercicio

¿Cómo ajustaría el código para que los datos estadísticos sean expresados en decimales de 2 cifras?

Gráficas

Dada la gran cantidad de datos es recomendable ayudarse continuamente de gráficas. En este punto se hará una gráfica combinada, donde se muestre la serie completa y el comportamiento mensual multianual mediante boxplots. Para las gráficas se va a emplear la librería ggplot, para esto primero se debe organizar la información en un data frame, el cual se introduce como primera variable y posteriormente se especifican los nombres de las columnas que serán graficadas. Mediante el signo “+” se indica el tipo de geometría que se utilizará (“geom_…”) y se le da formato al gráfico.

Para generar los boxplot mensuales se debe generar una columna que especifique el mes, para lo que se usa la función month y se especifica que se agrupara con esta columna.

# Plot de serie temporal y boxplot mensual multianual
datos.df <- data.frame(t = index(datos), datos)

p <- ggplot(datos.df, aes(t, datos)) + geom_line() +
  labs(title = "Precipitación en la Estación Nataga", x = "Años", y = "P (mm/día)") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"), text = element_text(size = 14))

datos.df[,1] <- month(index(datos)) # Mes de adquisición del dato

bp <- ggplot(datos.df, aes(t, datos, group = t)) +
  geom_boxplot(colors=t) + labs(x = "Meses", y = "P (mm/día") +
  scale_x_continuous(breaks = c(1:12)) +
  theme(text = element_text(size = 14))
tp <- ggarrange(p, bp, nrow = 2, ncol = 1, heights = c(1, 0.8))
tp

Análisis de Datos Atípicos

Los boxplot mostrados anteriormente son ilustrativos, en la presente sección se calcularán numéricamente de manera manual y automática. Estos diagramas se basan en los cuartiles, calculados con la función quantile. Y la distancia entre el primer y tercer cuartil, que se conoce como distancia intercuartílica (IR) que repesenta la diferencia entre el percentil 75 (Q3) y el percentil 25 (Q1) en un conjunto de datos. Mide la dispersión del 50% medio de los valores. Esta distancia se multiplica por un coeficiente (variable “c”) que toma valores entre 1.5 y 3, dada la gran variabilidad que los datos diarios de temperatura se recomienda utilizar 3 para esta variable. El resultado de esta multiplicación se resta al primer cuartil y se suma al tercero para obtener los limites por fuera de los cuales un dato se considera atípico.

c <- 3 # Coefiente para boxplots

dias.ppt <- datos
dias.ppt[dias.ppt == 0] <- NA # Solo se tienen en cuenta los dias con precipitación para el análisis

# Manual
cuart <- quantile(dias.ppt, na.rm = T) # Cuartiles
icr <- cuart[4] - cuart[2] # Distancia intercuartilica 
lims <- c(cuart[2]-c*icr, cuart[4]+c*icr) # Limites
names(lims) <- c("Inferior", "Superior")
lims
## Inferior Superior 
##      -33       51
atipicos.select.m <- datos < lims[1] | datos > lims[2] # TRUE cuando es atípico

R incluye en sus librerías por defecto una forma automática de hacer este análisis. Con la función boxplot.stats se obtienen los cuartiles en “stats”, el tamaño de la muestra en “n” y los valores atípicos en “out”. Con este último resultado solo basta verificar si los valores de la serie se encuentran dentro de la lista atípicos para marcarlos como tal, para esto se utiliza el comparador lógico “%in%”.

# Automático
stats <- boxplot.stats(as.numeric(dias.ppt), coef = c) # Boxplot automático
stats
## $stats
## [1]  0.2  3.0  7.0 15.0 51.0
## 
## $n
## [1] 5047
## 
## $conf
## [1] 6.733117 7.266883
## 
## $out
##   [1]  60.0  61.0 100.0  61.0  55.0  59.0  62.0  62.0  53.0  80.0  60.0 130.0
##  [13]  60.0  79.0  58.0  81.0  52.0  53.0  55.0  76.0  91.0  65.0 100.0  70.0
##  [25]  59.0  67.0  72.0  53.0  65.0  52.0  52.0  60.0  77.0  68.0  71.0  57.0
##  [37]  75.0  56.0  64.0  87.0  67.0  57.0  92.0  54.0  54.0  88.0  60.0  52.0
##  [49]  52.0  58.0  60.0  65.0 109.6  51.3  55.5  68.5  51.5  68.5  56.8  96.7
##  [61]  69.2  70.1  57.7  57.3  61.4  51.6  70.2  66.2  51.6  66.3  51.2  59.9
##  [73]  75.1  64.2  56.6  63.4 132.0  73.3  54.4  55.2  63.1  68.1  57.0  58.0
##  [85]  75.0  60.0  55.0 117.0  58.0  70.0  58.0  64.0  55.0  61.0  88.0  53.0
##  [97]  53.0  76.0  57.0  57.0  90.0  69.0  88.0  72.0  53.0  53.0  54.0  98.0
## [109]  55.0  68.0  58.0  68.0  56.0  53.0  74.0  65.0  78.0  59.0  63.0  59.0
## [121]  66.0  90.0
atipicos.select.a <- datos %in% stats$out # TRUE cuando es atípico

cat("\n¿Coinciden los resultados manuales con los automáticos? -", all(atipicos.select.a == atipicos.select.m, na.rm = T))
## 
## ¿Coinciden los resultados manuales con los automáticos? - TRUE
cat("\n", sum(atipicos.select.a, na.rm = T), "atípicos encontrados")
## 
##  122 atípicos encontrados

Nuevamente, es recomendable tener una ayuda visual, para esto se pueden graficar sobre la serie temporal los atípicos identificados.

atipicos <- data.frame(t = index(datos)[atipicos.select.a], datos[atipicos.select.a])

p1 <- p + geom_point(data = atipicos, aes(t, datos), color = "indianred4", size = 1)
p1

La eliminación de datos atípicos concluye con un proceso manual en el cual el profesional a cargo debe involucrar toda otra información disponible así como su propio criterio para decidir eliminar o no cada dato identificado como anómalo. Para esto se recomienda guardar un nuevo archivo de trabajo con la información necesaria para hacer este proceso, así como la gráfica generada previamente.

atipicos.csv <- as.data.frame(datos)
atipicos.csv[,2] <- atipicos.select.a

write.csv(atipicos.csv, "AnálisisExploratorio.csv")

Caracterización Final de la Serie

Luego de la eliminación de los datos atípicos (y de cualquier otro proceso que se haya realizado sobre la serie) se puede hacer la caracterización final de esta. Se pueden repetir los cálculos estadísticos, así como obtener información adicional como la de los valores mensuales, anuales o multianuales para aproximarse al comportamiento general de la zona de estudio.

En este punto, se recomienda limpiar nuevamente el espacio de trabajo, leer el archivo de datos sobre el que se trabajó y organizarlo nuevamente para trabajar como serie temporal.

rm(list  = ls())

datos <- read.csv("AnálisisExploratorio.csv")
datos <- as.xts(datos$datos, order.by = as.Date(datos$X))

Se recalculan los estadísticos y se verifica que la variación con los datos originales no sea significativa.

## Estadisticas
nombres <- c("Media", "Desv. Estandar", "Asimetria", "Curtosis", "Máximo", "Mínimo", "Mediana","Coef. Variacíón", "Número Total de Datos", "Número de Datos Faltantes", "Porcentaje de Faltantes")

# Calculo de las estadísticas
carc <- vector(length = length(nombres)) # Vector vacio para caracterización estadística
carc[1] <- mean(datos, na.rm = T)
carc[2] <- sd(datos, na.rm = T)
carc[3] <- skewness(datos, na.rm = T)
carc[4] <- kurtosis(datos, na.rm = T)
carc[5] <- max(datos, na.rm = T)
carc[6] <- min(datos, na.rm = T)
carc[7] <- median(datos, na.rm = T)
carc[8] <- carc[2] / carc[1]
carc[9] <- sum(!is.na(datos))
carc[10] <- sum(is.na(datos))
carc[11] <- carc[10] / length(datos)

carc <- data.frame("Estadística" = nombres, Valor = carc)
carc
##                  Estadística        Valor
## 1                      Media 5.454130e+00
## 2             Desv. Estandar 1.097241e+01
## 3                  Asimetria 3.510595e+00
## 4                   Curtosis 2.024767e+01
## 5                     Máximo 1.320000e+02
## 6                     Mínimo 0.000000e+00
## 7                    Mediana 0.000000e+00
## 8            Coef. Variacíón 2.011762e+00
## 9      Número Total de Datos 1.094400e+04
## 10 Número de Datos Faltantes 1.400000e+01
## 11   Porcentaje de Faltantes 1.277605e-03

Para agrupar mensual o anualmente se utiliza la función aggregate, en donde se introduce un indice bajo el cual se desea que se agrupen los datos, se especifica la manera de agruparlos y cualquier elemento extra que se considere. Para obtener este indice a partir de las fechas diarias se “corta” la fecha con la función cut y se especifica a que resolución temporal se desear cortar. Este proceso cambia el tipo de variable, así que se debe volver a definirla como una serie temporal (xts).

Cuando la serie tiene numerosos datos faltantes se debe especificar una función que descarte un mes o año cuando a este le falta cierto porcentaje de datos.

## Agrupación
mensual <- aggregate(datos, by = cut(index(datos), "month"), sum, na.rm = T)
mensual <- as.xts(mensual, order.by =as.Date(index(mensual)))
head(mensual)
##            [,1]
## 1992-01-01  130
## 1992-02-01  129
## 1992-03-01   56
## 1992-04-01  108
## 1992-05-01  252
## 1992-06-01   34
anual <- aggregate(datos, by = cut(index(datos), "year"), sum, na.rm = T)
anual <- as.xts(anual, order.by = as.Date(index(anual)))
head(anual)
##            [,1]
## 1992-01-01 1517
## 1993-01-01 2236
## 1994-01-01 2636
## 1995-01-01 1787
## 1996-01-01 2508
## 1997-01-01 1971

En precipitación se debe tener en cuenta que para pasar los datos de escala temporal diaria a mensual o anual, se suman pero para pasar a valores multianuales se promedia. Este tratamiento difiere de variable a variable (por ejemplo para Temperatura).

mensual.ma <- aggregate(mensual, by = month(index(mensual)), mean, na.rm = T)
mensual.ma
##             
## 1  174.08333
## 2  186.98000
## 3  234.84000
## 4  217.41667
## 5  178.68333
## 6  107.44333
## 7   82.52000
## 8   56.96000
## 9   79.04333
## 10 185.04000
## 11 266.19000
## 12 220.46667
anual.ma <- mean(anual, na.rm = T)
anual.ma
## [1] 1989.667

Finalmente, se puede graficar el comportamiento anual promedio de la zona. Incluso se pueden mostrar varias variables en una sola gráfica para complementar la información.

temp <- c(30.0, 28.6, 27.9, 28.1, 29.7, 30.2, 31.0, 31.5, 30.7, 28.8, 27.9, 28.6)
mensual.ma <- data.frame(m = 1:12, m.abb = month.abb, ppt = mensual.ma, temp = temp)
p0 <- ggplot(mensual.ma, aes(m, ppt)) + geom_col(fill = "darkslategray3") +
  scale_x_continuous(breaks = 1:12, labels = mensual.ma$m.abb) +
  labs(title = "Precipitación Mensual Multianual en Nátaga", subtitle = sprintf("Precipitación total = %.f mm", anual.ma), x = "Meses", y = "P (mm)") +
  theme(plot.title = element_text(hjust = 0.5))
p0

coef <- 0.2
p1 <- ggplot(mensual.ma, aes(m)) +
  geom_col(aes(y = ppt), fill = "darkslategray3") +
  geom_line(aes(y = temp/coef), color = "firebrick3") +
  scale_x_continuous(breaks = 1:12, labels = mensual.ma$m.abb) +
  scale_y_continuous(name = "Precipitación (mm)", sec.axis = sec_axis(~.*coef, name = "Temperatura (ºC)")) +
  labs(title = "Caracterización Mensual Multianual en Nátaga", subtitle = sprintf("Precipitación total = %.f mm \nTemperatura Media = %.1f ºC", anual.ma, mean(temp)), x = "Meses") +
  theme(plot.title = element_text(hjust = 0.5))
p1