Anexo_2

Seccion 2.1

Mapa de interpolacion con IDW

library(gstat)        
library(sp)           
library(raster)

Magdal <- shapefile("C:/Users/vale_/Desktop/UNAL/6to semestre/GB/magda2.shp")

Precipitacion <- shapefile("C:/Users/vale_/Desktop/UNAL/6to semestre/GB/prec2.shp")

Precipitacion@bbox <- Magdal@bbox

grd              <- as.data.frame(spsample(Precipitacion, "regular", n=319000))
names(grd)       <- c("X", "y")                                          
coordinates(grd) <- c("X", "y")                                     
gridded(grd)     <- TRUE                                            
fullgrid(grd)    <- TRUE                                            


proj4string(grd) <- proj4string(Precipitacion)                            


P.idw <- gstat::idw(rainfall ~ 1, Precipitacion, newdata=grd, idp=2.0)                                                        


r          <- raster(P.idw)

r

class      : RasterLayer 
dimensions : 744, 429, 319176  (nrow, ncol, ncell)
resolution : 358.7203, 358.7203  (x, y)
extent     : 904756.3, 1058647, 1479931, 1746819  (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 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
source     : memory
names      : var1.pred 
values     : 9.687549e-05, 66.17632  (min, max)

Magdal

class       : SpatialPolygonsDataFrame 
features    : 30 
extent      : 904750.7, 1058591, 1480000, 1746829  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 2
names       :         MUNIC, CODIGO 
min values  :     ALGARROBO,  47001 
max values  : ZONA BANANERA,  47980 

r.m <- raster::mask(r, Magdal)                                  

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

r.m

class      : RasterLayer 
dimensions : 744, 429, 319176  (nrow, ncol, ncell)
resolution : 358.7203, 358.7203  (x, y)
extent     : 904756.3, 1058647, 1479931, 1746819  (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 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
source     : memory
names      : var1.pred 
values     : 9.687549e-05, 66.17632  (min, max)

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.002) +                
  tm_legend(legend.outside=TRUE)

Seccion 2.2

Una mejor visualizacion del mapa

library(leaflet)                                                
library(RColorBrewer)                                            
pale <- colorNumeric(c("red", "orange", "yellow", "blue", "darkblue"), values(precip.mask),                                
                     na.color = "transparent")                  


leaflet() %>% addTiles() %>%                                     
  addRasterImage(r.m, colors = pale, opacity = 0.6) %>%         
  addLegend(pal = pale, values = values(r.m),                   
            title = "IDW Precipitación interpolada en Magdalena\n desde 03.06. a 07.06. en 2020 [mm]")

Seccion 2.3

Validacion de la interpolacion

Precipitacion

class       : SpatialPointsDataFrame 
features    : 772 
extent      : 904750.7, 1058591, 1480000, 1746829  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +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,             rain, rainfall 
min values  :     ALGARROBO,  47001,                0,        0 
max values  : ZONA BANANERA,  47980, 66.4547348022461,     66.5 

P <- Precipitacion                

IDW.out <- vector(length = length(P))

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

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

  abline(0,1)

par(OP)

Error Cuadratico Medio (RMSE)

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

[1] 5.129825

Validacion Cruzada/ Intervalos de Confianza (IC)

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

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

Zj <- as.matrix(apply(Zi, 1, sum, na.rm=T)/ n )                 

c1 <- apply(Zi,2,'-',Zj)                                        
c1 <- apply(c1^2, 1, sum, na.rm=T)                              

CI <- sqrt(1/(n*(n-1)) * c1)                                    


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

Seccion 2.4

Mapa

r <- raster(img.sig, layer="v")                                 
r.m <- raster::mask(r,Magdal)                                   


tm_shape(r.m) + tm_raster(n=7, title = "IDW\n95% intervalo de confianza \n(en mm)") +                                         
  tm_shape(P) + tm_dots(size = 0.09) +                           
  
  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 10 x64 (build 18362)

Matrix products: default

locale:
[1] LC_COLLATE=Spanish_Colombia.1252 
[2] LC_CTYPE=Spanish_Colombia.1252   
[3] LC_MONETARY=Spanish_Colombia.1252
[4] LC_NUMERIC=C                     
[5] LC_TIME=Spanish_Colombia.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets 
[6] methods   base     

other attached packages:
 [1] forcats_0.5.0      stringr_1.4.0      dplyr_1.0.0       
 [4] purrr_0.3.4        readr_1.3.1        tidyr_1.1.0       
 [7] tibble_3.0.1       ggplot2_3.3.1      tidyverse_1.3.0   
[10] sf_0.9-3           rgdal_1.5-8        RColorBrewer_1.1-2
[13] leaflet_2.0.3      tmap_3.0           raster_3.1-5      
[16] sp_1.4-2           gstat_2.0-6       

loaded via a namespace (and not attached):
 [1] nlme_3.1-144       fs_1.4.1           xts_0.12-0        
 [4] lubridate_1.7.8    httr_1.4.1         tools_3.6.3       
 [7] backports_1.1.7    R6_2.4.1           KernSmooth_2.23-16
[10] DBI_1.1.0          colorspace_1.4-1   withr_2.2.0       
[13] tidyselect_1.1.0   compiler_3.6.3     leafem_0.1.1      
[16] cli_2.0.2          rvest_0.3.5        xml2_1.2.5        
[19] scales_1.1.1       classInt_0.4-3     digest_0.6.25     
[22] rmarkdown_2.1      base64enc_0.1-3    dichromat_2.0-0   
[25] pkgconfig_2.0.3    htmltools_0.4.0    dbplyr_1.4.4      
[28] htmlwidgets_1.5.1  rlang_0.4.6        readxl_1.3.1      
[31] rstudioapi_0.11    FNN_1.1.3          generics_0.0.2    
[34] farver_2.0.3       zoo_1.8-8          jsonlite_1.6.1    
[37] crosstalk_1.1.0.1  magrittr_1.5       Rcpp_1.0.4.6      
[40] munsell_0.5.0      fansi_0.4.1        abind_1.4-5       
[43] lifecycle_0.2.0    stringi_1.4.6      leafsync_0.1.0    
[46] yaml_2.2.1         tmaptools_3.0      grid_3.6.3        
[49] blob_1.2.1         parallel_3.6.3     crayon_1.3.4      
[52] lattice_0.20-38    stars_0.4-1        haven_2.3.0       
[55] hms_0.5.3          knitr_1.28         pillar_1.4.4      
[58] spacetime_1.2-3    codetools_0.2-16   reprex_0.3.0      
[61] XML_3.99-0.3       glue_1.4.1         evaluate_0.14     
[64] modelr_0.1.8       png_0.1-7          vctrs_0.3.0       
[67] cellranger_1.1.0   gtable_0.3.0       assertthat_0.2.1  
[70] xfun_0.14          lwgeom_0.2-4       broom_0.5.6       
[73] e1071_1.7-3        rsconnect_0.8.16   class_7.3-15      
[76] viridisLite_0.3.0  intervals_0.15.2   units_0.6-6       
[79] ellipsis_0.3.1