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 utilizaron 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 dpto. del Huila, esta información fue descargada directamente de la plataforma DHIME.

Antes de iniciar con el tema en si, es recomendable empezar el código 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")

library(xts) # Series Temporales
library(lubridate) # Utilidades para manejo de fechas
library(moments) # Momentos estadísiticos
library(ggplot2) # Plot
library(ggpubr) # Plot

setwd("~/Diego/02.Profesionales/UPJ/DATOS/") # Carpeta de trabajo

1 Caracterización Inicial

Inicialmente se lee el archivo, para preservar las tildes y caracteres especiales se especifica el tipo de codificación como “UTF-8”. De la totalidad de la información del archivo, el interés se centra en las fechas y los valores. Con estas dos columnas se crea una variable tipo “xts”, las cuales se utilizan para el manejo de series temporales. El DHIME maneja las fechas en el 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.

datos <- read.csv("excel.csv.csv", encoding = "UTF-8") # Lee el archivo del DHIME
datos <- as.xts(datos$Valor, order.by = as.Date(datos$Fecha)) # Organiza los valores como serie temporal

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

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

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 Estacion Nataga", x = "Fechas", y = "mm/h") +
  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() + labs(x = "Meses", y = "mm/h") +
  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

2 Análisis de 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. 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

Las eliminación de 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") 
ggsave("Análisi exploratorio.png") # Guarda la gráfica

3 Caracterización Final de la Serie

Luego de la eliminación del 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 de diarios a mensuales o anuales se suman pero para pasar a valores multianuales se promedia. Este tratamiento difiere de variable a variable.

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