Identificación de atípicos y llenado de datos hidroclimatológicos con R
Camila García-Echeverri
21.02.2023
Este cuaderno de R muestra un ejercicio práctico del procesamiento de series de precipitación mensual a través del lenguaje R. La serie de datos a analizar fue descargada del portal DHIME y corresponde a los registros capturados por la estación El Delirio 21200130. El procedimiento a realizar es la identificación de atípicos con el método del boxplot, su eliminación y el llenado de los datos faltantes con el método racional deductivo.
Lo recomendado al iniciar a ejecutar/escribir un script es limpiar el entorno de trabajo. Adicionalmente, en R se recomienda iniciar con el llamado de las librerías a utilizarse y la definición de la carpeta o entorno de trabajo.
#Limpiar el entorno de trabajo
rm(list=ls())
#Cargar las librerias necesarias. En caso de no tenerlas instaladas instálelas usando:
# install.packages("xts") ; install.packages("lubridate") ; install.packages("rsoi")
library(xts)
library(lubridate)
library(rsoi)
setwd("Z:/My Drive/PERSONALES/TRABAJOS/DOCENCIA/CLASES R/CALIDAD DE DATOS")1 Caracterización inicial
Se realiza la lectura del archivo que contiene la serie de datos a analizar.
#Leer el archivo en formato .csv
datos <- read.csv("./21200130_DELIRIO_EL_P(mes).csv")
#El archivo tiene muchas columnas que para este caso no son de interés, así que selecciono las que incluyen la fecha y el valor.
datos <- datos[,17:18]
#convertir la serie en una variable temporal
datos <- as.xts(datos, order.by=as.Date(datos$Fecha))En este ejercicio vamos a analizar el periodo 1980-2010 (31 años), para lo cual debo recortar la serie a este periodo.
periodo_interes <- seq(as.Date("1980-01-01"), as.Date("2010-12-31"), "month")
datos <- datos[periodo_interes,]Los registros descargados del DHIME no incluyen filas vacias, es decir, que no incluyen fechas que por alguna razón no tienen registro. Por tanto, el número total de registros posibles no coincide con el número total de registros existentes. Teniendo en cuenta esta situación, se calcula el porcentaje de datos faltantes.
#Obtener el número de registros existentes
reg_existentes <- length(datos$Fecha)
#Obtener el número de registros posibles
fechas <- seq(as.Date(index(datos)[1]),as.Date(last(index(datos))), "month")
reg_posibles <- length(fechas)
#Calcular el procentaje de datos faltantes
print(paste0("El porcentaje de datos faltantes es: ", round((100 - reg_existentes/reg_posibles*100),1), " %"))## [1] "El porcentaje de datos faltantes es: 1.1 %"Vamos a crear una nueva variable que incluya las filas con datos vacios.
serie <- xts(order.by = fechas)
serie <- cbind(datos$Valor, serie)
rm(datos)2 Análisis de atípicos
Vamos a identificar los atípicos a través del diagrama de cajas y bigotes, los cuales serán construidos para cada uno de los meses del año. Para esto, lo primero que vamos a hacer es generar una columna en la variable serie que indique el mes al cual corresponde cada registro.
Adicionalmente, vamos a generar la variable de entrada a la función boxplot. Para crear 12 boxplot en un mismo gráfico lo más fácil es generar una variable tipo lista, donde cada elemento corresponda a cada mes.
serie$Mes <- month(index(serie))
boxp_input <- lapply(seq(1,12,1), function(x){
tmp <- which(serie$Mes==x)
tmp <- serie[tmp,1]
})Ahora, vamos a construir los boxplots mensuales.
#Generar la gráfica
boxplot(boxp_input, xlab="Mes", ylab="Precipitación (mm/mes)",
main="Estación El Delirio [21200130]", outcol=2)
Con la función boxplot() tambien se puede identificar el valor del atípico.
#obtener los valores atípicos por mes
valores_atipicos <- lapply(boxp_input, function(x){
tmp <- boxplot(x, plot=F)$out
})
valores_atipicos## [[1]]
## [1] 219.7
##
## [[2]]
## numeric(0)
##
## [[3]]
## numeric(0)
##
## [[4]]
## numeric(0)
##
## [[5]]
## numeric(0)
##
## [[6]]
## [1] 332.5 290.7
##
## [[7]]
## [1] 325.8 410.9
##
## [[8]]
## [1] 257.1
##
## [[9]]
## [1] 119.5
##
## [[10]]
## [1] 243.4
##
## [[11]]
## [1] 247.9
##
## [[12]]
## numeric(0)Contar el número de atípicos.
tmp <- unlist(lapply(valores_atipicos, function(x){
length(which(!is.na(x)))
}))
cat("Se encontraron", sum(tmp), "datos atípicos")## Se encontraron 9 datos atípicosEncontrar los meses que corresponden al valor atípico identificado.
#identificar cuales meses tuvieron atípicos
meses_con_atipicos <- which(tmp!=0)
rm(tmp)
#Encontrar la fecha exacta del registro identificado como atípico
registros_atipicos <- lapply(seq(1,12,1), function(x){
#aplico estas operaciones cuando en el mes se identificó más de un atípico
if(length(valores_atipicos[[x]])!=0){
tmp_save <- list()
for(i in 1:length(valores_atipicos[[x]])){
#aquí basicamente se esta buscando el valor atípico en el mes correspondiente y se está extrayendo la fecha de su registro
tmp <- which(valores_atipicos[[x]][i]==boxp_input[[x]])
tmp_save[[i]] <- boxp_input[[x]][tmp,]
}
}else{
#Es caso de que en el mes no hay atípicos no hago la busqueda y le asigno NA
tmp_save <- NA
}
if(class(tmp_save)=="list"){tmp_save <- do.call(rbind,tmp_save)}
return(tmp_save)
})
registros_atipicos## [[1]]
## Valor
## 1997-01-01 219.7
##
## [[2]]
## [1] NA
##
## [[3]]
## [1] NA
##
## [[4]]
## [1] NA
##
## [[5]]
## [1] NA
##
## [[6]]
## Valor
## 2002-06-01 332.5
## 2004-06-01 290.7
##
## [[7]]
## Valor
## 1994-07-01 325.8
## 1997-07-01 410.9
##
## [[8]]
## Valor
## 1991-08-01 257.1
##
## [[9]]
## Valor
## 2005-09-01 119.5
##
## [[10]]
## Valor
## 2006-10-01 243.4
##
## [[11]]
## Valor
## 2008-11-01 247.9
##
## [[12]]
## [1] NAEliminar los atípicos de las series mensuales
#Con las fechas de los atipicos lo que hago es ir a la serie original y ordenarle al algoritmo que a esa fecha le asigne un NA
series_sin_atipicos <- lapply(seq(1,12,1), function(x){
if(x %in% meses_con_atipicos){
tmp <- index(boxp_input[[x]]) %in% index(registros_atipicos[[x]])
tmp_serie <- boxp_input[[x]]
tmp_serie[tmp,] <- NA
}else{
tmp_serie <- boxp_input[[x]]
}
return(as.numeric(tmp_serie))
})3 LLenado de datos faltantes
Debido a que se identificaron atípicos, el porcentaje de datos a rellenar es mayor que el porcentaje de datos faltantes.
rellenar <- length(which(is.na(do.call(c,series_sin_atipicos))==T))
cat("El porcentaje de datos a rellenar es ", round(rellenar/reg_posibles*100,1), "%")## El porcentaje de datos a rellenar es 3.5 %Para realizar el llenado de faltantes utilizaremos el método racional deductivo. Para esto debemos organizar la serie por años.
#Lo organizo de tal forma que las filas corresponden a los meses y las columnas a los años
series_sin_atipicos <- do.call(rbind,series_sin_atipicos)
#esto lo hago para colocarle nombres a las columnas, cada columna corresponde a un año
colnames(series_sin_atipicos) <- unique(year(fechas))Vamos a retirar las columnas que contienen valores NA (faltantes) y clasificar los datos en dos variables. Una que contiene los años con datos faltantes y otra que contiene los años sin datos faltantes.
anos_con_faltantes <- apply(series_sin_atipicos, 2, function(x){
length(which(is.na(x)==T))
})
anos_sin_faltantes <- which(anos_con_faltantes==0)
anos_con_faltantes <- which(anos_con_faltantes!=0)
series_con_faltantes <- series_sin_atipicos[,anos_con_faltantes]
series_sin_faltantes <- series_sin_atipicos[,anos_sin_faltantes]Paso 1 del método racional deductivo: calcular la precipitación media anual de los años sin datos faltantes.
p_media_anual <- apply(series_sin_faltantes,2,mean)Paso 2 del método racional deductivo: calcular los porcentajes de cada registro de los años sin faltantes. El porcentaje se calcula respecto al total de cada año.
#calcular los porcentajes asociados a cada registro
porcentajes_meses <- lapply(1:dim(series_sin_faltantes)[2], function(x){
series_sin_faltantes[,x]/p_media_anual[x]*100
})
porcentajes_meses <- do.call(cbind, porcentajes_meses)
#verifico que todas las columnas sumen 1200
apply(porcentajes_meses,2,sum)## [1] 1200 1200 1200 1200 1200 1200 1200 1200 1200 1200 1200 1200 1200 1200 1200
## [16] 1200 1200 1200 1200 1200 1200Paso 3 del método racional deductivo: calculo el valor medio mensual de los porcentajes.
porc_media_mensual <- apply(porcentajes_meses,1,mean)Paso 4 del método racional deductivo: procedo a rellenar los faltantes aplicando la siguiente ecuación
\[\begin{equation} {P_{i}} = \frac{\sum (P)}{1200-\sum (S_{i})}*S_{i} \tag{3.1} \end{equation}\]Donde, i: cada uno de los meses desconocidos Pi: precipitación mensual desconocida en cada año incompleto. Si: porcentajes promedio de los meses desconocidos P:precipitaciones mensuales conocidas en años incompletos Si: porcentaje promedio del mes desconocido
#sumo los registros disponibles de cada año con faltantes
p_total_series_faltantes <- apply(series_con_faltantes, 2, sum, na.rm=T)
#en este ciclo recorro cada uno de los años con faltantes
#sumo los porcentajes medios mensuales de los meses que tienen faltantes
#busco los meses que no tienen registros (meses a rellenar) y con ellos itero para aplicar la ecuación
for(i in 1:dim(series_con_faltantes)[2]){
pos_na <- which(is.na(series_con_faltantes[,i])==T)
suma_porc_falt <- sum(porc_media_mensual[pos_na])
#Este ciclo se va a encargar de iterar cada uno de los datos por rellenar
for(f in 1:length(pos_na)){
pos_na_f <- pos_na[f]
#Aquí aplico la fórmula
llenado <- p_total_series_faltantes[i]/(1200-suma_porc_falt)*porc_media_mensual[pos_na_f]
series_con_faltantes[pos_na_f,i] <- llenado
}
}Lo último que queda es organizamos la serie rellenada y la guardamos en la ubicación deseada, en este caso voy a guardarlo en la misma carpeta de trabajo con un nombre diferente al archivo original.
serie_final <- matrix(nrow=12, ncol=length(unique(year(fechas))))
serie_final[,anos_sin_faltantes] <- series_sin_faltantes
serie_final[,anos_con_faltantes] <- series_con_faltantes
serie_final <- c(serie_final)
serie_final <- cbind(fechas, as.data.frame(serie_final))
colnames(serie_final) <- c("Mes", "P(mm)")
serie_final[,2] <- round(serie_final[,2],1)
write.csv(serie_final,"./21200130_DELIRIO_EL_P(mes)_rellenada.csv", row.names = F)