Interpolación

Ana María Montaño Hernández

precip <- raster(("C:/Users/user/Documents/Intro_to_R/chirps-v2.0.2020.04.6.tif"))

precip

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 +ellps=WGS84 +towgs84=0,0,0 
source     : C:/Users/user/Documents/Intro_to_R/chirps-v2.0.2020.04.6.tif 
names      : chirps.v2.0.2020.04.6 

(aoi <- st_read("C:/Users/user/Documents/Norte de Santander/54_NORTE DE SANTANDER/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\user\Documents\Norte de Santander\54_NORTE DE SANTANDER\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
Simple feature collection with 40 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -73.63379 ymin: 6.872201 xmax: -72.04761 ymax: 9.290847
geographic CRS: WGS 84
Simple feature collection with 40 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -73.63379 ymin: 6.872201 xmax: -72.04761 ymax: 9.290847
geographic CRS: WGS 84
First 10 features:
   DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR                           MPIO_CRSLC
1          54      54001     CÚCUTA                                 1972
2          54      54003     ABREGO                                 1806
3          54      54051  ARBOLEDAS                                 1835
4          54      54099  BOCHALEMA                                 1826
5          54      54109 BUCARASICA                                 1838
6          54      54125     CÁCOTA                                 1630
7          54      54128    CÁCHIRA                                 1911
8          54      54172  CHINÁCOTA                                 1775
9          54      54223  CUCUTILLA                                 1834
10         54      54239    DURANIA Ordenanza 12 del 27 de Marzo de 1911
   MPIO_NAREA MPIO_NANO         DPTO_CNMBR Shape_Leng Shape_Area
1   1133.2031      2017 NORTE DE SANTANDER  2.8636253 0.09289252
2   1382.4498      2017 NORTE DE SANTANDER  2.2393759 0.11336240
3    456.1490      2017 NORTE DE SANTANDER  1.1611936 0.03736204
4    170.2669      2017 NORTE DE SANTANDER  0.7565329 0.01394455
5    270.7909      2017 NORTE DE SANTANDER  0.8281364 0.02220518
6    138.9222      2017 NORTE DE SANTANDER  0.5561253 0.01136823
7    615.8199      2017 NORTE DE SANTANDER  1.6560289 0.05045724
8    166.2543      2017 NORTE DE SANTANDER  0.7572865 0.01361349
9    374.9105      2017 NORTE DE SANTANDER  0.8858476 0.03069885
10   175.0128      2017 NORTE DE SANTANDER  0.6020091 0.01433714
                         geometry
1  MULTIPOLYGON (((-72.4778 8....
2  MULTIPOLYGON (((-73.01687 8...
3  MULTIPOLYGON (((-72.73134 7...
4  MULTIPOLYGON (((-72.60265 7...
5  MULTIPOLYGON (((-72.95019 8...
6  MULTIPOLYGON (((-72.62101 7...
7  MULTIPOLYGON (((-73.04222 7...
8  MULTIPOLYGON (((-72.58771 7...
9  MULTIPOLYGON (((-72.79776 7...
10 MULTIPOLYGON (((-72.63625 7...

precip.crop <- raster::crop(precip, extent(aoi))

precip.mask <- mask(x = precip.crop, mask = aoi)

precip.mask

class      : RasterLayer 
dimensions : 49, 32, 1568  (nrow, ncol, ncell)
resolution : 0.05, 0.05  (x, y)
extent     : -73.65, -72.05, 6.849999, 9.299999  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : chirps.v2.0.2020.04.6 
values     : 2.060457, 38.14587  (min, max)

Plotear:

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

leaflet() %>% addTiles() %>%
  addRasterImage(precip.mask, colors = pal, opacity = 0.6) %>%
  addLegend(pal = pal, values = values(precip.mask),
    title = "CHIRPS precipitaciones en Norte de Santander desde 26.04 hasta 30.04 de 2020 [mm]")

ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Inf

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

class       : SpatialPointsDataFrame 
features    : 721 
extent      : -73.575, -72.075, 6.924999, 9.274999  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       : chirps.v2.0.2020.04.6 
min values  :      2.06045699119568 
max values  :      38.1458702087402 

names(precip.points) <- "Lluvia"

precip.points

class       : SpatialPointsDataFrame 
features    : 721 
extent      : -73.575, -72.075, 6.924999, 9.274999  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       :           Lluvia 
min values  : 2.06045699119568 
max values  : 38.1458702087402 

str(precip.points)

Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 721 obs. of  1 variable:
  .. ..$ Lluvia: num [1:721] 21.53 13.99 16.87 16.75 8.82 ...
  ..@ coords.nrs : num(0) 
  ..@ coords     : num [1:721, 1:2] -73 -73.1 -73.1 -73 -73.2 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "x" "y"
  ..@ bbox       : num [1:2, 1:2] -73.57 6.92 -72.07 9.27
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

plot(precip.mask, main= "Datos CHIRPS de precipitaciones desde el 26 al 30.04.2020 [mm]")
plot(aoi, add=TRUE)

ignoring all but the first attribute

points(precip.points$x, precip.points$y, col = "darkblue", cex = 0.4)

aoi$area = st_area(aoi)

aoi

Simple feature collection with 40 features and 10 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -73.63379 ymin: 6.872201 xmax: -72.04761 ymax: 9.290847
geographic CRS: WGS 84
First 10 features:
   DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR                           MPIO_CRSLC
1          54      54001     CÚCUTA                                 1972
2          54      54003     ABREGO                                 1806
3          54      54051  ARBOLEDAS                                 1835
4          54      54099  BOCHALEMA                                 1826
5          54      54109 BUCARASICA                                 1838
6          54      54125     CÁCOTA                                 1630
7          54      54128    CÁCHIRA                                 1911
8          54      54172  CHINÁCOTA                                 1775
9          54      54223  CUCUTILLA                                 1834
10         54      54239    DURANIA Ordenanza 12 del 27 de Marzo de 1911
   MPIO_NAREA MPIO_NANO         DPTO_CNMBR Shape_Leng Shape_Area
1   1133.2031      2017 NORTE DE SANTANDER  2.8636253 0.09289252
2   1382.4498      2017 NORTE DE SANTANDER  2.2393759 0.11336240
3    456.1490      2017 NORTE DE SANTANDER  1.1611936 0.03736204
4    170.2669      2017 NORTE DE SANTANDER  0.7565329 0.01394455
5    270.7909      2017 NORTE DE SANTANDER  0.8281364 0.02220518
6    138.9222      2017 NORTE DE SANTANDER  0.5561253 0.01136823
7    615.8199      2017 NORTE DE SANTANDER  1.6560289 0.05045724
8    166.2543      2017 NORTE DE SANTANDER  0.7572865 0.01361349
9    374.9105      2017 NORTE DE SANTANDER  0.8858476 0.03069885
10   175.0128      2017 NORTE DE SANTANDER  0.6020091 0.01433714
                         geometry             area
1  MULTIPOLYGON (((-72.4778 8.... 1132277146 [m^2]
2  MULTIPOLYGON (((-73.01687 8... 1382095614 [m^2]
3  MULTIPOLYGON (((-72.73134 7...  455945007 [m^2]
4  MULTIPOLYGON (((-72.60265 7...  170163472 [m^2]
5  MULTIPOLYGON (((-72.95019 8...  270683414 [m^2]
6  MULTIPOLYGON (((-72.62101 7...  138836998 [m^2]
7  MULTIPOLYGON (((-73.04222 7...  615651563 [m^2]
8  MULTIPOLYGON (((-72.58771 7...  166142537 [m^2]
9  MULTIPOLYGON (((-72.79776 7...  374725056 [m^2]
10 MULTIPOLYGON (((-72.63625 7...  174908727 [m^2]

(border_sf <- aoi %>%
   summarise(area = sum(area)/1000000))

Simple feature collection with 1 feature and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -73.63379 ymin: 6.872201 xmax: -72.04761 ymax: 9.290847
geographic CRS: WGS 84
            area                       geometry
1 21846.28 [m^2] POLYGON ((-72.4556 7.553288...

(border <- as(border_sf, "Spatial"))

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : -73.63379, -72.04761, 6.872201, 9.290847  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       :             area 
value       : 21846.2845105745 

p.sf.magna2 <- st_transform(st_as_sf(precip.points), crs= 3116)

NSantander.sf.magna <- st_transform(aoi, crs= 3116)

precip2 <- as(p.sf.magna2, "Spatial")

precip2$precipitaciones <- round(precip2$Lluvia, 1)

precip2

class       : SpatialPointsDataFrame 
features    : 721 
extent      : 1055435, 1221276, 1257812, 1517605  (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 
variables   : 2
names       :           Lluvia, precipitaciones 
min values  : 2.06045699119568,             2.1 
max values  : 38.1458702087402,            38.1 

(NSantander2 <- as(NSantander.sf.magna, "Spatial"))

class       : SpatialPolygonsDataFrame 
features    : 40 
extent      : 1048947, 1224322, 1252020, 1519363  (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 
variables   : 10
names       : DPTO_CCDGO, MPIO_CCDGO,        MPIO_CNMBR,                          MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO,         DPTO_CNMBR,     Shape_Leng,       Shape_Area,             area 
min values  :         54,      54001,            ABREGO,                                1555,   44.94540102,      2017, NORTE DE SANTANDER, 0.455016349208, 0.00368621836307, 44907945.6872762 
max values  :         54,      54874, VILLA DEL ROSARIO, Ordenanza 9 de Noviembre 25 de 1948, 2680.07859017,      2017, NORTE DE SANTANDER,  4.08629755008,   0.220098910874,  2678750166.7455 

precip2@bbox <- NSantander2@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,]

ptrain <- spTransform(ptos_train, crs(precip.mask))
ptest <- spTransform(ptos_test, crs(precip.mask))

Plotear

lplot <- leaflet(data = precip2) %>% # data = original body - to get the zoom right
  addProviderTiles("CartoDB.Positron") %>% 
  addRasterImage(precip.mask, colors = pal, opacity = 0.6) %>%
  addCircleMarkers(data = ptrain, # first group
                   radius = 1,
                   fillOpacity = .7,
                   stroke = FALSE,
                   popup = ~htmlEscape(precipitaciones),
                   color = pal(ptos_train$precipitaciones), # using already created palette
                   clusterOptions = markerClusterOptions(),
                   group = "Training") %>% 
  addCircleMarkers(data = ptest, # second group
                   radius = 10,
                   fillOpacity = .7,
                   stroke = FALSE,
                   popup = ~htmlEscape(precipitaciones),
                   color = pal(ptos_test$precipitaciones), # using already created palette
                   clusterOptions = markerClusterOptions(),
                   group = "Test") %>% 
  addLegend(position = "bottomright",
            values = ~precipitaciones,
            opacity = .7,
            pal = pal, # palette declared previously
            title = "Precipitación") %>% 
  leaflet::addLayersControl(overlayGroups = c("Training", "Test"),
                   options = layersControlOptions(collapsed = FALSE)) %>% 
  addResetMapButton()

lplot

Interpolación

1. Polígonos Thiessen

th <- as(dirichlet(as.ppp(ptos_train)), "SpatialPolygons")

crs(th) <- crs(ptos_train)
crs(NSantander2) <- crs(ptos_train)

crs(th)

CRS arguments:
 +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 

crs(ptos_train)

CRS arguments:
 +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 

th.z <- over(th, ptos_train, fn=mean)

th.spdf <- SpatialPolygonsDataFrame(th, th.z)

th.clp <- raster::intersect(NSantander2, th.spdf)

tm_shape(th.clp) + 
  tm_polygons(col="precipitaciones", palette="RdBu", midpoint=30.0,
              title="Polígonos Thiessen\n Precipitación predicha\n [en mm])") +
  tm_legend(legend.outside=TRUE)

2. IDW

grd              <- as.data.frame(spsample(ptos_train, "regular", n=500000))
# You need to figure out what is the expected size of the output grd
names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE  # Create SpatialPixel object
fullgrid(grd)    <- TRUE  # Create SpatialGrid object

# Add P's projection information to the empty grid
proj4string(grd) <- proj4string(ptos_train)

# Interpolate the grid cells using a power value of 2 (idp=2.0)
P.idw <- gstat::idw(precipitaciones ~ 1, ptos_train, newdata=grd, idp=2.0)

[inverse distance weighted interpolation]

# Convert to raster object then clip to AOI
r       <- raster(P.idw)

r

class      : RasterLayer 
dimensions : 900, 556, 500400  (nrow, ncol, ncell)
resolution : 288.7443, 288.7443  (x, y)
extent     : 1060841, 1221383, 1257800, 1517670  (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     : 2.100869, 38.09815  (min, max)

NSantander2

class       : SpatialPolygonsDataFrame 
features    : 40 
extent      : 1048947, 1224322, 1252020, 1519363  (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 
variables   : 10
names       : DPTO_CCDGO, MPIO_CCDGO,        MPIO_CNMBR,                          MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO,         DPTO_CNMBR,     Shape_Leng,       Shape_Area,             area 
min values  :         54,      54001,            ABREGO,                                1555,   44.94540102,      2017, NORTE DE SANTANDER, 0.455016349208, 0.00368621836307, 44907945.6872762 
max values  :         54,      54874, VILLA DEL ROSARIO, Ordenanza 9 de Noviembre 25 de 1948, 2680.07859017,      2017, NORTE DE SANTANDER,  4.08629755008,   0.220098910874,  2678750166.7455 

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

r.m

class      : RasterLayer 
dimensions : 900, 556, 500400  (nrow, ncol, ncell)
resolution : 288.7443, 288.7443  (x, y)
extent     : 1060841, 1221383, 1257800, 1517670  (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     : 2.100869, 38.09815  (min, max)

Plotear

tm_shape(r.m) + 
  tm_raster(n=10,palette = "RdBu", auto.palette.mapping = FALSE,
            title="Distancia Inversa Ponderada\n Precipitación prevista [en mm]") + 
  tm_shape(ptos_train) + tm_dots(size=0.05, col ="black") +
  tm_legend(legend.outside=TRUE)

The argument auto.palette.mapping is deprecated. Please use midpoint for numeric data and stretch.palette for categorical data to control the palette mapping.

P <- ptos_train

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

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

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

  abline(0,1)

par(OP)

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

[1] 2.371934

Cross-validation

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

class       : SpatialPointsDataFrame 
features    : 100 
extent      : 1060831, 1204606, 1268823, 1512058  (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 
variables   : 2
names       :           Lluvia, precipitaciones 
min values  : 2.60495686531067,             2.6 
max values  : 38.1458702087402,            38.1 

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

# Remove a point then interpolate (do this n times for each point)
st <- stack()
for (i in 1:n){
  Z1 <- gstat::idw(precipitaciones~1, P[-i,], newdata=grd, idp=2.0)
  st <- addLayer(st,raster(Z1,layer=1))
  # Calculated pseudo-value Z at j
  Zi[,i] <- n * img$var1.pred - (n-1) * Z1$var1.pred
}

# Jackknife estimator of parameter Z at location j
Zj <- as.matrix(apply(Zi, 1, sum, na.rm=T) / n )

# Compute (Zi* - Zj)^2
c1 <- apply(Zi,2,'-',Zj)            # Compute the difference
c1 <- apply(c1^2, 1, sum, na.rm=T ) # Sum the square of the difference

# Compute the confidence interval
CI <- sqrt( 1/(n*(n-1)) * c1)

# Create (CI / interpolated value) raster
img.sig   <- img
img.sig$v <- CI /img$var1.pred 

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

# Plot the map
tm_shape(r.m) + tm_raster(n=7 ,title="IDW\n95% intervalo de confianza \n(en mm)") +
  tm_shape(P) + tm_dots(size=0.2) +
  
  tm_legend(legend.outside=TRUE)

Cambiar o revisar los intervalos donde muestran el error.

3. Kriging

f.1 <- as.formula(precipitaciones ~ X + Y) 
P$X <- coordinates(P)[,1]
P$Y <- coordinates(P)[,2]

# Compute the sample variogram; note that the f.1 trend model is one of the
# parameters passed to variogram(). This tells the function to create the 
# variogram on the de-trended data.
var.smpl <- variogram(f.1, P, cloud = FALSE, cutoff=100000, width=8990)

# Compute the variogram model by passing the nugget, sill and range values
# to fit.variogram() via the vgm() function.
dat.fit  <- fit.variogram(var.smpl, fit.ranges = TRUE, fit.sills = TRUE,
                          vgm(psill=3, model="Gau", range=150000, nugget=0.2))

dat.fit

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

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

[using universal kriging]

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

r.m

class      : RasterLayer 
dimensions : 900, 556, 500400  (nrow, ncol, ncell)
resolution : 288.7443, 288.7443  (x, y)
extent     : 1060841, 1221383, 1257800, 1517670  (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     : 1.414373, 36.52093  (min, max)

tm_shape(r.m) + 
  tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE, midpoint = NA,
            title="Universal Kriging\nPrecipitación predicha \n(en mm)") +
  tm_shape(P) + tm_dots(size=0.02, col="black") +
  tm_legend(legend.outside=TRUE)

The argument auto.palette.mapping is deprecated. Please use midpoint for numeric data and stretch.palette for categorical data to control the palette mapping.

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

leaflet() %>% addTiles() %>%
  addRasterImage(r.m, colors = pal, opacity = 0.6) %>%
  addLegend(pal = pal, values = values(r.m),
    title = "Precipitación interpolada con Kriging en Norte de Santander desde 26.04 al 30.04.2020 [mm]")

ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -InfSome values were outside the color scale and will be treated as NASome values were outside the color scale and will be treated as NA

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

tm_shape(r.m) + 
  tm_raster(n=7, palette ="Reds",
            title="Interpolación kriging\nMapa de varianza \n(en mm cuadrados)") +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, NSantander2)

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, col="black") +
  tm_legend(legend.outside=TRUE)