Anexo 4.Interpolación por método Kriging de datos de precipitación en Boyacá
Marly Daniela Quiroga Buitrago
24/06/2020
En este cuaderno se muestra uno de los métodos de interpolación espacial aplicado a un conjunto de datos de precipitación correspondientes al departamento de Boyacá.
Kriging
Ajustar el modelo de variograma
Primero, necesitamos crear un modelo de variograma. Tenga en cuenta que el modelo de variograma se calcula sobre los datos de tendencia. Esto se implementa en el siguiente fragmento de código pasando el modelo de tendencia de primer orden (definido en un fragmento de código anterior como objeto de fórmula f.1) a la función variograma.
Es importante subir el mapa de puntos de precipitacion y el archivo que represente nuestra área de interés:
library(raster)## Loading required package: splibrary(gstat)Boyaca<-shapefile("C:/Users/Brian/Desktop/Daniela/geomatica/Informe final/boyaca2.shp")Precipitacion <- shapefile('C:/Users/Brian/Desktop/Daniela/geomatica/Informe final/precipitacion2.shp')Precipitacion@bbox <-Boyaca@bboxAjustar el modelo de variograma `
# Definir la ecuación polinómica de primer orden.
f.1 <- as.formula(lluvia ~ X + Y)
# Calcular el variograma de muestra; tenga en cuenta que el modelo de tendencia f.1 es uno de los parámetros pasados a variogram (). Esto le dice a la función que cree el variograma en los datos de tendencia.
var.smpl <- variogram(f.1, P, cloud = FALSE, cutoff=100000, width=8990)
# Calcule el modelo de variograma pasando los valores de nugget, sill y range a fit.variogram () a través de la función vgm ().
dat.fit <- fit.variogram(var.smpl, fit.ranges = TRUE, fit.sills = TRUE,
vgm(psill=20, model="Gau", range=150000, nugget=0.0))
## el modelo puede ser Exp, Gau, Sph, Mat Vea el paquete gstat - función vgm - para obtener los nombres de estos modelos Revise la teoría de kriging para comprender su significado.dat.fit## model psill range
## 1 Nug 4.516668 0.00
## 2 Gau 281.040415 43045.33# El siguiente diagrama nos permite evaluar el ajuste
plot(var.smpl, dat.fit, xlim=c(0,130000))
Generar superficie Kriged
Luego, use el modelo de variograma dat.fit para generar una superficie interpolada kriged. La función krige nos permite incluir el modelo de tendencia, lo que nos evita tener que reducir la tendencia de los datos, corregir los residuos y luego combinar los dos rásteres. En cambio, todo lo que tenemos que hacer es pasar krige la fórmula de tendencia f.1.
# Definir el modelo de tendencia
f.1 <- as.formula(lluvia ~ X + Y)
# Realice la interpolación de krige (tenga en cuenta el uso del modelo de variograma creado en el paso anterior)
dat.krg <- krige( f.1, P, grd, dat.fit)## [using universal kriging]# Convertir superficie kriged en un objeto ráster para el recorte
r <- raster(dat.krg)
r.m <- raster::mask(r, Boyaca)r.m## class : RasterLayer
## dimensions : 531, 600, 318600 (nrow, ncol, ncell)
## resolution : 500.2415, 500.2415 (x, y)
## extent : 934866.4, 1235011, 1006726, 1272354 (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 +units=m +no_defs
## source : memory
## names : var1.pred
## values : 3.369428, 97.49588 (min, max)# Trazar el mapa
library(tmap)
tm_shape(r.m) +
tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE,midpoint = NA,
title="Kriging universal \nPrecipitación prevista \n (en mm)") +
tm_shape(P) + tm_dots(size=0.2) +
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.
## class : SpatialPolygonsDataFrame
## features : 123
## extent : -74.66496, -71.94885, 4.655196, 7.055557 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 9
## names : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR, MPIO_CRSLC, MPIO_NAREA, MPIO_NANO, DPTO_CNMBR, Shape_Leng, Shape_Area
## min values : 15, 15001, ALMEIDA, 1537, 25.35271459, 2017, BOYACÃ, 0.204497978002, 0.00206924119091
## max values : 15, 15897, ZETAQUIRÃ, Ordenanza 8 de Diciembre 4 de 1965, 1513.60104233, 2017, BOYACÃ, 2.40391520946, 0.123609862522Una mejor visualización de la superficie interpolada:
library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("red", "orange", "yellow", "blue", "darkblue"), values(precipitacion.mask),
na.color = "transparent")
leaflet() %>% addTiles() %>%
addRasterImage(r.m, colors = pal, opacity = 0.6) %>%
addLegend(pal = pal, values = values(r.m),
title = "Interpolación Kriging de la precipitación en Boyacá\n del 26.04 a 30.04 in 2020 [mm]")## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NAGenerar los mapas de varianza e intervalo de confianza
El objeto dat.krg almacena no solo los valores interpolados, sino también los valores de varianza. Estos se pueden pasar al objeto ráster para la asignación de la siguiente manera:
r <- raster(dat.krg, layer="var1.var")
r.m <- raster::mask(r, Boyaca)
tm_shape(r.m) +
tm_raster(n=7, palette ="Reds",
title="Interpolación Kriging\nMapa de varianza \n(en mm cuadrado)") +tm_shape(P) + tm_dots(size=0.2) +
tm_legend(legend.outside=TRUE)
Un mapa más fácil de interpretar es el mapa del intervalo de confianza del 95% que se puede generar a partir del objeto de varianza de la siguiente manera (los valores del mapa deben interpretarse como el número de mm por encima y por debajo de la cantidad de lluvia estimada).
r <- sqrt(raster(dat.krg, layer="var1.var")) * 1.96
r.m <- raster::mask(r, Boyaca)
tm_shape(r.m) +
tm_raster(n=7, palette ="Reds",
title="Interpolación Kriging \n95% Mapa de CI \n(en mm)") +tm_shape(P) + tm_dots(size=0.2) +
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 8.1 x64 (build 9600)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Spanish_Colombia.1252 LC_CTYPE=Spanish_Colombia.1252
## [3] LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Colombia.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] RColorBrewer_1.1-2 leaflet_2.0.3 tmap_3.0 gstat_2.0-6
## [5] raster_3.0-12 sp_1.4-1
##
## loaded via a namespace (and not attached):
## [1] zoo_1.8-7 xfun_0.12 sf_0.9-3 lattice_0.20-38
## [5] colorspace_1.4-1 htmltools_0.4.0 stars_0.4-1 viridisLite_0.3.0
## [9] spacetime_1.2-3 yaml_2.2.1 base64enc_0.1-3 XML_3.99-0.3
## [13] rlang_0.4.5 e1071_1.7-3 DBI_1.1.0 lifecycle_0.2.0
## [17] stringr_1.4.0 munsell_0.5.0 htmlwidgets_1.5.1 leafsync_0.1.0
## [21] codetools_0.2-16 evaluate_0.14 knitr_1.28 crosstalk_1.1.0.1
## [25] parallel_3.6.3 class_7.3-15 leafem_0.1.1 xts_0.12-0
## [29] Rcpp_1.0.3 KernSmooth_2.23-16 scales_1.1.0 classInt_0.4-3
## [33] lwgeom_0.2-1 jsonlite_1.6.1 abind_1.4-5 FNN_1.1.3
## [37] farver_2.0.3 png_0.1-7 digest_0.6.25 stringi_1.4.6
## [41] tmaptools_3.0 grid_3.6.3 rgdal_1.4-8 tools_3.6.3
## [45] magrittr_1.5 dichromat_2.0-0 rmarkdown_2.1 R6_2.4.1
## [49] units_0.6-5 intervals_0.15.2 compiler_3.6.3