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/. Por favor realice este procedimiento de manera previa y guarde el archivo en su carpeta de trabajo.

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

A partir de aquí usted puede facilitar la transcripción del código copiando las líneas que se encuentran en los recuadro grises. Los recuadros blancos indican los resultados que usted debe visualizar al ejecutar de manera correcta las instrucciones.

rm(list = ls()) # Limpia el espacio de trabajo

# Si es la primera vez que realiza el procedimiento, por favor elimine el símbolo # de todas las líneas que dicen install.packages

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

library(xts) # Series Temporales
library(lubridate) # Utilidades para manejo de fechas
library(moments) # Momentos estadísiticos
library(RColorBrewer) # colores para gráficos
library(ggplot2) # Plot
library(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.

setwd("C:\\Users\\ASUS\\Downloads") # Carpeta de trabajo

Caracterización Inicial

Nota Importante: este tipo de herramientas (Rsutdio, Python) están pensados para escritura en inglés. En este idioma no hay un uso frecuente de tildes y caracteres especiales, así que para preservar las tildes y caracteres especiales se debe especificar el tipo de codificación como “UTF-8”.

Ahora si, entremos en materia. Para el desarrollo de sus actividades académicas seguramente será muy útil relacionar información ambiental e hidrológica con fechas (series de tiempo). Dado esto, el interés de este ejercicio de clase se centrará 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.

La mayoría de fuentes oficiales de información ambiental 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("nataga.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

Cuando trabajamos información hidroclimatológica (precipitación, caudal, temperatura, velocidad del viento, humedad realtiva, brillo solar, entre otras), se debe tener especial cuidado en la fuente, dado que el IDEAM, por ejemplo, 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.

Recuerde que en hidrología el valor cero (0) es diferente del dato vacio (NAN). El primero indica, por ejemplo en precipitación que ese día no llovió (0 mm), el segundo indica que ese día no se realizó la medición de la variable.

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

# 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, donde se muestre la serie completa de datos de Precipitación en el tiempo. Para las siguientes 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.

# Plot de serie temporal 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))
p

A estos gráficos se les puede cambiar el color y el tipo de línea. Probemos con líneas azules!.

# Plot de serie temporal mensual multianual
datos.df <- data.frame(t = index(datos), datos)
p <- ggplot(datos.df, aes(t, datos)) + geom_line(colour = "blue") +
  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))
p

Ejercicio:

  1. Identificamos cuál fue el comando usado para realizar el cambio de color?

  2. Genere la misma gráfica en color: verde, rosa, naranja

  3. Exporte la gráfica usando la función export que esta disponible en el menú Plots del lado izquierdo de la pantalla.

y si ahora probamos con el grosor y tipo de la línea. Los comandos que nos ayudan con estas funciones son: lwd (para definir el ancho) y linetype (para definir el tipo de línea que queremos crear).

# Plot de serie temporal mensual multianual
datos.df <- data.frame(t = index(datos), datos)
p <- ggplot(datos.df, aes(t, datos)) + geom_line(colour = "blue",lwd = 2, linetype = 2) +
  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))
p

Ejercicio:

1.Revise la página https://r-charts.com/es/r-base/tipos-lineas/ y realice diferentes pruebas de ancho y tipo de línea.

BoxPlot

Bien ahora hablemos un poco más del análisis estadístico mediante el uso de boxplot. Recordemos…los gráficos de caja (boxplot) fueron originalmente desarrollados por M.E. Spear, permiten conocer cómo se distribuyen los datos dentro de una variable. A diferencia de los histogramas que requieren un tamaño de muestra de al menos 30 casos para ser útiles, los gráficos de caja pueden ser construidos con tan solo 5 casos y aportan más detalles acerca de las colas de la distribución. Representan la información que se observa en la figura:

Mediana. Valor que deja a la mitad de los casos por encima y a la otra mitad por debajo.

Primer Cuartil (Q1). El 25% de los casos se encuentran por debajo de este valor.

Tercer Cuartil (Q3). El 75% de los casos se encuentran por encima de este valor.

Rango Intercuartílico (RIC). Es la diferencia entre el tercer y el primer cuartil.

Límites Superior o Inferior (Ls o Li). Ls contiene los casos por encima de Q3 más 1,5 veces el rango intercuartílico o Li por debajo de Q1 – 1,5xRIC (Estilo de Tukey). Cuando los valores no son posibles en lugar de emplear la aproximación anterior se escogen los valores máximo o mínimo de la muestra (Estilo de Spears).

Los valores atípicos son aquellos que están más a allá de los límites inferior y superior. Cuando los valores atípicos están más allá de 3 veces el RIC en lugar del 1.5 son denominados valores extremos.

Para generar los gráficos de boxplot en una escala mensual 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 boxplot mensual multianual
datos.df <- data.frame(t = index(datos), datos)
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))
bp

Ejercicio

Juguemos con los gráficos de boxplot.

  1. Resaltar la media (mean) dentro del gráfico usando la función stat_summary(Color). Para el ejemplo la podrá observar en color rojo, haga cambios para 4 colores diferentes.
# Plot boxplot mensual multianual
datos.df <- data.frame(t = index(datos), datos)
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))+
  stat_summary(fun=mean, geom="point", color="red")
bp

  1. cambiar el color del borde y resaltar los datos atípicos a través de la función geom_boxplot(color=“red”, outlier.colour=“red”).
# Plot boxplot mensual multianual
datos.df <- data.frame(t = index(datos), datos)
datos.df[,1] <- month(index(datos)) # Mes de adquisición del dato
bp <- ggplot(datos.df, aes(t, datos, group = t)) +
  geom_boxplot(color="red", outlier.colour="red") + labs(x = "Meses", y = "mm/h") +
  scale_x_continuous(breaks = c(1:12)) +
  theme(text = element_text(size = 14))
bp

  1. y si queremos omitir los datos atípicos… qué se les ocurre?
# Plot boxplot mensual multianual
datos.df <- data.frame(t = index(datos), datos)
datos.df[,1] <- month(index(datos)) # Mes de adquisición del dato
bp <- ggplot(datos.df, aes(t, datos, group = t)) +
  geom_boxplot(color=4, outlier.colour=NA) + labs(x = "Meses", y = "mm/h") +
  scale_x_continuous(breaks = c(1:12)) +
  theme(text = element_text(size = 14))
bp

  1. Por qué no tenemos datos atípicos de los valores bajos? Si usaramos otra variable u otro sistema de unidades, esto podría cambiar. Explique su respuesta.

Pero y si quisieramos tener boxplot con otras formas… no cuadrados.

# Plot boxplot mensual multianual
datos.df <- data.frame(t = index(datos), datos)
datos.df[,1] <- month(index(datos)) # Mes de adquisición del dato
bp <- ggplot(datos.df, aes(t, datos, group = t)) +
  geom_boxplot(notch=TRUE, notchwidth = 0.5,
               color="blue", alpha = 0.7,
               outlier.colour="red", outlier.size=0.7) + labs(x = "Meses", y = "mm/h") +
  scale_x_continuous(breaks = c(1:12)) +
  theme(text = element_text(size = 14))
bp

  1. Qué otras formas de boxplot encuentran?
# Plot boxplot mensual multianual
datos.df <- data.frame(t = index(datos), datos)
datos.df[,1] <- month(index(datos)) # Mes de adquisición del dato
bp <- ggplot(datos.df, aes(t, datos, group = t)) +
  geom_violin(trim=FALSE, fill="lightblue")+
  geom_boxplot(width=0.1) + labs(x = "Meses", y = "mm/h") +
  scale_x_continuous(breaks = c(1:12)) +
  theme(text = element_text(size = 14))
bp

Series de Tiempo - Agrupadas

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.

Nota: si realiza el proceso para temperatura, esta no se deberá sumar en ninguna agregación.

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. En la última gráfica se puede ver la combinación de datos de precipitación y temperatura. Los datos de temperatura fueron creados estadísticamente y no corresponden a los datos reales generados en la zona.

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