En este cuaderno presenta ejemplos para el rellenado de series de datos, se utilizan valores de precipitación diaria en la estación Nataga en el departamento del Huila, Colombia. Las metodologías utilizadas para el rellenado son media móvil y métodos espaciales, específicamente el Inverso Ponderado de la Distancia (IDW, por sus siglas en ingles). La información necesaria para ejecutar este procedimiento se encuentra en este enlace.

El código se recomienda iniciarlo limpiando el espacio de trabajo, llamando todas las librerías necesarias y seleccionando el directorio de trabajo.

rm(list = ls())
library(xts) # Series Temporales
library(lubridate) # Utilidades para manejo de fechas
library(terra) # Información Geográfica
library(raster)
library(phylin) # IDW

setwd("~/Diego/02.Profesionales/UPJ/DATOS/Rellenado & Inteterpolacion/")

1 Media Movil

La media móvil es útil cuando se cuentan con ventanas cortas sin información, desde el punto de vista de su implementación cuenta con la ventaja de su sencillez y facilidad de aplicar tanto en lenguajes de programación como en herramientas como Excel.

Se inicia cargando la información que se desea completar, es recomendable trabajarla como serie de tiempo, convirtiéndola en xts, sin embargo, no es necesario para este proceso. Posteriormente se hace uso de la función rollapply. La ventana de tiempo seleccionada es a criterio del profesional, para precipitación 7 días es una buena opción. A la función se le especifica que aplique el promedio (mean), que omita los valores vacios (na.rm = T) y la indicación fill = c(“extend”, NA) reduce los errores en los extremos de la serie donde no se puede aplicar la ventana de 7 días. Posterior a aplicar la media móvil, se rellanan las casillas con datos faltantes de la serie original con estos nuevos datos.

Cuando existen periodos consecutivos de datos faltantes mas largos que la ventana definida no se podrá rellenar. En este caso se puede combinar dos metodos de rellenado como los que se verán en la siguiente sección.

datos <- read.csv("Nataga.csv") # Lee los datos
head(datos) # Los muestra para ver como está el formato 
##        Fecha Valor
## 1 1992-01-01     0
## 2 1992-01-02     0
## 3 1992-01-03     0
## 4 1992-01-04     3
## 5 1992-01-05     0
## 6 1992-01-06    15
datos <- as.xts(datos[,2], order.by = as.Date(datos[,1])) # Los convierte en serie temporal

datos.rellenados <- rollapply(datos, 7, mean, na.rm = T, fill = c("extend", NA)) # Aplica la media movil

datos[is.na(datos)] <- datos.rellenados[is.na(datos)] # Extrae los valores de la media movil para rellenar las casillas 'NA' de los datos originales
cat(sprintf("\nLos valores faltantes antes de la media movil eran %i y despues son %i.\n", sum(is.na(datos)), sum(is.na(datos.rellenados))))
## 
## Los valores faltantes antes de la media movil eran 6 y despues son 6.

2 Inverso Ponderado de la Distancia (IDW)

Para este ejemplo se utilizarán datos sintéticos, debido a que no se tiene información de otras estaciones. Se utilizarán cinco estaciones, Nataga y cuatro que la rodean para la interpolación, estas estaciones se encuentran en un archivo tipo shapefile, llamado Estaciones.shp, para leerlo se utiliza la función vect del paquete terra.

En los datos sintéticos se asumió que la estación Natagá no presentaba registro por lo que se dividieron las estaciones entre las que tienen valores y las que no. Mediante la función geom se extraen las coordenadas tanto de las estaciones con las que se va a interpolar como a la que se va a interpolar. Luego se emplea la función idw en donde inicialmente se le definen los valores conocidos, luego las coordenadas asociadas a estos valores y finalmente las coordenadas a las que se desea interpolar.

estaciones <- vect("Estaciones.shp") # Lee el shp de estaciones
as.data.frame(estaciones[,2:5]) # Muestra un extrato de la tabla de atributos
##     CODIGO                       nombre     CATEGORIA                TECNOLOGIA
## 1 21057090     PUENTE ITAIBE [21057090]  Limnimétrica              Convencional
## 2 21050290         TESALIA 2 [21050290] Pluviométrica              Convencional
## 3 21050090     NATAGA  - AUT [21050090] Pluviométrica Automática con Telemetría
## 4 21050220          SAN LUIS [21050220] Pluviométrica              Convencional
## 5 21080110 SAN JOSE HACIENDA [21080110] Pluviográfica              Convencional
plot(estaciones, 2) # Plotea la ubicación de las estaciones

datos.sinteticos <- c(22, 17, NA, 16, 9) # Datos sintéticos 

stn.valores <- terra::geom(estaciones[!is.na(datos.sinteticos),])[,c("x", "y")] # Coordenadas de estaciones con valores
stn.valores <- as.data.frame(stn.valores)

stn.na <- terra::geom(estaciones[is.na(datos.sinteticos),])[,c("x", "y")]  # Coordenadas de estaciones sin valores
stn.na <- data.frame(x = stn.na[1], y = stn.na[2])

datos.interpolados <- idw(datos.sinteticos[!is.na(datos.sinteticos)], stn.valores, stn.na, p = 1) # Aplicación IDW
cat(sprintf("\nEl valor interpolado en la estación faltante es %.1fmm\n", datos.interpolados[1])) 
## 
## El valor interpolado en la estación faltante es 16.0mm

Es análogo el proceso cuando se tiene una red de puntos y se desea crear un campo de precipitaciones. En vez de interpolar al punto de la estación faltantes se interpola a cada uno de los puntos de la red. Posteriormente se puede crear un raster con estas interpolaciones. Después de aplicar la interpolación se organiza un data frame con las coordenadas y los valores interpolados, después se indica que las columnas denominadas “x” y “y” corresponden con las coordenadas, se indica que se trata de una malla, se crea el raster y se aclara que sistemas de coordenadas utiliza.

datos.sinteticos <- c(22, 17, 16, 16, 9) # Datos sintéticos 

stn <- terra::geom(estaciones)[,c("x", "y")] # Coordenadas de las estaciones
stn <- as.data.frame(stn)

grid <- vect("Grid.shp") # shp de la red de puntos
grid <- terra::geom(grid)[,c("x", "y")] # coordenadas de los puntos
grid <- as.data.frame(grid)

datos.interpolados <- idw(datos.sinteticos, stn, grid, p = 1) # Aplicación IDW

grid <- data.frame(grid, valores = datos.interpolados[,1]) 
coordinates(grid) <- ~x+y # Se especifica que 'x' y 'y' son las coordenadas
gridded(grid) <- T # Se indica que se trata de una red de puntos
grid <- raster(grid, "valores") # Se crea el raster
projection(grid) <- crs("+init=epsg:3116") # Se le asigna un sistemas de coordenadas

writeRaster(grid, "CampoPrecipitación.tif", overwrite = T) # Se escribe el raster
plot(grid)