Introducción

En este cuaderno se exponen dos tecnicas de interpolación, estas son: IDW (Distancia Inversa Ponderada) y OK (kriging Ordinario). Ambas tecnicas se usan para tener una superficie continua del carbono organico del suelo (SOC) a una profundidad de 15 - 30 cm a partir de muestras obtenidas de SoilGrids 250m.

Configuración

Primero limpiremos la memeoria

rm(list=ls())

Cargar la librerias necesarias para el trabajo. Antes de cargarlas, verifica que estas hayan sido instaladas previamente.

library(sp)
library(terra)
library(sf)
library(stars)
library(gstat)
library(automap)
library(leaflet)
library(leafem)
library(ggplot2)
library(dplyr)
library(curl)
h <- new_handle()
handle_setopt(h, http_version = 2)

Leer los datos de entrada

Con el fin de imitar los datos del mundo real, vamos a leer la capa SOC que descargamos de ISRIC hacinedo uso de la libreria terra.

archivo <- ('C:/Users/LENOVO/Desktop/Geomatica/R/SOC15.tif')
(soc <- rast(archivo))
## class       : SpatRaster 
## dimensions  : 728, 816, 1  (nrow, ncol, nlyr)
## resolution  : 0.002261029, 0.00239011  (x, y)
## extent      : -76.4, -74.555, 4.06, 5.8  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : SOC15.tif 
## name        : SOC15 
## min value   :     0 
## max value   :  2921

Ahora cambiaremos las unidades de los datos, las pasaremos de dg/kg a g/kg utilizando su factor de conversion el cual es 10.

soc.perc <-  soc/10

Transformación de este CRS al CRS WGS84

geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(soc.perc, geog))
## class       : SpatRaster 
## dimensions  : 750, 796, 1  (nrow, ncol, nlyr)
## resolution  : 0.002319123, 0.002319123  (x, y)
## extent      : -76.4, -74.55398, 4.060658, 5.8  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        :    SOC15 
## min value   :   0.0000 
## max value   : 287.2864

Convirtamos la capa SpatRaster en un stars object

stars.soc = st_as_stars(geog.soc)

Ahora relizamos un mapa

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      stars.soc,
      opacity = 0.8,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 8:130)
    ) 
#
m

Muestreo del mundo real

A continuación obtendremos una muestra de aproximadamente 500 sitios del mundo real utilizando una muestra localizada aleatoriamente.

set.seed(123456)

# Muestreo aleatoreo de 500 puntos
(samples <- spatSample(geog.soc, 500, "random", as.points=TRUE))
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 500, 1  (geometries, attributes)
##  extent      : -76.39884, -74.55514, 4.061817, 5.79884  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : SOC15
##  type        : <num>
##  values      : 201.5
##                72.28
##                20.74

Debemos describir las principales caracteristicas del objeto muestra

Ahora, tenemos que convertir el objeto SpatVector en un objeto feature simple

(muestras <- sf::st_as_sf(samples))
nmuestras <- na.omit(muestras)

Visualicemos las muestras

longit <- st_coordinates(muestras)[,1]
latit <- st_coordinates(muestras)[,2]
soc <- muestras$SOC15
id <- seq(1,500,1) 
(sitios <- data.frame(id, longit, latit, soc))

Ahora borraremos los valores NA

sitios <- na.omit(sitios)
head(sitios)
##   id    longit    latit       soc
## 1  1 -76.32695 5.662012 201.52055
## 2  2 -76.19244 5.703756  72.27789
## 3  3 -74.62703 5.508950  20.73886
## 4  4 -75.63121 5.251527  47.95502
## 5  5 -76.18780 5.272399 124.67467
## 6  6 -75.13724 5.629544  76.56056

Visualizemos las muestras

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      stars.soc,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 8:130)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions())
m

Interpolación

Creación de objeto´gstat

Para interpolar, primero necesitamos crear un objeto de la clase gstat.

Un objeto gstat contiene toda la información necesaria para realizar la interpolación espacial.

Vamos a utilizar tres parámetros de la función gstat: - Fórmula: la “fórmula” de predicción que especifica las variables dependientes e independientes - Data: los datos de calibración - Modelo: el modelo de variograma

Interpolación IDW

Para interpolar utilizando este metodo creamos el siguinte objeto gstat, especificando unicamente la formula y los datos:

g1 = gstat(formula = SOC15  ~  1, data = nmuestras)

Ahora que el modelo de interpolación g1 está definido, podemos utilizar la función predict para interpolar.

Vamos a crear un objeto raster con valores de celda iguales a 1:

rrr = aggregate(geog.soc, 4)

¿Que es rrr?

rrr
## class       : SpatRaster 
## dimensions  : 188, 199, 1  (nrow, ncol, nlyr)
## resolution  : 0.009276493, 0.009276493  (x, y)
## extent      : -76.4, -74.55398, 4.056019, 5.8  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        :    SOC15 
## min value   :   0.0000 
## max value   : 255.0888

Definiremos nuevos valores para rrr

values(rrr) <-1

Definiremos nuevos nombres para rrr

names(rrr) <- "valor"

Ahora rrr es:

rrr
## class       : SpatRaster 
## dimensions  : 188, 199, 1  (nrow, ncol, nlyr)
## resolution  : 0.009276493, 0.009276493  (x, y)
## extent      : -76.4, -74.55398, 4.056019, 5.8  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : valor 
## min value   :     1 
## max value   :     1

Covertimos rrr en un objeto stars

stars.rrr = st_as_stars(rrr)

El siguiente codigo interpola los valores SOC según el modelo definido en g1 y el raster definido en stars.rrr:

z1 = predict(g1, stars.rrr)
## [inverse distance weighted interpolation]

¿Que es Z1?

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                 Min.  1st Qu.   Median     Mean  3rd Qu.     Max.  NA's
## var1.pred  0.9032976 34.64183 48.23325 51.72116 60.39388 242.9888     0
## var1.var          NA       NA       NA      NaN       NA       NA 37412
## dimension(s):
##   from  to offset       delta                       refsys x/y
## x    1 199  -76.4  0.00927649 +proj=longlat +datum=WGS8... [x]
## y    1 188    5.8 -0.00927649 +proj=longlat +datum=WGS8... [y]

En el objeto z1 hay dos atributos, solo el primero tiene información. A este lo renombraremos como “soc”:

z1 = z1["var1.pred",,]
names(z1) = "soc"

Necesitrampos una paleta de colores

paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), domain = 10:100, na.color = "transparent")

El raster SOC interpolado lo podemos observar a continuación:

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.6,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 11:55)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal=paleta, values= z1$soc,
    title = "Interpolación IDW de SOC [g/kg]"
    )
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
m  # Print the map

Interpolación OK

Primeramente calcularemos y examinaremos el variograma utilizando la función Variogram.

La función requiere dos argumentos: formula y datos

El siguiente codigo calcula el variograma del objeto “muestras”, sin covariables.

v_emp_ok = variogram(SOC15 ~ 1, data=nmuestras)

Ilustermos el variograma:

plot(v_emp_ok)

Utilizaremos la función autofitVariogram para ajustar un modelo de variograma a nuestro variograma.

v_mod_ok = autofitVariogram(SOC15 ~ 1, as(nmuestras, "Spatial"))

La función elige el tipo de modelo que mejor se ajusta y también afina sus parámetros.

El modelo ajustado puede representarse mediante un plot:

plot(v_mod_ok)

El objeto resultante posee una lista con varios componentes,el variograma y el modelo de variograma ajustado. El componente $var_model del objeto resultante contiene el modelo de variograma actual:

v_mod_ok$var_model
##   model      psill    range kappa
## 1   Nug   20.52881  0.00000   0.0
## 2   Ste 1171.51991 79.38399   0.3

Los componenetes resultantes a partir del ploteo del objeto v_mod_ok representan lo siguiente:

  • Model: El modelo de variograma ajustado.
  • Kappa o lag: Distancia entre puntos
  • Nugget: Valor del variograma con lag 0; errores en las mediciones.
  • sill: Valor en el que el semivariograma se aplana por primera vez. El nivel máximo de semivarianza.
  • Range: El punto en el que el semivariograma alcanza el sill en el eje x.

Ahora, el modelo puede pasarse a la función gstat para asi poder continuar con la interpolación OK:

g2 = gstat(formula = SOC15 ~ 1, model = v_mod_ok$var_model, data = nmuestras)
z2= predict(g2, stars.rrr)
## [using ordinary kriging]

Nuevamente vamos a cambiar el nombre del atributo var1.pred:

z2 = z2["var1.pred",,]
names(z2) = "soc"

Las prediciones de Kriging ordinario se puede ver a continuación.

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 11:55)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal = paleta, values= z2$soc,
    title = "Interpolación OK de SOC [g/kg]"
    )
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
m  # Print the map

Evaluación de los resultados

Evaluación cualitativa

Otra forma de visualizar los tres resultados de interpolación:

colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), domain = 10:100, na.color = "transparent")
m <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(stars.soc, opacity = 0.8, colorOptions = colores, group="Mundo real") %>%
  addGeoRaster(z1, colorOptions = colores, opacity = 0.8, group= "IDW")  %>%
  addGeoRaster(z2, colorOptions = colores, opacity = 0.8, group= "OK")  %>%
  # Add layers controls
  addLayersControl(
    overlayGroups = c("Mundo real", "IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
    addLegend("bottomright", pal = paleta, values= z1$soc,
    title = "Carbono organico del suelo [g/kg]"
)
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
m  # Print the map

Validación cruzada

Vamos a comparar los datos de los dos metodos de interpolación, esto lo haermos a partir de la validación Leave-One-Out, para esto: Eliminamos un punto de los datos de calibración, hacemos una predicción para dicho punto y obtenemos una tabla con un valor observado y un valor predicho para todos los puntos.

Validación cruzada Leave-One-Out:

cv1 = gstat.cv(g1)
cv2 = gstat.cv(g2)

los resultados de gstat.cv contiene los siguintes atributos: - var1.pred — Valor previsto - var1.var — Variance (solo para Kriging) - observed — valor observado - residual — Valor observado - previsto - zscore—Z - score (solo para Kriging) - fold—Cross - ID de validación

cv1 = na.omit(cv1)
cv1
## class       : SpatialPointsDataFrame 
## features    : 500 
## extent      : -76.39884, -74.55514, 4.061817, 5.79884  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 6
## names       :        var1.pred, var1.var,         observed,          residual, zscore, fold 
## min values  : 17.2576789621635,       NA,                0, -83.6916712396385,     NA,    1 
## max values  :  184.56223330443,       NA, 246.571563720703,  118.695562119791,     NA,  500

Convertiremos el objeto cv1

cv1 = st_as_sf(cv1)

Ahora trazemos los residuales

sp::bubble(as(cv1[, "residual"], "Spatial"))

Ahora calcularemos los indices de precisión de la predicción, como el error cuadratico medio (RMSE)

sqrt(sum((cv1$var1.pred - cv1$observed)^2) / nrow(cv1))
## [1] 21.26027

repetiremos el proceso para los resultados de OK

Conversión:

cv2 = st_as_sf(cv2)

calculo de RSME para OK

sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))
## [1] 18.94742

Con base en los resultados obtenidos a partir de la validación cruzada para los dos metodos de interpolación usados, podemos deducir que el metodo IDW tiene menor precisicion que el metodo OK, ya que el error cuadratico medio de este, 21.26027, es mas alto que el presentado en el Kriging ordinario, 18.94742.

Referencias

Este trabajo se realizó tomando como guía el cuaderno: Lizarazo, I. 2023, Spatial interpolation of soil organic carbon.

Cite este trabajo así: Cely, N., 2023. Interpolación espacial del carbono organico del suelo. Disponible en https://rpubs.com/ncely/1056095

sessionInfo()
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8  LC_CTYPE=Spanish_Colombia.utf8   
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.utf8    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] curl_5.0.0    dplyr_1.1.2   ggplot2_3.4.2 leafem_0.2.0  leaflet_2.1.2
##  [6] automap_1.1-9 gstat_2.1-1   stars_0.6-1   abind_1.4-5   sf_1.0-13    
## [11] terra_1.7-29  sp_1.6-0     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.10        lattice_0.20-45    FNN_1.1.3.2        png_0.1-8         
##  [5] class_7.3-20       zoo_1.8-12         digest_0.6.31      utf8_1.2.3        
##  [9] R6_2.5.1           plyr_1.8.8         evaluate_0.21      e1071_1.7-13      
## [13] highr_0.10         pillar_1.9.0       rlang_1.1.1        rstudioapi_0.14   
## [17] raster_3.6-20      jquerylib_0.1.4    rmarkdown_2.21     rgdal_1.6-6       
## [21] htmlwidgets_1.6.2  munsell_0.5.0      proxy_0.4-27       compiler_4.2.2    
## [25] xfun_0.39          pkgconfig_2.0.3    base64enc_0.1-3    htmltools_0.5.5   
## [29] tidyselect_1.2.0   tibble_3.2.1       intervals_0.15.3   codetools_0.2-18  
## [33] reshape_0.8.9      fansi_1.0.4        spacetime_1.3-0    withr_2.5.0       
## [37] grid_4.2.2         jsonlite_1.8.4     lwgeom_0.2-13      gtable_0.3.3      
## [41] lifecycle_1.0.3    DBI_1.1.3          magrittr_2.0.3     units_0.8-2       
## [45] scales_1.2.1       KernSmooth_2.23-20 cli_3.6.1          cachem_1.0.8      
## [49] farver_2.1.1       bslib_0.4.2        ellipsis_0.3.2     xts_0.13.1        
## [53] generics_0.1.3     vctrs_0.6.2        tools_4.2.2        glue_1.6.2        
## [57] crosstalk_1.2.0    parallel_4.2.2     fastmap_1.1.1      yaml_2.3.7        
## [61] colorspace_2.1-0   classInt_0.4-9     knitr_1.43         sass_0.4.6