Anexo 4:Kriging

P$X <- coordinates(P)[,1]
P$Y <- coordinates(P)[,2]

f.1 <- as.formula(rain ~ X + Y) 
var.smpl <- variogram(f.1, P, cloud = FALSE, cutoff=100000, width=8990)
dat.fit  <- fit.variogram(var.smpl, fit.ranges = TRUE, fit.sills = TRUE,
                          vgm(psill=3, model="Mat", range=150000, nugget=0.0))

dat.fit

plot(var.smpl, dat.fit, xlim=c(0,130000))

f.1 <- as.formula(rain ~ X + Y) 
dat.krg <- krige( f.1, P, grd, dat.fit)

[using universal kriging]

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

r.m

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     : 20.9284, 168.282  (min, max)

tm_shape(r.m) + 
  tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE, 
            title="Universal Kriging\nPredicted precipitation \n(in mm)") +
  tm_shape(P) + tm_dots(size=0.02) +
  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.

r   <- raster(dat.krg, layer="var1.var")
r.m <- raster::mask(r, valle2)

tm_shape(r.m) + 
  tm_raster(n=7, palette ="Reds",
            title="Kriging Interpolation\nVariance map \n(in squared mm)") +tm_shape(P) + tm_dots(size=0.02) +
  tm_legend(legend.outside=TRUE)

r   <- sqrt(raster(dat.krg, layer="var1.var")) * 1.96
r.m <- raster::mask(r, valle2)

tm_shape(r.m) + 
  tm_raster(n=7, palette ="Reds",
            title="Kriging Interpolation\n95% CI map \n(in mm)") +tm_shape(P) + tm_dots(size=0.02) +
  tm_legend(legend.outside=TRUE)