Anexo 3: Kriging y IDW

rm(list=ls())

library(sf)
library(rgdal)
library(rasterVis)
library(dplyr)
library(raster)
library(gstat)
library(geoR)
library(berryFunctions)

(valle <- st_read("C:/Users/LUISA CARRION/Downloads/76_VALLE_DEL_CAUCA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

Reading layer `MGN_MPIO_POLITICO' from data source 
  `C:\Users\LUISA CARRION\Downloads\76_VALLE_DEL_CAUCA\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 42 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -77.54977 ymin: 3.091239 xmax: -75.70724 ymax: 5.047394
Geodetic CRS:  WGS 84
Simple feature collection with 42 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -77.54977 ymin: 3.091239 xmax: -75.70724 ymax: 5.047394
Geodetic CRS:  WGS 84
First 10 features:
   DPTO_CCDGO MPIO_CCDGO   MPIO_CNMBR
1          76      76001         CALI
2          76      76036    ANDALUCÍA
3          76      76041 ANSERMANUEVO
4          76      76054      ARGELIA
5          76      76111         BUGA
6          76      76113 BUGALAGRANDE
7          76      76122   CAICEDONIA
8          76      76126       CALIMA
9          76      76130   CANDELARIA
10         76      76147      CARTAGO
                             MPIO_CRSLC MPIO_NAREA MPIO_NANO
1                                  1536  563.04276      2017
2      Ordenanza 38 de Abril 25 de 1921  110.46042      2017
3                  Ordenanza 29 de 1925  305.45118      2017
4                  Ordenanza 15 de 1956   90.79604      2017
5                                  1555  825.86513      2017
6                                  1791  396.78132      2017
7  Ordenanza 21 del 20 de Abril de 1923  166.98369      2017
8      Ordenanza 49 de Junio 23 de 1939  793.49323      2017
9                                  1797  296.46056      2017
10                                 1863  248.16005      2017
        DPTO_CNMBR Shape_Leng  Shape_Area                       geometry
1  VALLE DEL CAUCA  1.1636697 0.045733211 MULTIPOLYGON (((-76.59175 3...
2  VALLE DEL CAUCA  0.6651004 0.008984851 MULTIPOLYGON (((-76.22406 4...
3  VALLE DEL CAUCA  0.9460442 0.024871077 MULTIPOLYGON (((-76.01558 4...
4  VALLE DEL CAUCA  0.4538261 0.007391017 MULTIPOLYGON (((-76.14316 4...
5  VALLE DEL CAUCA  2.0063716 0.067157337 MULTIPOLYGON (((-76.31608 3...
6  VALLE DEL CAUCA  1.0336063 0.032279294 MULTIPOLYGON (((-76.15131 4...
7  VALLE DEL CAUCA  0.6465900 0.013590500 MULTIPOLYGON (((-75.8539 4....
8  VALLE DEL CAUCA  1.5441977 0.064484981 MULTIPOLYGON (((-76.51747 4...
9  VALLE DEL CAUCA  0.8709866 0.024086066 MULTIPOLYGON (((-76.30455 3...
10 VALLE DEL CAUCA  0.8955580 0.020211481 MULTIPOLYGON (((-75.94518 4...

(precip <- raster("C:/Users/LUISA CARRION/Downloads/chirps-v2.0.2021.03.6.tif"))

class      : RasterLayer 
dimensions : 2000, 7200, 14400000  (nrow, ncol, ncell)
resolution : 0.05, 0.05  (x, y)
extent     : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : chirps-v2.0.2021.03.6.tif 
names      : chirps.v2.0.2021.03.6 

names(precip) <- "X2019.07"

(precip.mask <- mask(x = precip, mask = valle))

class      : RasterLayer 
dimensions : 2000, 7200, 14400000  (nrow, ncol, ncell)
resolution : 0.05, 0.05  (x, y)
extent     : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : X2019.07 
values     : 17.53901, 182.8078  (min, max)

(precip.points <- rasterToPoints(precip.mask, spatial = TRUE))

class       : SpatialPointsDataFrame 
features    : 677 
extent      : -77.475, -75.725, 3.124999, 5.024999  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :         X2019.07 
min values  : 17.5390129089355 
max values  : 182.807846069336 

names(precip.points) <- "rain"

(p.sf.magna <- st_transform(st_as_sf(precip.points), crs=32618))

Simple feature collection with 677 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 224942.3 ymin: 345532.7 xmax: 419590.3 ymax: 555525.6
Projected CRS: WGS 84 / UTM zone 18N
First 10 features:
       rain                  geometry
1  59.29984 POINT (380829.3 555525.6)
2  65.32109 POINT (380820.3 549997.7)
3  68.39664 POINT (380811.4 544469.8)
4  65.94273 POINT (386355.6 544461.1)
5  64.68098   POINT (375257.7 538951)
6  66.65178   POINT (380802.5 538942)
7  64.36810 POINT (386347.2 538933.3)
8  71.65340 POINT (369703.3 533432.4)
9  58.66259 POINT (375248.5 533423.1)
10 63.25397 POINT (380793.7 533414.1)

(valle.sf.magna <- st_transform(valle, crs=32618))

Simple feature collection with 42 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 216628.4 ymin: 341813.7 xmax: 421558.5 ymax: 558001.9
Projected CRS: WGS 84 / UTM zone 18N
First 10 features:
   DPTO_CCDGO MPIO_CCDGO   MPIO_CNMBR
1          76      76001         CALI
2          76      76036    ANDALUCÍA
3          76      76041 ANSERMANUEVO
4          76      76054      ARGELIA
5          76      76111         BUGA
6          76      76113 BUGALAGRANDE
7          76      76122   CAICEDONIA
8          76      76126       CALIMA
9          76      76130   CANDELARIA
10         76      76147      CARTAGO
                             MPIO_CRSLC MPIO_NAREA MPIO_NANO
1                                  1536  563.04276      2017
2      Ordenanza 38 de Abril 25 de 1921  110.46042      2017
3                  Ordenanza 29 de 1925  305.45118      2017
4                  Ordenanza 15 de 1956   90.79604      2017
5                                  1555  825.86513      2017
6                                  1791  396.78132      2017
7  Ordenanza 21 del 20 de Abril de 1923  166.98369      2017
8      Ordenanza 49 de Junio 23 de 1939  793.49323      2017
9                                  1797  296.46056      2017
10                                 1863  248.16005      2017
        DPTO_CNMBR Shape_Leng  Shape_Area                       geometry
1  VALLE DEL CAUCA  1.1636697 0.045733211 MULTIPOLYGON (((323192.2 39...
2  VALLE DEL CAUCA  0.6651004 0.008984851 MULTIPOLYGON (((364148.2 46...
3  VALLE DEL CAUCA  0.9460442 0.024871077 MULTIPOLYGON (((387393.9 54...
4  VALLE DEL CAUCA  0.4538261 0.007391017 MULTIPOLYGON (((373225.2 52...
5  VALLE DEL CAUCA  2.0063716 0.067157337 MULTIPOLYGON (((353893.4 44...
6  VALLE DEL CAUCA  1.0336063 0.032279294 MULTIPOLYGON (((372238.9 47...
7  VALLE DEL CAUCA  0.6465900 0.013590500 MULTIPOLYGON (((405258.7 48...
8  VALLE DEL CAUCA  1.5441977 0.064484981 MULTIPOLYGON (((331546.4 44...
9  VALLE DEL CAUCA  0.8709866 0.024086066 MULTIPOLYGON (((355088.2 38...
10 VALLE DEL CAUCA  0.8955580 0.020211481 MULTIPOLYGON (((395186.3 53...

precip2 <- as(p.sf.magna, 'Spatial')

valle2 <- as(valle.sf.magna, 'Spatial')

precip2@bbox <- valle2@bbox

train_index <- sample(1:nrow(precip2), 0.7 * nrow(precip2))
test_index <- setdiff(1:nrow(precip2), train_index)
ptos_train <- precip2[train_index, ]
ptos_test  <- precip2[test_index,]

ptos_train

class       : SpatialPointsDataFrame 
features    : 473 
extent      : 224942.3, 419590.3, 345549.4, 549997.7  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
variables   : 1
names       :             rain 
min values  : 18.0643348693848 
max values  : 182.807846069336 

Interpolación por IDW

library(gstat)

grd   <- as.data.frame(spsample(precip2, "regular", n=319000))

names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE  
fullgrid(grd)    <- TRUE

x <- CRS(SRS_string='EPSG:32618')

wkt <- comment(x)

(y <- CRS(SRS_string = wkt))

CRS arguments:
 +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 

sp::proj4string(grd) <- y

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

[inverse distance weighted interpolation]

r <- raster(P.idw)

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

class      : RasterLayer 
dimensions : 580, 550, 319000  (nrow, ncol, ncell)
resolution : 372.6692, 372.6692  (x, y)
extent     : 216727.2, 421695.3, 341814.9, 557963.1  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
source     : memory
names      : var1.pred 
values     : 18.29223, 182.779  (min, max)

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  [mm]")

Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs
Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum World Geodetic System 1984 in Proj4 definition
Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs
Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum World Geodetic System 1984 in Proj4 definition

Interpolacioón por Kriging

p.sf.magna

Simple feature collection with 677 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 224942.3 ymin: 345532.7 xmax: 419590.3 ymax: 555525.6
Projected CRS: WGS 84 / UTM zone 18N
First 10 features:
       rain                  geometry
1  59.29984 POINT (380829.3 555525.6)
2  65.32109 POINT (380820.3 549997.7)
3  68.39664 POINT (380811.4 544469.8)
4  65.94273 POINT (386355.6 544461.1)
5  64.68098   POINT (375257.7 538951)
6  66.65178   POINT (380802.5 538942)
7  64.36810 POINT (386347.2 538933.3)
8  71.65340 POINT (369703.3 533432.4)
9  58.66259 POINT (375248.5 533423.1)
10 63.25397 POINT (380793.7 533414.1)

(the_crs <- crs(p.sf.magna))

CRS arguments:
 +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 

p.sf.magna$x <- st_coordinates(p.sf.magna)[,1]
p.sf.magna$y <- st_coordinates(p.sf.magna)[,2]

geoprec <- as.geodata(cbind(p.sf.magna$x,p.sf.magna$y,p.sf.magna$rain)) 
vario <- variog(geoprec, max.dist=300000)

variog: computing omnidirectional variogram

fit <- variofit(vario)

variofit: covariance model used is matern 
variofit: weights used: npairs 
variofit: minimisation function used: optim 

Warning in variofit(vario) :
  initial values not provided - running the default search

variofit: searching for best initial value ... selected values:
              sigmasq   phi         tausq kappa
initial.value "3762.34" "155430.18" "0"   "0.5"
status        "est"     "est"       "est" "fix"
loss value: 17963712484.39 

plot(vario) ; lines(fit)

d <- sapply(1:nrow(p.sf.magna), function(i) min(berryFunctions::distance(
                          p.sf.magna$x[i], p.sf.magna$y[i], p.sf.magna$x[-i], p.sf.magna$y[-i])))
hist(d/1000, breaks=20, main="distance to closest gauge [km]")

res <- 1000 # 0.5 km, since stations are 4.6 km apart on average 
grid <- expand.grid(seq(min(p.sf.magna$x),max(p.sf.magna$x),res),
                    seq(min(p.sf.magna$y),max(p.sf.magna$y),res))
krico <- krige.control(type.krige="OK", obj.model=fit)
krobj <- krige.conv(geoprec, locations=grid, krige=krico)

krige.conv: model with constant mean
krige.conv: Kriging performed using global neighbourhood 

grid_sf <- sf::st_as_sf(grid, coords=1:2, crs=sf::st_crs(valle.sf.magna))
isinp <- sapply(sf::st_within(grid_sf, valle.sf.magna), length) > 0
krobj2 <- krobj
krobj2$predict[!isinp] <- NA

(k_raster <- raster::rasterFromXYZ(cbind(grid, temp=krobj2$predict),
crs=y))

class      : RasterLayer 
dimensions : 210, 195, 40950  (nrow, ncol, ncell)
resolution : 1000, 1000  (x, y)
extent     : 224442.3, 419442.3, 345032.7, 555032.7  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
source     : memory
names      : temp 
values     : 17.16097, 182.7659  (min, max)

leaflet() %>% addTiles() %>%
  addRasterImage(k_raster$temp, colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(r.m$var1.pred),
    title = "Precipitaciones abril 2021 [mm]")

Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs
Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum World Geodetic System 1984 in Proj4 definition
Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs
Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum World Geodetic System 1984 in Proj4 definition
Warning in colors(.) :
  Some values were outside the color scale and will be treated as NA

P <- ptos_train
IDW.out <- vector(length = length(P))
for (i in 1:length(P)) {
  IDW.out[i] <- gstat::idw(rain ~ 1, P[-i,], P[i,], idp=2.0)$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]
[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]
[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]
[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]
[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]
[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]

# Plot the differences
OP <- par(pty="s", mar=c(4,3,0,0))
  plot(IDW.out ~ P$rain, asp=1, xlab="Observed", ylab="Predicted", pch=16,
       col=rgb(0,0,0,0.5))
  abline(lm(IDW.out ~ P$rain), col="red", lw=2,lty=2)

  abline(0,1)

par(OP)

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

[1] 10.75762

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

class       : SpatialPointsDataFrame 
features    : 47 
extent      : 241621.9, 414044.5, 345558.1, 533389.7  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
variables   : 1
names       :             rain 
min values  : 21.0774211883545 
max values  : 160.437210083008 

img <- gstat::idw(rain~1, P, newdata=grd, idp=3.0)

[inverse distance weighted interpolation]

n   <- length(P)
Zi  <- matrix(nrow = length(img$var1.pred), ncol = n)

st <- stack()
for (i in 1:n){
  Z1 <- gstat::idw(rain~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}

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 

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

tm_shape(r.m) + tm_raster(n=7,title="IDW\n95% confidence interval \n(in mm)") +
  tm_shape(P) + tm_dots(size=0.02) +
  
  tm_legend(legend.outside=TRUE)