1. Introducción

En este cuaderno veremos dos formas de interpolación espacial: Inverse Distance Weighted (IDW), y Ordinary Kriging (OK), utilizaremos estas técinas para obtener una interpretación e interpolación de los datos de Arcilla (clay) en el suelo a una profundidad de 15-30 cm a partir de las muestras obtenidas en SoilGrids250m.

2. Configuración

El primer paso es limpiar la memoria:

rm(list=ls())

Ahora, vamos a cargar las librerías necesarias para este cuaderno, previamente cargadas en la consola:

library(terra)
## terra 1.9.27
library(sf)
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
library(sp)
## Warning: package 'sp' was built under R version 4.6.1
library(stars)
## Warning: package 'stars' was built under R version 4.6.1
## Cargando paquete requerido: abind
library(gstat)
## Warning: package 'gstat' was built under R version 4.6.1
library(automap)
## Warning: package 'automap' was built under R version 4.6.1
library(RColorBrewer)
library(leaflet)
library(leafem)
## Warning: package 'leafem' was built under R version 4.6.1

Revisamos qué documentos podemos usar:

list.files()
##  [1] "clay_igh_15_30.tif"            "clay_stder.gpkg"              
##  [3] "IDW_clay_stder.tif"            "interpolacionespacial.html"   
##  [5] "interpolacionespacial.nb.html" "interpolacionespacial.Rmd"    
##  [7] "interpolacionespacial_files"   "muestreo.gpkg"                
##  [9] "municipios_col.gpkg"           "MunicipiosSantander.gpkg"     
## [11] "OK_clay_stder.tif"             "puntos.gpkg"                  
## [13] "rsconnect"                     "Soilgrids.nb.html"            
## [15] "Soilgrids.Rmd"

3. Lectura de datos

Leemos las muestras usando la librería sf:

samples <- sf::st_read("clay_stder.gpkg")
## Reading layer `clay_stder' from data source 
##   `C:\Users\isabe\OneDrive\Documentos\GB2\Proyecto7\clay_stder.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 2000 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.99655 ymin: 6.006384 xmax: -72.00211 ymax: 7.998851
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

También leemos la lista de municipios de nuestro departamento, Santander:

munic <- sf::st_read("MunicipiosSantander.gpkg")
## Reading layer `col_adm2' from data source 
##   `C:\Users\isabe\OneDrive\Documentos\GB2\Proyecto7\MunicipiosSantander.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 87 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.55179 ymin: 5.705701 xmax: -72.5161 ymax: 8.120555
## Geodetic CRS:  WGS 84

4. Exploración de los datos de entrada

Ahora vamos a realizar un histograma de nuestras muestras:

summary(samples)
##  clay_igh_15_30             geom     
##  Min.   :  0.0   POINT        :2000  
##  1st Qu.:343.5   epsg:NA      :   0  
##  Median :370.7   +proj=long...:   0  
##  Mean   :358.7                       
##  3rd Qu.:393.5                       
##  Max.   :488.0
hist(samples$clay)

Redondeamos los valores a dos dígitos:

samples$clay = round(samples$clay,2)

Para la visualización de las muestras, definiremos primero la paleta de colores a usar:

pal <- colorNumeric(c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"), domain = samples$clay)

Después ejecutamos:

leaflet() %>%
  addPolygons(
    data = munic,
    color = "gray",
    opacity = 1,
    weight = 1,
    fillOpacity = 0.2) %>%
 addCircleMarkers(
    data = samples,
    radius= 1.5, 
    label = ~clay,
    color = ~pal(clay),
    fillOpacity = 1,
    stroke = F
  ) %>%
  addLegend(
    data = samples,
    pal = pal,
    values = ~clay,
    position = "bottomleft",
    title = "CLAY:",
    opacity = 0.9) %>%
  addProviderTiles("OpenStreetMap")

5. Interpolación

Debemos crear un objeto gstat, para empezar a interpolar la arcilla (clay) con el método de IDW:

g1 = gstat(formula = clay ~1, data = samples)

Ahora usamos la librería Terra para crear objetos raster con valores de celda iguales a 1:

(rrr<- terra::rast(xmin = -74.55179, xmax = -72.00211, ymin = 5.705701, ymax = 8.120555, nrows = 370, ncols = 329, vals = 1, crs = "epsg:4326"))
## class       : SpatRaster
## size        : 370, 329, 1  (nrow, ncol, nlyr)
## resolution  : 0.007749787, 0.006526632  (x, y)
## extent      : -74.55179, -72.00211, 5.705701, 8.120555  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s)   : memory
## name        : lyr.1
## min value   :     1
## max value   :     1
stars.rrr <- stars::st_as_stars(rrr)
z1 = predict(g1, stars.rrr)
## [inverse distance weighted interpolation]
z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
## Warning in (function (..., deparse.level = 1) : number of columns of result is
## not a multiple of vector length (arg 1)
##                Min.  1st Qu.   Median    Mean  3rd Qu.     Max.    NAs     NA's
## var1.pred  1.334276 341.5757 365.7062 356.876 376.1191 484.6013      0 1.334276
## var1.var         NA       NA       NA     NaN       NA       NA 121730 0.000000
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -74.55   0.00775 WGS 84 [x]
## y    1 370  8.121 -0.006527 WGS 84 [y]

Para el ploteo creamos una paleta de colores:

library("terra")
library("leaflet")

Ahora ejecutamos para visualizar:

r_final <- terra::rast(z1[1])

pal_raster <- colorNumeric(
palette = "YlGn",
domain = terra::values(r_final),
na.color = "transparent"
)

leaflet() %>%
addProviderTiles("OpenStreetMap") %>%
addRasterImage(
r_final,
colors = pal_raster, opacity = 0.7 ) %>% addPolygons( data = munic, color = "gray", weight = 1, fillOpacity = 0.2 ) %>% addCircleMarkers( data = samples, radius = 1.5, label = ~clay, color = ~pal(clay), fillOpacity = 1, stroke = FALSE ) %>% addLegend( pal = pal_raster, values = terra::values(r_final), position = "bottomleft", title = "Soil Clay [%]" )
a <- which(is.na(z1[[1]]))
z1[[1]][a] = 0.0
z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
## Warning in (function (..., deparse.level = 1) : number of columns of result is
## not a multiple of vector length (arg 1)
##                Min.  1st Qu.   Median    Mean  3rd Qu.     Max.    NAs     NA's
## var1.pred  1.334276 341.5757 365.7062 356.876 376.1191 484.6013      0 1.334276
## var1.var         NA       NA       NA     NaN       NA       NA 121730 0.000000
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -74.55   0.00775 WGS 84 [x]
## y    1 370  8.121 -0.006527 WGS 84 [y]
names(z1) = "clay"

5.1 Interpolación con Krigin Ordinary

Ahora vamos a calcular y visualizar el variograma:

v_emp_ok = variogram(clay ~1, data= samples)
plot(v_emp_ok)

v_mod_ok = automap::autofitVariogram(clay ~1, as (samples, "Spatial"))
## Warning in automap::autofitVariogram(clay ~ 1, as(samples, "Spatial")): Some models where removed for being either NULL or having a negative sill/range/nugget, 
##  set verbose == TRUE for more information
plot(v_mod_ok)

v_mod_ok$var_model

El variograma que tenemos ahora podemos convertirlo en una función gstat:

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

Creamos un subconjunto de atributos de valores predichos y lo renombramos:

a<- which (is.na(z2[[1]]))
z2[[1]][a]=0.0
names(z2) = "clay"

Y en el siguiente ploteo vemos las predicciones:

paleta <- colorNumeric(
  palette = c("orange", "yellow", "cyan", "green"),
  domain = c(0, 100)
)
m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
    z2,
    opacity = 0.7,
    colorOptions = colorOptions(
      palette = c("orange", "yellow", "cyan", "green"),
      domain = range(z2$clay, na.rm = TRUE)
    )
  ) %>%
  addCircleMarkers(
    data = samples,
    radius = 1.5,
    label = ~clay,
    color = ~paleta(clay),
    fillOpacity = 1,
    stroke = FALSE
  ) %>%
  addLegend(
    "bottomright",
    pal = paleta,
    values = c(0,100),
    title = "OK CLAY interpolation [%]"
  )
## Warning in paleta(clay): Some values were outside the color scale and will be
## treated as NA
## Warning in paleta(clay): Some values were outside the color scale and will be
## treated as NA
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m

6. Evaluación de resultados

Esta es otra perspectiva de los resultados de interpolación:

colores <- colorOptions(
  palette = c("orange", "yellow", "cyan", "green"),
  domain = c(0,500),
  na.color = "transparent"
)
m <- leaflet() %>%
  addTiles()%>%
  addGeoRaster(z1, 
               colorOptions = colores,
               opacity = 0.8,
               group = "IDW"
               )%>%
  addGeoRaster(z2,
               colorOptions = colores,
               opacity = 0.8,
               group = "OK"
               )%>%
  addLayersControl(
    overlayGroups = c("IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  addLegend(
    "bottomright",
    pal = paleta,
    values =c(0,100),
    title = "Soil Clay [%]"
  )
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m

6.1 Validación cruzada

Una vez tenemos los gráficos de interpolación del contenido de arcilla en el suelo y su representación en mapa, podemos comparar los resultados, esto se realiza con la función gsta.cv y se realiza haciendo una predicción punto por punto excluyendo a cada uno de los datos de calibración, por ello es conveniente ejecutar dicha función desde la consola:

## sp::bubble(cv2, "residual")
## sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))
## [1] 45.533

Aplicamos los comandos anteriores ahora con la interpolación IDW, que en este caso será cv1:

## sp::bubble(cv1, "residual")
## sqrt(sum((cv1$var1.pred - cv1$observed)^2) / nrow(cv1))
## [1] 47.20467

Una vez obtenido el RMSE (error cuadrático medio), el cual indica lo lejos que están las predicciones de los valores reales, podemos objetivamente concluir que el método de Ordinary Kriging fue el que más precisión en la predicción de los valores.

7 Guardando los resultados de la interpolación

write_stars(
  z1, dsn = "./IDW_clay_stder.tif", layer = 1
)
write_stars(
  z2, dsn = "./OK_clay_stder.tif", layer = 1
)

8 Conclusiones

En este cuaderno realizamos la lectura de los datos descargados desde SoilGrids respecto al contenido de arcilla (clay) en los suelos de Santander, seleccionamos 2000 puntos de muestra, luego obtuvimos su histograma e hicimos la interpolación y su ploteo. Con la validación cruzada y el RMSE pudimos llegar a la conlcusión de que en este caso el método OK fue más preciso que el método IDW.

Bibliografía

Cite this work as follows: Lizarazo, I. 2025. Spatial interpolation. Available at: https://rpubs.com/ials2un/Spatial_Interpolation

sessionInfo()
## R version 4.6.0 (2026-04-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## 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    
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] leafem_0.2.5       leaflet_2.2.3      RColorBrewer_1.1-3 automap_1.1-20    
##  [5] gstat_2.1-6        stars_0.7-3        abind_1.4-8        sp_2.2-3          
##  [9] sf_1.1-1           terra_1.9-27      
## 
## loaded via a namespace (and not attached):
##  [1] generics_0.1.4          sass_0.4.10             class_7.3-23           
##  [4] KernSmooth_2.23-26      lattice_0.22-9          magrittr_2.0.5         
##  [7] digest_0.6.39           evaluate_1.0.5          grid_4.6.0             
## [10] fastmap_1.2.0           plyr_1.8.9              jsonlite_2.0.0         
## [13] e1071_1.7-17            reshape_0.8.10          DBI_1.3.0              
## [16] crosstalk_1.2.2         scales_1.4.0            codetools_0.2-20       
## [19] jquerylib_0.1.4         cli_3.6.6               rlang_1.2.0            
## [22] units_1.0-1             intervals_0.15.5        base64enc_0.1-6        
## [25] cachem_1.1.0            yaml_2.3.12             otel_0.2.0             
## [28] FNN_1.1.4.1             raster_3.6-32           tools_4.6.0            
## [31] parallel_4.6.0          dplyr_1.2.1             ggplot2_4.0.3          
## [34] spacetime_1.3-3         png_0.1-9               vctrs_0.7.3            
## [37] R6_2.6.1                zoo_1.8-15              proxy_0.4-29           
## [40] lifecycle_1.0.5         classInt_0.4-11         leaflet.providers_3.0.0
## [43] htmlwidgets_1.6.4       pkgconfig_2.0.3         pillar_1.11.1          
## [46] bslib_0.11.0            gtable_0.3.6            glue_1.8.1             
## [49] Rcpp_1.1.1-1.1          tidyselect_1.2.1        tibble_3.3.1           
## [52] xfun_0.57               rstudioapi_0.18.0       knitr_1.51             
## [55] farver_2.1.2            htmltools_0.5.9         rmarkdown_2.31         
## [58] xts_0.14.2              compiler_4.6.0          S7_0.2.2