Interpolación IDW para datos de precipitación en Vichada

1. Introducción

Este cuaderno corresponde al tecer anexo del informe final correspondiente a la asignatura Geomática Básica, en el se ilustra el metodo de la distancia inversa ponderada para los datos de precipitación

2. Uso del método de la distancia inversa ponderada para la interpolación de los datos de precipitación

Primero, se crea una cuadricula vacia , con un total de 100000 celdas y se le asigna información de proyección

grd              <- as.data.frame(spsample(ptos_train, "regular", n=100000))
names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE  
fullgrid(grd)    <- TRUE  
proj4string(grd) <- proj4string(ptos_train)

Posteriormente, se realiza la terea de interpolación con ayuda de la función idw de la libreria gstat .Para la interpolación de los datos , será empleada una potencia de 3 (idp=3.0) , luego el resultado de esta interpolación se convierte a un objeto ráster y se recorta a la extensión de la zona de estudio (En este caso a la extensión del objeto “Vichada 2”)

P.idw <- gstat::idw(precipitacion ~ 1, ptos_train, newdata=grd, idp=3.0)

[inverse distance weighted interpolation]

r       <- raster(P.idw)
(r.m   <- raster::mask(r, Vichada2))

class      : RasterLayer 
dimensions : 313, 320, 100160  (nrow, ncol, ncell)
resolution : 1249.511, 1249.511  (x, y)
extent     : 1333284, 1733127, 799447.3, 1190544  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
source     : memory
names      : var1.pred 
values     : 0.5035295, 84.58405  (min, max)

Ahora, se plotea el resultado de la interpolación.Primero se realiza el gráfico con tmap y posteriormente se hace uso de la libreria leaflet para lograr una mejor visualización

tm_shape(r.m) + tm_raster(n=10,palette = "RdBu", midpoint = 43, title="Interpolación IDW \n para precipitación \n(mm)") + 
  tm_shape(ptrain) + tm_dots(size = 0.03, col = "#c4f0f5") +
  tm_legend(legend.outside=TRUE, title.size=1.2, text.size= 0.8)

library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("#B2182B", "#FFE4A9", "#afc9a9", "#86bbdb", "#2e0091"), values(r.m),
  na.color = "transparent")

leaflet() %>% addTiles() %>%addRasterImage(r.m, colors = pal, opacity = 0.6) %>%addLegend(pal = pal, values = values(r.m), title = "Interpolación IDW para precipitación en Vichada del 26.04 al 30.04 en 2020 [mm]")

2.1 Validación Cruzada

Ahora, se realizará una tarea de validación (Leave-one-out cross-validation) para medir el error en los valores interpolados y usar esta información para perfeccionar la tarea de interopolación , pues , a través de esta tarea de validación se puede conocer cual potencia se debe usar para tener una mejor aproximación a los datos

Esta tarea de validación ,consiste en extraer un objeto de la muestra y usar los valores de los objetos restantes para tratar de estimar su valor .Este proceso se realiza n veces , es decir se realiza una vez por cada objeto de la muestra.

P <- ptrain
IDW.out <- vector(length = length(P))
for (i in 1:length(P)) {IDW.out[i] <- gstat::idw(precipitacion ~ 1, P[-i,], P[i,], idp=3.0)$var1.pred}

El siguiente gráfico , corresponde a un diagrama de dispersión entre los datos observados, es decir los valores reales del punto que fue omitido en cada caso y los datos de precipitación pronósticados o estimados

OP <- par(pty="s", mar=c(4,3,0,0))
  plot(IDW.out ~ P$precipitacion, asp=1, xlab="Observado", ylab="Pronosticado", pch=20,col=rgb(0,0,0,0.5))
abline(lm(IDW.out ~ P$precipitacion), col="red", lw=2,lty=3)

abline(0,1)

También , se puede elegir la potencia en función del error cuadratico medio, pues aquella que tiene un valor menor en este parámetro se considera un mejor estimador ,por tanto , los valores estimados estarán mas cercanos a los valores observados

sqrt( sum((IDW.out - P$precipitacion)^2) / length(P))

[1] 2.234104

2.2 Intervalos de Confianza para la interpolación IDW

Otra medida que puede ser empleada para conocer la presición del metodo de interpolación aplicado es el cálculo de los intervalos de confianza con ayuda de la técnica de validación cruzada para cada pixel. Si, el rango de valores es muy alto no podemos estar seguros del valor que fue estimado, en cambio si el rango de valores para el intervalo de confianza es mas bajo se tendra mas seguridad del valor obtenido en la interpolación pues nos indica que este ultimo no es tan sensible a la presencia o ausencia de un punto de muestra

Este proceso se ilustrará para el 5% de los datos , pues el cálculo de intervalo de confianza para cada valor estimado puede tardarse mucho

index = sample(1:nrow(ptos_train), 0.05 * nrow(ptos_train))
(P <- ptos_train[index, ])

class       : SpatialPointsDataFrame 
features    : 81 
extent      : 1339270, 1695656, 821295.6, 1183910  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 4
names       :         MUNIC, CODIGO,            lluvia, precipitacion 
min values  :      CUMARIBO,  99001, 0.674332559108734,           0.7 
max values  : SANTA ROSALÍA,  99773,  80.6951599121094,          80.7 

Se crea la superficie interpolada

img <- gstat::idw(precipitacion~1, P, newdata=grd, idp=3.0)
n   <- length(P)
Zi  <- matrix(nrow = length(img$var1.pred), ncol = n)

Luego,se aplica el método de validación cruzada

st <- stack()
for (i in 1:n){
  Z1 <- gstat::idw(precipitacion~1, P[-i,], newdata=grd, idp=3.0)
  st <- addLayer(st,raster(Z1,layer=1))   
  Zi[,i] <- n * img$var1.pred - (n-1) * Z1$var1.pred}

Se calcula el intervalo de confianza

Se crea un raster a partir del intervalo de confianza y la predicción de la variable

img.sig   <- img
img.sig$v <- CI /img$var1.pred 

Posteriormente , se recorta este ráster a la extensión del departamento de Vichada

r1 <- raster(img.sig, layer="v")
r.m1 <- raster::mask(r1, Vichada2)

Se plotea el mapa correspondiente a los intervalos de confianza

tm_shape(r.m1) + tm_raster(palette = "Reds" , n=7,title="IDW\nIntervalo de confianza (95%) \n(mm)") +
  tm_shape(P) + tm_dots(size=0.02) +
  
  tm_legend(legend.outside=TRUE)

CRS object has comment, which is lost in output