Códigos para la interpolación

1.Introducción

Este cuaderno corresponde al anexo número 4 del informa final de la asignatura Geomática Básica.En el se llevarán a cabo tareas de interpolación para los datos de precipitación en el departamento de Vichada mediante el ajuste a un polinomio de primer orden y el metodo Kriging

2.Polinomio de Primer orden

Se selecciona el 70% de los datos de precipitación con la función sample , luego se llama “P” a este conjunto de datos

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,]
P <- ptos_train

Se plantea la ecuación de primer orden de la forma \(Precipitación= Intercepto+a{X}+b{X}\)

f.1 <- as.formula(precipitacion ~ X + Y) 

Se agregan X y Y al objeto P

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

Se realiza en modelo de regresión

lm.1 <- lm( f.1, data=P)

Se usa la salida del modelo de regresión para realizar la interpolación

dat.1st <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.1, newdata=grd))) 

Se recorta a la extensión del área de estudio

r   <- raster(dat.1st)
r.m1 <- raster::mask(r, Vichada2)

Se plotea el mapa con la interpolación obtenida

tm_shape(r.m1) + 
  tm_raster(n=10, palette="RdBu", midpoint = 43, 
            title="Interpolación con polinomio de 1er orden \n para precipitacion \n(mm)") +
  tm_shape(P) + tm_dots(size=0.02) +
  tm_legend(legend.outside=TRUE)

3.KRIGING

3.1 Variograma Experimental

Se usa un ajuste de primer orden para buscar eliminar una posible tendencia en los datos con el objetivo que la media y la variación de la precipitación sea constante.Al eliminar la tendencia deja unos residuos para utilizar en la interpolación .

Si se realiza la operación de elevar al cuadrado la diferencia de entre 2 valores de precipitación y luego dividirlos en 2 , se obtendrán unos puntos que al representarlos conforman el semivariograma

f.1 <- as.formula(precipitacion ~ X + Y) 
var.smpl <- variogram(f.1, P, cloud = FALSE, cutoff=200000, width=10000)

plot(var.smpl)

3.2Variograma de Ajuste a los Datos

Con ayuda de la libreria automap y la función autofitVariogram() se pueden obtener los parámetros de un posible modelo que se ajuste a los datos experimentales

Modelo <-  autofitVariogram(f.1, P)
summary(Modelo)

Experimental variogram:

Fitted variogram model:

Sums of squares betw. var. model and sample var.[1] 0.004027835

plot(Modelo, pch=16, col="black")

dat.fit  <- fit.variogram(var.smpl, fit.ranges = TRUE, fit.sills = TRUE, vgm(psill=68.711235, model="Gau", range=46602.27   , nugget=4))
plot(var.smpl, dat.fit, xlim=c(0,150000))

3.3 Superficie de Interpolación

Se define el modelo de tendencia

f.1 <- as.formula(precipitacion ~ X + Y) 

Se realiza la interpolación con la función krige()

dat.krg <- krige( f.1, P, grd, dat.fit)

Se crea un objeto ráster y se recorta a la extensión del área de interés

r1 <- raster(dat.krg)
r.k <- raster::mask(r1, Vichada2)
r.k

class      : RasterLayer 
dimensions : 313, 320, 100160  (nrow, ncol, ncell)
resolution : 1249.355, 1249.355  (x, y)
extent     : 1333231, 1733024, 798736.5, 1189785  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +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     : 0.2016994, 79.89707  (min, max)

Se plotea la salida de la interpolación

tm_shape(r.k) + 
  tm_raster(n=10, palette="RdBu", midpoint = NA,
            title="Universal Kriging\npara precipitacion \n(mm)") +
  tm_shape(P) + tm_dots(size=0.02, col="white") +
  tm_legend(legend.outside=TRUE)

CRS object has comment, which is lost in output

Se realiza un mapeo de la carianza y los intervalos de confianza del 95%

r2   <- raster(dat.krg, layer="var1.var")
r.j <- raster::mask(r2, Vichada2)

tm_shape(r.j) + 
  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)

CRS object has comment, which is lost in output

Los valores del mapa de intervalos de confianza representan las posibles subestimaciones o sobreestimaciones en los valores de precipitación

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

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

CRS object has comment, which is lost in output

sessionInfo()

R version 4.0.0 (2020-04-24)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.6 LTS

Matrix products: default
BLAS:   /usr/lib/atlas-base/atlas/libblas.so.3.0
LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8        LC_COLLATE=C.UTF-8    
 [5] LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8    LC_PAPER=C.UTF-8       LC_NAME=C             
 [9] LC_ADDRESS=C           LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

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

other attached packages:
[1] gstat_2.0-6    tmap_3.0       automap_1.0-14 raster_3.1-5   sp_1.4-2      

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.4.6       lattice_0.20-41    FNN_1.1.3          png_0.1-7          class_7.3-16      
 [6] zoo_1.8-8          digest_0.6.25      packrat_0.5.0      R6_2.4.1           plyr_1.8.6        
[11] evaluate_0.14      e1071_1.7-3        pillar_1.4.4       rlang_0.4.6        rstudioapi_0.11   
[16] rmarkdown_2.3      rgdal_1.5-10       stringr_1.4.0      htmlwidgets_1.5.1  munsell_0.5.0     
[21] compiler_4.0.0     xfun_0.15          pkgconfig_2.0.3    tmaptools_3.0      base64enc_0.1-3   
[26] htmltools_0.5.0    tidyselect_1.1.0   tibble_3.0.1       intervals_0.15.2   codetools_0.2-16  
[31] XML_3.99-0.3       reshape_0.8.8      spacetime_1.2-3    viridisLite_0.3.0  crayon_1.3.4      
[36] dplyr_1.0.0        sf_0.9-4           grid_4.0.0         jsonlite_1.6.1     lwgeom_0.2-5      
[41] lifecycle_0.2.0    DBI_1.1.0          magrittr_1.5       units_0.6-7        scales_1.1.1      
[46] KernSmooth_2.23-16 stringi_1.4.6      leafsync_0.1.0     leaflet_2.0.3      ellipsis_0.3.1    
[51] xts_0.12-0         generics_0.0.2     vctrs_0.3.1        RColorBrewer_1.1-2 tools_4.0.0       
[56] dichromat_2.0-0    leafem_0.1.1       glue_1.4.1         purrr_0.3.4        crosstalk_1.1.0.1 
[61] abind_1.4-5        parallel_4.0.0     yaml_2.2.1         colorspace_1.4-1   stars_0.4-1       
[66] classInt_0.4-3     knitr_1.28