En este cuaderno se muestra uno de los métodos de interpolación espacial aplicado a un conjunto de datos de precipitación correspondientes al departamento de Boyacá.

IDW

Muchos paquetes comparten los mismos nombres de funciones. Esto puede ser un problema cuando estos paquetes se cargan en una misma sesión R. Por ejemplo, la función de intersección está disponible en los paquetes base, spatstat y raster, todos los cuales se cargan en esta sesión actual. Para asegurarse de que se selecciona la función adecuada, es una buena idea introducir el nombre de la función con el nombre del paquete como en raster :: intersect ().

Este consejo se usará en el siguiente fragmento de código cuando se llame a la función idw que está disponible tanto en spatstat como en gstat. Tenga en cuenta que la función dirichlet (como la mayoría de las funciones en el paquete spatsat) requiere que el objeto de punto esté en formato ppp, de ahí la sintaxis en línea as.ppp (Precipitacion).

La salida de IDW es una raster. Esto requiere que primero creemos una cuadrícula ráster vacía, luego interpolemos los valores de precipitación a cada celda de la cuadrícula sin muestrear. Se utilizará un valor de potencia IDW de 2 (idp = 2.0).

library(gstat)# Use la rutina idw de gstat
library(sp)# Utilizado para la función spsample
library(raster)

Es importante subir el mapa de puntos de precipitacion y el archivo que represente nuestra área de interés:

Boyaca<-shapefile("C:/Users/Brian/Desktop/Daniela/geomatica/Informe final/boyaca2.shp")
Precipitacion <- shapefile('C:/Users/Brian/Desktop/Daniela/geomatica/Informe final/precipitacion2.shp')

Debemos asegurarnos de que las dos extensiones coincidan:

Precipitacion@bbox <-Boyaca@bbox
# # Crear una cuadrícula vacía donde n es el número total de celdas
grd              <- as.data.frame(spsample(Precipitacion, "regular", n=319000))
# # Necesita averiguar cuál es el tamaño esperado del grd de salida
names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE  # Crear objeto SpatialPixel
fullgrid(grd)    <- TRUE  # Crear objeto SpatialGrid

# Agregar información de proyección de Precipitacion a la cuadrícula vacía
proj4string(grd) <- proj4string(Precipitacion)

# Interpolar las celdas de la cuadrícula utilizando un valor de potencia de 2 (idp = 2.0)
P.idw <- gstat::idw(lluvia ~ 1, Precipitacion, newdata=grd, idp=2.0)
## [inverse distance weighted interpolation]
# Convertir a objeto ráster y luego recortar a AOI
r       <- raster(P.idw)
r
## class      : RasterLayer 
## dimensions : 531, 600, 318600  (nrow, ncol, ncell)
## resolution : 500.2415, 500.2415  (x, y)
## extent     : 934929.4, 1235074, 1006622, 1272251  (xmin, xmax, ymin, ymax)
## crs        : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## source     : memory
## names      : var1.pred 
## values     : 4.701615, 93.04158  (min, max)
Boyaca
## class       : SpatialPolygonsDataFrame 
## features    : 123 
## extent      : 934932.9, 1235253, 1006592, 1272398  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## variables   : 2
## names       :      MUNIC, CODIGO 
## min values  :    ALMEIDA,  15001 
## max values  : ZETAQUIRÁ,  15897
r.m   <- raster::mask(r, Boyaca)
r.m
## class      : RasterLayer 
## dimensions : 531, 600, 318600  (nrow, ncol, ncell)
## resolution : 500.2415, 500.2415  (x, y)
## extent     : 934929.4, 1235074, 1006622, 1272251  (xmin, xmax, ymin, ymax)
## crs        : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## source     : memory
## names      : var1.pred 
## values     : 4.701615, 93.04158  (min, max)

Ahora trazamos el mapa:

library(tmap)
tm_shape(r.m) + 
  tm_raster(n=10,palette = "RdBu", auto.palette.mapping = FALSE,
            title="Distancia inversa ponderada \nPrecipitación prevista \n(en mm)") + 
  tm_shape(Precipitacion) + tm_dots(size=0.2) +
  tm_legend(legend.outside=TRUE)
## Warning: The argument auto.palette.mapping is deprecated. Please use midpoint
## for numeric data and stretch.palette for categorical data to control the palette
## mapping.

Una mejor visualización de la superficie interpolada:

## class       : SpatialPolygonsDataFrame 
## features    : 123 
## extent      : -74.66496, -71.94885, 4.655196, 7.055557  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 9
## names       : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR,                         MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO, DPTO_CNMBR,     Shape_Leng,       Shape_Area 
## min values  :         15,      15001,    ALMEIDA,                               1537,   25.35271459,      2017,    BOYACÁ, 0.204497978002, 0.00206924119091 
## max values  :         15,      15897, ZETAQUIRÁ, Ordenanza 8 de Diciembre 4 de 1965, 1513.60104233,      2017,    BOYACÁ,  2.40391520946,   0.123609862522
library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("red", "orange", "yellow", "blue", "darkblue"), values(precipitacion.mask),
  na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(r.m, colors = pal, opacity = 0.6) %>%
  addLegend(pal = pal, values = values(r.m),
    title = "IDW Precipitación interpolada en Boyacá\n del 26.04 to 30.04 en 2020 [mm]")

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

La elección de la función de poder puede ser subjetiva. Para ajustar la elección del parámetro de potencia, puede realizar una rutina de validación de omisión, que 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, para medir el error en los valores interpolados.

Precipitacion
## class       : SpatialPointsDataFrame 
## features    : 83 
## extent      : 934932.9, 1235253, 1006592, 1272398  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## variables   : 3
## names       :          MUNIC, CODIGO,           lluvia 
## min values  :        ALMEIDA,  15022, 4.69908825556437 
## max values  : VILLA DE LEYVA,  15861, 93.0493799845378
P <- Precipitacion
IDW.out <- vector(length = length(P))
for (i in 1:length(P)) {IDW.out[i] <- gstat::idw(lluvia ~ 1, P[-i,], P[i,], idp=2)$var1.pred}
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]

Trazar las diferencias

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

  par(OP)

El error cuadrático medio (RMSE) se puede calcular desde IDW.out de la siguiente manera:

# Calcular RMSE
sqrt( sum((IDW.out - P$lluvia)^2) / length(P))
## [1] 10.52289

Intervalos de Confianza para la interpolación IDW

Además de generar una superficie interpolada, puede crear un mapa de intervalo de confianza del 95% del modelo de interpolación. Aquí crearemos un mapa de IC del 95% a partir de una interpolación IDW que utiliza un parámetro de potencia de 2 (idp = 2.0).

# Implementación de una técnica de navaja para estimar un intervalo de confianza en cada punto no muestreado.

# Crear la superficie interpolada
img <- gstat::idw(lluvia~1, P, newdata=grd, idp=2.0)
## [inverse distance weighted interpolation]
n   <- length(P)
Zi  <- matrix(nrow = length(img$var1.pred), ncol = n)

# Elimina un punto y luego interpola (haz esto n veces para cada punto)
st <- stack()
for (i in 1:n){
  Z1 <- gstat::idw(lluvia~1, P[-i,], newdata=grd, idp=2.0)
  st <- addLayer(st,raster(Z1,layer=1))
# Seudovalor Z calculado en j
  Zi[,i] <- n * img$var1.pred - (n-1) * Z1$var1.pred
}
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
#Estimador de navaja del parámetro Z en la ubicación j
Zj <- as.matrix(apply(Zi, 1, sum, na.rm=T) / n )

# Calcular (Zi* - Zj)^2
c1 <- apply(Zi,2,'-',Zj)# Calcular la diferencia
c1 <- apply(c1^2, 1, sum, na.rm=T ) # Suma el cuadrado de la diferencia

# Calcular el intervalo de confianza
CI <- sqrt( 1/(n*(n-1)) * c1)

# Crear ráster (CI / valor interpolado)
img.sig   <- img
img.sig$v <- CI /img$var1.pred 
# # Recorte el raster de confianza a AOI
r <- raster(img.sig, layer="v")
r.m <- raster::mask(r, Boyaca)

# Trazar el mapa
tm_shape(r.m) + tm_raster(n=7,title="IDW\nIntervalo de confianza del 95% \n(en mm)") +
  tm_shape(P) + tm_dots(size=0.2) +
  
  tm_legend(legend.outside=TRUE)

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 8.1 x64 (build 9600)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.1252  LC_CTYPE=Spanish_Colombia.1252   
## [3] LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] RColorBrewer_1.1-2 leaflet_2.0.3      tmap_3.0           raster_3.0-12     
## [5] sp_1.4-1           gstat_2.0-6       
## 
## loaded via a namespace (and not attached):
##  [1] zoo_1.8-7          xfun_0.12          sf_0.9-3           lattice_0.20-38   
##  [5] colorspace_1.4-1   htmltools_0.4.0    stars_0.4-1        viridisLite_0.3.0 
##  [9] spacetime_1.2-3    yaml_2.2.1         base64enc_0.1-3    XML_3.99-0.3      
## [13] rlang_0.4.5        e1071_1.7-3        DBI_1.1.0          lifecycle_0.2.0   
## [17] stringr_1.4.0      munsell_0.5.0      htmlwidgets_1.5.1  leafsync_0.1.0    
## [21] codetools_0.2-16   evaluate_0.14      knitr_1.28         crosstalk_1.1.0.1 
## [25] parallel_3.6.3     class_7.3-15       leafem_0.1.1       xts_0.12-0        
## [29] Rcpp_1.0.3         KernSmooth_2.23-16 scales_1.1.0       classInt_0.4-3    
## [33] lwgeom_0.2-1       jsonlite_1.6.1     abind_1.4-5        FNN_1.1.3         
## [37] farver_2.0.3       png_0.1-7          digest_0.6.25      stringi_1.4.6     
## [41] tmaptools_3.0      grid_3.6.3         rgdal_1.4-8        tools_3.6.3       
## [45] magrittr_1.5       dichromat_2.0-0    rmarkdown_2.1      R6_2.4.1          
## [49] units_0.6-5        intervals_0.15.2   compiler_3.6.3