1 Introducción

Antes de usar una serie de precipitación en hidrología es indispensable verificar su calidad. Los registros de estaciones pluviométricas pueden tener:

Inhomogeneidades: cambios en la serie que no se deben al clima, sino a causas artificiales (cambio de ubicación del pluviómetro, cambio de instrumento, crecimiento de vegetación o construcciones alrededor, cambio de observador).
Valores atípicos: datos muy alejados del comportamiento típico, que pueden ser errores (digitación, unidades) o eventos extremos reales.
Datos faltantes (fallas del observador o del equipo).

En este documento se desarrolla, paso a paso, un flujo de trabajo en R para:

Cargar y organizar series diarias de tres estaciones (Carrizal, Balsayaco y Papayal).
Agregar las series a escala mensual y anual controlando los datos faltantes.
Evaluar la homogeneidad con la curva de doble masa y la prueba de Pettitt.
Detectar valores atípicos con diagramas de caja (boxplots) y decidir cuáles eliminar.

Datos necesarios. Tres archivos CSV en la carpeta de trabajo: CARRIZAL.csv, BALSAYACO.csv y PAPAYAL.csv. Cada uno debe tener dos columnas: Fecha (formato AAAA-MM-DD) y Valor (precipitación diaria en mm).

2 Paquetes

library(ggplot2)    # gráficos
library(ggpubr)     # ggarrange(): organizar varios gráficos en una figura
library(xts)        # series de tiempo con índice de fechas
library(lubridate)  # manejo de fechas: year(), month(), ...
library(trend)      # pruebas de tendencia y cambio: pettitt.test()

Paquete	Para qué se usa aquí
`xts`	Estructura de datos para series de tiempo: una matriz cuyo índice son fechas. Permite unir series por fecha automáticamente.
`lubridate`	Extraer año (`year()`) o mes (`month()`) de una fecha.
`trend`	Contiene `pettitt.test()`.
`ggplot2` / `ggpubr`	Gráficos y composición de figuras.

3 Carga y organización de los datos

3.1 Lectura de los archivos

La idea es leer cada archivo y unirlo en una sola tabla donde cada columna sea una estación y cada fila una fecha.

files <- c("CARRIZAL.csv", "BALSAYACO.csv", "PAPAYAL.csv")

series_diarias <- NULL
for (file in files) {
  diaria <- read.csv(file)

  fechas <- as.Date(diaria$Fecha)
  diaria <- as.xts(diaria$Valor, order.by = fechas)
  series_diarias <- cbind(series_diarias, diaria)
}

¿Qué hace este ciclo?

read.csv(file) lee el archivo como un data.frame.
as.Date(diaria$Fecha) convierte el texto de la columna Fecha en objetos de tipo fecha.
as.xts(valores, order.by = fechas) crea la serie de tiempo.
cbind() une la nueva serie con las anteriores. En xts, cbind alinea por fecha: si una estación no tiene dato en un día donde otra sí, se rellena con NA. Así las estaciones no necesitan tener períodos idénticos.

3.2 Completar el calendario

Si una estación no tiene filas para días donde no hubo medición, esas fechas simplemente no aparecen en el archivo. Para trabajar con una serie continua (un dato por día, aunque sea NA) se construye la secuencia completa de fechas y se une a la serie.

fecha_inicial <- min(index(series_diarias))
fecha_final   <- max(index(series_diarias))
fechas        <- seq(fecha_inicial, fecha_final, by = "day")

series_diarias <- cbind(
  series_diarias,
  xts(order.by = fechas)   # serie vacía solo con el calendario completo
)

colnames(series_diarias) <- sub("\\.csv", "", files)   # nombres de columnas

xts(order.by = fechas) crea una serie sin columnas pero con todas las fechas; al hacer cbind con ella se agregan las fechas faltantes con NA. Finalmente se renombran las columnas quitando la extensión .csv con sub().

head(series_diarias)

##            CARRIZAL BALSAYACO PAPAYAL
## 2000-01-01        2       0.0       4
## 2000-01-02        0       1.8       0
## 2000-01-03        8       4.2       0
## 2000-01-04        0       0.0       3
## 2000-01-05        0       0.3       5
## 2000-01-06        3       0.3      19

summary(series_diarias)

##      Index               CARRIZAL         BALSAYACO         PAPAYAL       
##  Min.   :2000-01-01   Min.   :  0.000   Min.   : 0.000   Min.   :  0.000  
##  1st Qu.:2006-09-07   1st Qu.:  0.000   1st Qu.: 0.200   1st Qu.:  0.000  
##  Median :2013-05-15   Median :  3.000   Median : 2.000   Median :  0.000  
##  Mean   :2013-05-15   Mean   :  6.753   Mean   : 5.301   Mean   :  5.315  
##  3rd Qu.:2020-01-21   3rd Qu.:  9.000   3rd Qu.: 7.200   3rd Qu.:  8.500  
##  Max.   :2026-09-28   Max.   :110.000   Max.   :85.600   Max.   :150.000  
##                       NA's   :372       NA's   :293      NA's   :443

4 Agregación mensual y anual

Muchos análisis se hacen con totales mensuales o anuales. Pero sumar ignorando los faltantes subestima la precipitación real: un mes con 25 días sin dato no puede compararse con un mes completo. Por eso se define una regla de calidad.

4.1 Regla de agregación

agregar <- function(x) {
  if (sum(is.na(x)) > 0.2 * length(x)) {
    return(NA)
  } else {
    return(sum(x, na.rm = TRUE))
  }
}

La función agregar() recibe los datos de un período (un mes o un año) y:

Si más del 20 % de los datos son NA, devuelve NA (el período se considera no confiable).
En caso contrario, suma los datos disponibles ignorando los NA.

El umbral del 20 % es un criterio de trabajo; otros autores o normas usan valores distintos (por ejemplo, 10 % o 25 %). Es una decisión del analista que debe documentarse.

4.2 Serie mensual

series_mensuales <- aggregate(series_diarias,
                              by = cut(index(series_diarias), "month"),
                              agregar)
series_mensuales <- as.xts(series_mensuales, as.Date(index(series_mensuales)))
head(series_mensuales)

##            CARRIZAL BALSAYACO PAPAYAL
## 2000-01-01      236     112.3     195
## 2000-02-01      199     159.7     164
## 2000-03-01      277     135.9      67
## 2000-04-01      412     196.5     129
## 2000-05-01      793     415.7     184
## 2000-06-01      616     286.2     141

cut(index(...), "month") clasifica cada fecha en su mes (devuelve el primer día del mes como etiqueta).
aggregate() aplica agregar() a cada grupo.
La última línea reconvierte el índice a fechas para que el objeto siga siendo una serie de tiempo.

4.3 Serie anual

series_anuales <- aggregate(series_diarias,
                            by = year(index(series_diarias)),
                            agregar)

# El índice resultante son años (2001, 2002, ...). Se convierte a fecha (1 de enero)
series_anuales <- as.xts(series_anuales,
                         as.Date(paste0(index(series_anuales), "-01-01")))
head(series_anuales)

##            CARRIZAL BALSAYACO PAPAYAL
## 2000-01-01   3908.0    2229.0    1927
## 2001-01-01   3908.2    2043.4      NA
## 2002-01-01   3371.0    1825.3    1654
## 2003-01-01   2992.0    1715.6    2001
## 2004-01-01   2704.0    1983.5    1363
## 2005-01-01   3194.9    2027.8      NA

Nota. El índice de la agregación anual es un número entero (el año). Si se pasara directamente a as.Date(), R lo interpretaría como “días desde 1970” y las fechas serían incorrectas; por eso se construye la fecha con paste0(año, "-01-01").

5 Análisis de homogeneidad

Una serie es homogénea cuando sus variaciones se deben únicamente a la variabilidad climática. Se presentan dos métodos complementarios: uno gráfico (curva de doble masa) y uno estadístico (prueba de Pettitt).

5.1 Curva de doble masa

5.1.1 Fundamento

La curva de doble masa grafica la precipitación acumulada de la estación en análisis contra la precipitación acumulada de una estación de referencia (o del promedio de varias) que se supone homogénea y correlacionada con la primera.

Si ambas estaciones registran el mismo régimen de lluvias, sus acumulados crecen de forma proporcional y los puntos siguen una línea recta. Un cambio de pendiente sostenido indica que en algún momento la relación entre las estaciones cambió, es decir, que la serie en análisis (o la de referencia) presenta una inhomogeneidad.

5.1.2 Preparación de los datos

series_homogeneidad <- series_diarias[, c("CARRIZAL", "BALSAYACO")] |>
  na.omit() |>
  as.data.frame()

series_homogeneidad <- cumsum(series_homogeneidad)
head(series_homogeneidad)

Se seleccionan las dos estaciones.
na.omit() elimina todos los días en que alguna de las dos no tiene dato. Esto es importante: los acumulados solo son comparables si se construyen con los mismos días en ambas estaciones.
as.data.frame() convierte a tabla para usar ggplot2.
cumsum() calcula el acumulado de cada columna.

5.1.3 Gráfico

p <- ggplot(series_homogeneidad, aes(BALSAYACO, CARRIZAL)) +
  geom_line(linewidth = 0.7) +
  labs(
    x = "Prec. Acumulada Referencia",
    y = "Prec. Acumulada Análisis",
    title = "Curva Doble Masa",
    subtitle = "Análisis: Carrizal\nReferencia: Balsayaco"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

p

¿Cómo interpretarla?

Una curva casi recta: no hay evidencia de inhomogeneidad.
Un quiebre claro de pendiente: revisar la historia de la estación (metadatos) para la fecha correspondiente. Si se confirma un cambio en las condiciones de medición, los datos de un tramo se pueden corregir multiplicándolos por la razón entre pendientes.

Limitación. El método es visual y subjetivo; pequeños cambios de pendiente son difíciles de detectar a ojo. Por eso se complementa con una prueba estadística.

5.2 Prueba de Pettitt

5.2.1 Fundamento

La prueba de Pettitt es una prueba no paramétrica (no supone que los datos sigan una distribución normal) para detectar un cambio abrupto en la media de una serie en un punto desconocido.

H₀ (hipótesis nula): no hay punto de cambio; la serie es homogénea.
H₁ (hipótesis alternativa): existe un punto de cambio.

pettitt.test() devuelve:

statistic (K): el valor del estadístico.
p.value: la probabilidad de observar un K igual o mayor si H₀ fuera cierta.
estimate: la posición (índice) dentro de la serie donde ocurre el cambio más probable.

Criterio de decisión: si p.value < 0.05 (nivel de significancia del 5 %), se rechaza H₀ y se concluye que hay un cambio significativo.

Importante. pettitt.test() no admite valores faltantes, por eso se aplica na.omit() antes. Como consecuencia, estimate es la posición dentro de la serie sin NA, y para recuperar la fecha hay que usar el índice de esa misma serie sin NA (index(x)[estimate]), no el de la serie completa.

5.2.2 Serie diaria

x_carrizal <- na.omit(series_diarias$CARRIZAL)
pettit_diario_carrizal <- pettitt.test(as.numeric(x_carrizal))
pettit_diario_carrizal

## 
##  Pettitt's test for single change-point detection
## 
## data:  as.numeric(x_carrizal)
## U* = 1513938, p-value = 1.265e-07
## alternative hypothesis: two.sided
## sample estimates:
## probable change point at time K 
##                            1717

x_balsayaco <- na.omit(series_diarias$BALSAYACO)
pettit_diario_balsayaco <- pettitt.test(as.numeric(x_balsayaco))
pettit_diario_balsayaco

## 
##  Pettitt's test for single change-point detection
## 
## data:  as.numeric(x_balsayaco)
## U* = 1331892, p-value = 7.369e-06
## alternative hypothesis: two.sided
## sample estimates:
## probable change point at time K 
##                            5504

Ahora se grafican las series y, solo si la prueba es significativa, se agrega una línea vertical en la fecha del cambio:

p_diario_carrizal <- data.frame(fecha = index(series_diarias),
                                valor = as.numeric(series_diarias$CARRIZAL)) |>
  ggplot(aes(fecha, valor)) +
  geom_line(color = "#77F") +
  labs(x = "Fecha", y = "mm", title = "Carrizal (diaria)") +
  theme_minimal()

if (pettit_diario_carrizal$p.value < 0.05) {
  p_diario_carrizal <- p_diario_carrizal +
    geom_vline(xintercept = index(x_carrizal)[pettit_diario_carrizal$estimate],
               linewidth = 0.9)
}

p_diario_balsayaco <- data.frame(fecha = index(series_diarias),
                                 valor = as.numeric(series_diarias$BALSAYACO)) |>
  ggplot(aes(fecha, valor)) +
  geom_line(color = "#77F") +
  labs(x = "Fecha", y = "mm", title = "Balsayaco (diaria)") +
  theme_minimal()

if (pettit_diario_balsayaco$p.value < 0.05) {
  p_diario_balsayaco <- p_diario_balsayaco +
    geom_vline(xintercept = index(x_balsayaco)[pettit_diario_balsayaco$estimate],
               linewidth = 0.9)
}

ggarrange(p_diario_carrizal, p_diario_balsayaco, nrow = 2)

5.2.3 Serie mensual

x_carrizal <- na.omit(series_mensuales$CARRIZAL)
pettit_mensual_carrizal <- pettitt.test(as.numeric(x_carrizal))
pettit_mensual_carrizal

## 
##  Pettitt's test for single change-point detection
## 
## data:  as.numeric(x_carrizal)
## U* = 6654, p-value = 0.0002318
## alternative hypothesis: two.sided
## sample estimates:
## probable change point at time K 
##                              76

x_balsayaco <- na.omit(series_mensuales$BALSAYACO)
pettit_mensual_balsayaco <- pettitt.test(as.numeric(x_balsayaco))
pettit_mensual_balsayaco

## 
##  Pettitt's test for single change-point detection
## 
## data:  as.numeric(x_balsayaco)
## U* = 2574, p-value = 0.5425
## alternative hypothesis: two.sided
## sample estimates:
## probable change point at time K 
##                             205

p_mensual_carrizal <- data.frame(fecha = index(series_mensuales),
                                 valor = as.numeric(series_mensuales$CARRIZAL)) |>
  ggplot(aes(fecha, valor)) +
  geom_line(color = "#77F") +
  labs(x = "Fecha", y = "mm/mes", title = "Carrizal (mensual)") +
  theme_minimal()

if (pettit_mensual_carrizal$p.value < 0.05) {
  p_mensual_carrizal <- p_mensual_carrizal +
    geom_vline(xintercept = index(x_carrizal)[pettit_mensual_carrizal$estimate],
               linewidth = 0.9)
}

p_mensual_balsayaco <- data.frame(fecha = index(series_mensuales),
                                  valor = as.numeric(series_mensuales$BALSAYACO)) |>
  ggplot(aes(fecha, valor)) +
  geom_line(color = "#77F") +
  labs(x = "Fecha", y = "mm/mes", title = "Balsayaco (mensual)") +
  theme_minimal()

if (pettit_mensual_balsayaco$p.value < 0.05) {
  p_mensual_balsayaco <- p_mensual_balsayaco +
    geom_vline(xintercept = index(x_balsayaco)[pettit_mensual_balsayaco$estimate],
               linewidth = 0.9)
}

ggarrange(p_mensual_carrizal, p_mensual_balsayaco, nrow = 2)

6 Análisis de valores atípicos

Un valor atípico es una observación que se aleja de forma marcada del resto. En hidrología hay que tener cuidado: un atípico no es necesariamente un error. Una tormenta extrema es un dato real y valioso, especialmente en diseño hidrológico. La estrategia es detectar automáticamente y decidir manualmente.

6.1 Criterio del diagrama de caja

El diagrama de caja marca como atípico todo valor que quede fuera del intervalo:

\[[\,Q_1 - 1.5\,\mathrm{RIC},\;\; Q_3 + 1.5\,\mathrm{RIC}\,]\]

donde $Q_1$ y $Q_3$ son el primer y tercer cuartil y $\mathrm{RIC} = Q_3 - Q_1$ es el rango intercuartílico. En R, boxplot() devuelve (además del gráfico) una lista cuyo elemento out contiene los valores atípicos.

Trabajaremos con la estación Papayal.

p_papayal <- data.frame(fecha = index(series_diarias),
                        valor = as.numeric(series_diarias$PAPAYAL)) |>
  ggplot(aes(fecha, valor)) +
  geom_line(color = "#77F") +
  theme_minimal()

p_papayal

6.2 Primer intento: con todos los datos (incluidos los ceros)

atipicos <- as.numeric(series_diarias$PAPAYAL) |>
  boxplot(main = "Papayal: todos los días")

length(atipicos$out)   # cuántos valores marca como atípicos

## [1] 564

Ubicamos en el tiempo los atípicos identificados. Se usa %in% para buscar las fechas cuyo valor está en la lista de atípicos:

serie_atipicos <- series_diarias$PAPAYAL[series_diarias$PAPAYAL %in% atipicos$out]
serie_atipicos <- data.frame(fecha = index(serie_atipicos),
                             valor = as.numeric(serie_atipicos))

p_papayal +
  geom_point(data = serie_atipicos, aes(fecha, valor)) +
  labs(x = "Fecha", y = "mm", title = "Atípicos considerando todos los días")

Problema: en una serie diaria de lluvia la mayoría de los días son secos (valor 0). Esto hace que $Q_1$, la mediana y con frecuencia $Q_3$ sean cercanos a cero, el RIC sea muy pequeño y, por tanto, el límite superior quede muy bajo. El resultado es que casi cualquier día lluvioso moderado aparece como “atípico”.

6.3 Segundo intento: solo días con lluvia

Lo razonable es preguntarse si un día es atípico entre los días en que llovió. Para eso se excluyen los ceros (y los NA).

papayal <- na.omit(series_diarias$PAPAYAL)
papayal_lluvia <- papayal[papayal > 0]

atipicos <- as.numeric(papayal_lluvia) |>
  boxplot(main = "Papayal: solo días con lluvia (> 0 mm)")

length(atipicos$out)

## [1] 270

serie_atipicos <- papayal_lluvia[papayal_lluvia %in% atipicos$out]
serie_atipicos <- data.frame(fecha = index(serie_atipicos),
                             valor = as.numeric(serie_atipicos))

p_papayal +
  geom_point(data = serie_atipicos, aes(fecha, valor)) +
  labs(x = "Fecha", y = "mm", title = "Atípicos considerando solo días con lluvia")

La cantidad de puntos marcados disminuye notablemente y ahora sí representa eventos realmente elevados respecto a los días lluviosos.

6.4 Revisión manual

Se exporta la lista de atípicos a un archivo CSV para revisarla fuera de R:

write.csv(serie_atipicos, "Atipicos_papayal.csv", row.names = FALSE)

atipicos_revisados <- read.csv("Atipicos_papayal_revisados.csv")
atipicos_revisados$fecha <- as.Date(atipicos_revisados$fecha)
atipicos_revisados$X <- atipicos_revisados$X == "X"   # TRUE = eliminar

atipicos_revisados

atipicos_revisados$X == "X" convierte la columna en lógica: TRUE para los datos que se van a eliminar y FALSE para los que se aceptan (las celdas vacías no son iguales a "X").

6.5 Eliminación de los atípicos rechazados

series_diarias_procesadas <- series_diarias$PAPAYAL
series_diarias_procesadas[atipicos_revisados$fecha[atipicos_revisados$X]] <- NA
series_diarias_procesadas <- data.frame(fecha = index(series_diarias_procesadas),
                                        valor = as.numeric(series_diarias_procesadas))

atipicos_revisados$atipico <- factor(
  ifelse(atipicos_revisados$X, "Removido", "Aceptado"),
  levels = c("Removido", "Aceptado")
)

Se toman las fechas marcadas con TRUE y en esas posiciones de la serie se asigna NA. No se borra la fila: se conserva el calendario continuo y el dato queda como faltante.
Se crea un factor con dos categorías para colorear el gráfico.

ggplot(series_diarias_procesadas, aes(fecha, valor)) +
  geom_line(color = "#77F") +
  geom_point(data = atipicos_revisados, aes(fecha, valor, color = atipico)) +
  scale_color_manual(values = c("#C33", "#FC3"), name = "Atípicos") +
  labs(
    x = "Fecha (días)",
    y = "mm",
    title = "Precipitación Diaria y Análisis de Atípicos",
    subtitle = "Estación Papayal"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

La línea azul corresponde a la serie ya depurada; los puntos rojos son los atípicos removidos (por lo que ya no forman parte de la línea) y los amarillos los atípicos aceptados como eventos reales.

7 Análisis de atípicos por mes

La lluvia tiene estacionalidad: en un clima bimodal o monomodal, un día de 60 mm puede ser normal en la temporada de lluvias, pero excepcional en la temporada seca. Un único umbral para todo el año puede pasar por alto atípicos de los meses secos y marcar de más en los meses húmedos. Por eso se repite el análisis mes a mes.

serie_atipicos_completa <- NULL

for (i in 1:12) {
  # Días con lluvia del mes i (de todos los años)
  serie_mes <- papayal_lluvia[month(index(papayal_lluvia)) == i]

  # plot = FALSE: calcula las estadísticas del boxplot sin dibujarlo
  atipicos <- boxplot(as.numeric(serie_mes), plot = FALSE)

  serie_atipicos <- serie_mes[serie_mes %in% atipicos$out]

  serie_atipicos_completa <- rbind(
    serie_atipicos_completa,
    data.frame(fecha = index(serie_atipicos),
               valor = as.numeric(serie_atipicos))
  )
}

nrow(serie_atipicos_completa)   # total de atípicos con el criterio mensual

## [1] 261

Para cada mes se filtran los días con lluvia de ese mes (en todos los años del registro), se calculan los límites propios de ese mes y se acumulan los atípicos en una sola tabla con rbind().

Primero comparamos la distribución mensual con un diagrama de caja por mes:

lluvia_mes <- data.frame(
  fecha = index(papayal_lluvia),
  valor = as.numeric(papayal_lluvia),
  mes   = factor(month(index(papayal_lluvia)), levels = 1:12,
                 labels = month.abb)
)

ggplot(lluvia_mes, aes(mes, valor)) +
  geom_boxplot(fill = "#CCF", outlier.colour = "#C33") +
  labs(x = "Mes", y = "mm/día",
       title = "Precipitación en días con lluvia, por mes",
       subtitle = "Estación Papayal") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

Y luego los atípicos mensuales sobre la serie completa:

p_papayal +
  geom_point(data = serie_atipicos_completa, aes(fecha, valor), color = "#C33") +
  labs(x = "Fecha", y = "mm",
       title = "Atípicos detectados mes a mes",
       subtitle = "Estación Papayal") +
  theme(plot.title = element_text(hjust = 0.5))

8 Rellenado de datos faltantes

Después de depurar los atípicos, si el proyecto lo requiere se hace rellenado de datos faltantes. Algunos análisis posteriores requieren una serie completa, por lo que se prueban dos métodos sencillos de rellenado.

sprintf("Hay %i faltantes en la serie original", sum(is.na(series_diarias_procesadas)))

## [1] "Hay 444 faltantes en la serie original"

8.1 Método 1: media móvil

Cada dato faltante se reemplaza por el promedio de una ventana de 7 días centrada en esa fecha. La idea es que la lluvia de un día se parece a la de los días cercanos.

serie_diaria_rellenada <- series_diarias_procesadas[,2]
serie_diaria_rellenada[is.na(serie_diaria_rellenada)] <- rollmean(serie_diaria_rellenada, 7, fill = "extend")

sprintf("Hay %i faltantes en la serie rellenada por media movil",
        sum(is.na(serie_diaria_rellenada)))

## [1] "Hay 0 faltantes en la serie rellenada por media movil"

¿Qué hace cada parte?

rollmean(x, 7, fill = "extend") calcula la media móvil con ventana de 7 datos. Con fill = "extend" los extremos de la serie (donde la ventana no cabe completa) se rellenan repitiendo el valor calculado más cercano, de modo que el resultado tenga la misma longitud que la serie original.
[is.na(serie_diaria_rellenada)] selecciona solo las posiciones donde faltaba el dato. Los datos observados no se tocan: la media móvil se usa únicamente para rellenar.
El último sprintf cuenta cuántos faltantes quedan después del rellenado.

8.2 Método 2: media general

Todos los faltantes se reemplazan por un único valor: el promedio de todos los datos disponibles de la serie.

serie_diaria_rellenada <- series_diarias_procesadas[,2]
serie_diaria_rellenada[is.na(serie_diaria_rellenada)] <- mean(serie_diaria_rellenada, na.rm = TRUE)

sprintf("Hay %i faltantes en la serie rellenada por media general",
        sum(is.na(serie_diaria_rellenada)))

## [1] "Hay 0 faltantes en la serie rellenada por media general"

Aquí mean(..., na.rm = TRUE) calcula el promedio ignorando los NA, y ese valor se asigna a todas las posiciones faltantes. Este método siempre deja la serie completa, pero tiene inconvenientes importantes:

Ignora la estacionalidad: asigna el mismo valor a un día de época seca y a uno de época lluviosa.
Como la serie diaria tiene muchos ceros, el promedio es un valor pequeño pero positivo, de modo que inventa lluvia en días que probablemente eran secos.
Reduce la variabilidad de la serie, lo que afecta estadísticos como la varianza y los extremos.

8.3 Método 3: Rellenado espacial con interpolación IDW

Cuando hay varias estaciones cercanas, un dato faltante en una de ellas se puede estimar a partir de lo que registraron las demás ese mismo día. A diferencia de la media móvil (que usa los días vecinos de la misma estación), aquí se usa la información de las estaciones vecinas en el mismo día.

8.3.1 Fundamento: ponderación por el inverso de la distancia (IDW)

En el método IDW (Inverse Distance Weighting) el valor estimado en un punto $x_0$ es un promedio ponderado de los valores observados $z_i$ en las estaciones con datos:

\[\hat{z}(x_0) = \frac{\sum_{i=1}^{n} w_i\, z_i}{\sum_{i=1}^{n} w_i}, \qquad w_i = \frac{1}{d_i^{\,p}}\]

donde $d_i$ es la distancia entre $x_0$ y la estación $i$, y $p$ es la potencia. Las estaciones más cercanas pesan más. Con $p$ grande la influencia de las estaciones lejanas se reduce rápidamente; con $p$ pequeño (como $p = 1$, usado aquí) la interpolación es más suave.

8.3.2 Paquetes y archivos

library(phylin)   # función idw()
library(terra)    # vect() y geom(): lectura de datos espaciales

Se necesitan dos archivos:

descargaDhime.csv: descarga de precipitación diaria de varias estaciones, con las columnas CodigoEstacion, Fecha y Valor.
Espacializacion.gpkg: un GeoPackage con una capa llamada Estaciones (puntos) que tiene un campo CODIGO con el código de cada estación.

8.3.3 Lectura y organización de los datos

datos <- read.csv("~/02. Profesionales/DOCENCIA/HIDROLOGIA-UNAL/PREPARACION/SIG/Espacializacion/descargaDhime.csv")
datos$Fecha <- as.Date(datos$Fecha)

cods <- unique(datos$CodigoEstacion)

series_diarias <- NULL
for (cod in cods) {
  prec_i <- datos[datos$CodigoEstacion == cod, c("Fecha", "Valor")]
  prec_i <- as.xts(prec_i$Valor, order.by = prec_i$Fecha)
  series_diarias <- cbind(series_diarias, prec_i)
}

colnames(series_diarias) <- cods

Aquí el archivo trae todas las estaciones juntas en formato largo (una fila por estación y fecha). El ciclo recorre cada código de estación, extrae su serie y la une con cbind(), que alinea las series por fecha, como se hizo al inicio del documento.

Igual que antes, se completa el calendario para tener un dato (o NA) por cada día:

fecha_inicial <- min(index(series_diarias))
fecha_final   <- max(index(series_diarias))
fechas        <- seq(fecha_inicial, fecha_final, by = "day")

series_diarias <- cbind(
  series_diarias,
  xts(order.by = fechas)   # serie vacía solo con el calendario completo
)

head(colnames(series_diarias))   # revisar cómo quedaron los nombres

## [1] "X16050170" "X16055010" "X16050060"

Nombres de columnas. Los códigos de estación son números, y cbind() sobre objetos xts convierte los nombres numéricos en nombres válidos de R anteponiendo una “X” (por ejemplo X21205010). El rellenado de más abajo depende de esto: con substring(2) se quita esa primera letra para recuperar el código original. Verifique con el head(colnames(...)) de arriba que sus nombres tienen el prefijo.

8.3.4 Diagnóstico inicial

summary(as.data.frame(series_diarias))

##    X16050170        X16055010        X16050060     
##  Min.   : 0.000   Min.   : 0.000   Min.   : 0.000  
##  1st Qu.: 0.000   1st Qu.: 0.000   1st Qu.: 0.000  
##  Median : 0.000   Median : 0.000   Median : 3.000  
##  Mean   : 4.512   Mean   : 2.961   Mean   : 7.134  
##  3rd Qu.: 2.275   3rd Qu.: 1.700   3rd Qu.: 9.000  
##  Max.   :70.500   Max.   :76.000   Max.   :79.000  
##  NA's   :31       NA's   :526      NA's   :2181

La fila NA's de cada columna indica cuántos días faltan por estación. Guarde estos números para compararlos con el resultado final.

8.3.5 Ubicación de las estaciones

estaciones <- vect("~/02. Profesionales/DOCENCIA/HIDROLOGIA-UNAL/PREPARACION/SIG/Espacializacion/Espacializacion.gpkg", "Estaciones")
estaciones

##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 3, 20  (geometries, attributes)
##  extent      : 4953414, 4960650, 2474399, 2486283  (xmin, xmax, ymin, ymax)
##  source      : Espacializacion.gpkg (Estaciones)
##  coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377) 
##  names       :    CODIGO          NOMBRE       CATEGORIA   TECNOLOGIA ESTADO
##  type        :     <num>           <chr>           <chr>        <chr>  <chr>
##  values      : 1.605e+07 RIO DE ORO [16~   Pluviométrica Convencional Activa
##                1.605e+07 BROTARE [16050~   Pluviométrica Convencional Activa
##                1.606e+07 AEROPUERTO AGU~ Climatológica ~ Convencional Activa
##  FECHA_INSTALACION FECHA_SUSPENSION ALTITUD LATITUD LONGITUD (and 10 more)
##              <chr>            <chr>   <int>   <num>    <num>              
##         15/08/1976               NA    1200   8.293   -73.39              
##         15/06/1960               NA    1545   8.401   -73.42              
##         15/03/1973               NA    1435   8.315   -73.36

vect() lee la capa Estaciones del GeoPackage como un objeto espacial (SpatVector). IDW necesita las coordenadas de cada estación para calcular distancias.

Sistema de coordenadas. Las distancias solo son interpretables si las coordenadas están en un sistema proyectado (en metros). Si la capa está en grados (latitud/longitud), las distancias quedan distorsionadas. Se puede revisar con crs(estaciones, describe = TRUE).

8.3.6 Rellenado día por día

La regla que se aplica a cada día (cada fila de la tabla) es:

Si no falta ninguna estación, no hay nada que hacer.
Si faltan dos o más estaciones, se descarta el día: con tan pocas estaciones con datos la interpolación sería poco confiable.
Si falta exactamente una estación, su valor se estima por IDW a partir de las demás.

for (i in 1:nrow(series_diarias)) {
  if (!any(is.na(series_diarias[i, ])))
    next
  else if (sum(is.na(series_diarias[i, ])) >= 2)
    next

  # Códigos de las estaciones con dato y de la estación sin dato en este día
  tiene_datos <- colnames(series_diarias)[!is.na(series_diarias[i, ])] |>
    substring(2)
  sin_datos <- colnames(series_diarias)[is.na(series_diarias[i, ])] |>
    substring(2)

  # Valores observados ese día
  datos_dia <- as.data.frame(series_diarias[i, !is.na(series_diarias[i, ])])

  # Coordenadas, en el mismo orden que los valores
  estaciones_datos <- geom(estaciones[match(tiene_datos, estaciones$CODIGO), ])[, c("x", "y")]

  estaciones_na <- geom(estaciones[match(sin_datos, estaciones$CODIGO), ])[, c("x", "y")] |>
    as.matrix() |>
    t()

  datos_interpolados <- idw(datos_dia, estaciones_datos, estaciones_na, p = 1) # Aplicación IDW

  series_diarias[i, is.na(series_diarias[i, ])] <- as.numeric(datos_interpolados)
}