1. Introducción

En este notebook se aplican técnicas de interpolación espacial, dentro de las que destacan IDW y Kriging Ordinary (OK), estas permiten estimar diferentes variables que permiten evaluar a un terreno, en este caso, se estimará el contenido de carbono orgánico en el departamento del Chocó con una profundidad de 15 a 30 cm usando IDW que asigna valores según la proximidad de las muestras y OK que se encarga de incorporar la estructura espacial de los datos para mayor precisión. La base de datos usada para el desarrollo de este proyecto pertenece a SOIL GRIDS.

2. Configuración del notebook

Limpieza de la memoria del programa usando:

rm(list=ls())

Configuración de las librerías a usar en el cuaderno:

library(sp)
library(terra)
## terra 1.8.54
library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(stars)
## Cargando paquete requerido: abind
library(gstat)
library(automap)
library(leaflet)
library(leafem)
library(ggplot2)
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:terra':
## 
##     intersect, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(curl)
## Using libcurl 8.14.1 with Schannel

Haciendo uso del paquete curl, este código permite realizar solicitudes HTTP de manera eficiente y flexible, usando handle 2 es posible hacer que las solicitudes a la web sean rápidas y eficientes, evitando sobrecargar los servidores.

h <- new_handle()
handle_setopt(h, http_version = 2)

3. Lectura de los datos

El archivo anteriormente descargado de SOIL GRIDS, corresponde a los valores de SOC de 15-30 cm de profundidad, poseen un tamaño de celda de 250m*250m; estas se cuardan en archivo geopackge, y se convierte a archivo tipo raster.

archivo <- ("soc_igh_15_30.tif")
(soc <-  rast(archivo))
## class       : SpatRaster 
## size        : 2100, 866, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8689500, -8473000, 441250, 966250  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source      : soc_igh_15_30.tif 
## name        : soc_igh_15_30

Los datos anteriormente nombrados, vienen en unidades distintas a las aceptadas por el sistema internacional de medidas, así, debe hacerse la siguiente conversión, dividimos en 10 para conservar las cifras decimales significativas:

soc.perc <-  soc/10

Transformación de los datos de CRS a el sistema WGS84 CRS:

geog = "+proj=longlat + datum= WGS84"
(geog.soc = project(soc.perc, geog))
## class       : SpatRaster 
## size        : 2130, 980, 1  (nrow, ncol, nlyr)
## resolution  : 0.002214188, 0.002214188  (x, y)
## extent      : -78.0065, -75.83659, 3.963751, 8.679971  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_igh_15_30 
## min value   :      9.175247 
## max value   :    344.548981

Convertir la capa SpatRaster a un objeto stars

stars.soc = st_as_stars(geog.soc)

Presentación de los datos visualmente mediante un gráfico:

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

4. Muestreo de algunos datos del departamento

Simulación pseudoaleatoria de toma de muestras de algunos puntos dentro del departamento.

set.seed(123456)

#Selección de 500 puntos aleatorios dentro del departamento.

(samples <-  spatSample(geog.soc, 500, "random", as.points=TRUE))
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 500, 1  (geometries, attributes)
##  extent      : -78.00539, -75.85541, 3.969287, 8.672222  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : soc_igh_15_30
##  type        :         <num>
##  values      :         23.87
##                        197.4
##                           NA

SpatVector a simple feature

(muestras <- sf::st_as_sf(samples))

Se omiten celdas sin valores y luego se visualizan los datos limpios:

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

Remover los valores NA

sitios <- na.omit(sitios)
head(sitios)

Visualizar 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 de los datos

5. Creación de objetos gstat

Para poder realizar una interpolación, es necesario crear un objeto de tipo gstat, esto es posible haciendo uso de la librería con el mismo nombre. Un objeto gstat contiene la información necesaria para realizar

6. Interpolación IDW

Determina los valores de celda mediante la combinación ponderada linealmente de un conjunto de puntos de una muestra. La ponderación es una función de la distancia inversa. La superficie interpolada debe ser la de una variable dependiente de la localización. Este método se basa en que la influencia de la variable que se está cartografiando disminuye con la distancia desde su ubicación en la muestra.

Para estimar el contenido de carbono orgánico del suelo (SOC) en el departamento del Chocó, se emplea la interpolación espacial mediante el método de Inverse Distance Weighting IDW. Para ello, se utiliza la función gstat, que permite definir el modelo de interpolación con base en los datos de calibración obtenidos en diferentes puntos de muestreo. El modelo se construye de la siguiente manera:

g1 = gstat(formula = soc_igh_15_30 ~ 1, data = nmuestras)
# a simple copy
rrr = aggregate(geog.soc, 4)
rrr
## class       : SpatRaster 
## size        : 533, 245, 1  (nrow, ncol, nlyr)
## resolution  : 0.008856751, 0.008856751  (x, y)
## extent      : -78.0065, -75.83659, 3.959323, 8.679971  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_igh_15_30 
## min value   :      10.57474 
## max value   :     332.98460

definir los nuevos valores

values(rrr) <-1
names(rrr) <- "valor"
rrr
## class       : SpatRaster 
## size        : 533, 245, 1  (nrow, ncol, nlyr)
## resolution  : 0.008856751, 0.008856751  (x, y)
## extent      : -78.0065, -75.83659, 3.959323, 8.679971  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : valor 
## min value   :     1 
## max value   :     1
stars.rrr = st_as_stars(rrr)
## [inverse distance weighted interpolation]
z1 = predict(g1, stars.rrr)
## [inverse distance weighted interpolation]
z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean  3rd Qu.     Max.   NA's
## var1.pred  13.06132 59.54735 102.3513 111.0713 159.2068 294.8568      0
## var1.var         NA       NA       NA      NaN       NA       NA 130585
## dimension(s):
##   from  to offset     delta                       refsys x/y
## x    1 245 -78.01  0.008857 +proj=longlat +datum=WGS8... [x]
## y    1 533   8.68 -0.008857 +proj=longlat +datum=WGS8... [y]

`

z1 = z1["var1.pred",,]
z1
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##                Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
## var1.pred  13.06132 59.54735 102.3513 111.0713 159.2068 294.8568
## dimension(s):
##   from  to offset     delta                       refsys x/y
## x    1 245 -78.01  0.008857 +proj=longlat +datum=WGS8... [x]
## y    1 533   8.68 -0.008857 +proj=longlat +datum=WGS8... [y]
names(z1) = "soc"
paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), domain = 13:300, na.color = "transparent")
m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 13:300)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal=paleta, values= z1$soc,
    title = "IDW SOC interpolation [%]"
    )
# Imprimir el gráfico 
m

7. Interpolación OK

Se define como una técnica geoestadística avanzada que se utilisa para estimar balores en ubicaciones no muestradas basándose en los valores conocidos de lugares cercanos. Es decir, se basa en la premisa de que los puntos conocidos que son cercanos a cierto espacio, son los que mejor se asemejan al valor real. Kriging Ordinary no solo considera la distancia, ya que también estudia el patrón de variabilidad espacial de la variable de interés, que en este caso es el carbono orgánico en el departamento del Chocó, Colombia.

8. Seriovariograma

Este modelo permite la cuantificación de la autocorrelación espacial de los datos, asignando los pesos adecuados a las predicciones. Usando la función variogram, especificamos la variable dependiente y los datos de calibración:

v_emp_ok = variogram(soc_igh_15_30 ~ 1, data=nmuestras)
plot(v_emp_ok)

Luego, se busca un modelo matemático que se ajuste de la mejor manera al semivariograma. Usando la funión autofitVariogram del paquete automap, se identifica con precisión del modelo que mejor se ajusta a los parámetros obtenidos anteriormente.

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

v_mod_ok$var_model

Luego de obtener el modelo ajustado al semivariograma, se usa la función gstat para ejecutar la interpolación mediante Kriging Ordinary. En este método se estiman los valores de contenido de carbono orgánico en las ubicaciones del terreno que no se muestraron.

## [using ordinary kriging]
g2 = gstat(formula = soc_igh_15_30 ~ 1, model = v_mod_ok$var_model, data = nmuestras)
z2= predict(g2, stars.rrr)
## [using ordinary kriging]
z2 = z2["var1.pred",,]
names(z2) = "soc"

Para generar el mapa de interpolaci+on OK con el contenido de carbono orgánico en el suelo a una profundidad de 15 a 30 cm, se plantea el siguiente código:

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 13:300)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal = paleta, values= z2$soc,
    title = "OK SOC interpolation [%]"
    )
#Imprimir el gráfico
m

9. Evaluación de resultados

9.1 Análisis cualitativo

Para analizar las predicciones de carbono orgánico en el suelo en el Chocó usando los métodos de interpolación IDW (Inverse Distance Weighting) y OK (Kriging Ordinary) se usa un gráfico interactivo usando la función leaflet con las siguientes convensiones: Naranja; para valores bajos de SOC. Amarillo - Cyan; identifican valores intermedios de SOC. Verde; para altos contenidos de SOC en el territorio. En el gráfico interactivo puede compararse la base de datos descargada en SOIL GRIDS, el método IDW y OK, permitiendo hacer una descripción visual de las diferencias encontradas en cada una de las interpolaciones:

colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), domain = 13:300, na.color = "transparent")
m <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(stars.soc, opacity = 0.8, colorOptions = colores, group="RealWorld") %>%
  addGeoRaster(z1, colorOptions = colores, opacity = 0.8, group= "IDW")  %>%
  addGeoRaster(z2, colorOptions = colores, opacity = 0.8, group= "OK")  %>%
  # Add layers controls
  addLayersControl(
    overlayGroups = c("RealWorld", "IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
    addLegend("bottomright", pal = paleta, values= z1$soc,
    title = "Soil organic carbon [%]"
)
#Imprimir el gráfico
m

Al observar el gráfico anterior, es posible afirmar que la interpolación Kriging Ordinary posee un mayor acercamiento a los datos reales encontrados en SOIL GRIDS, afirmando que la consigna de darle mayor peso a los datos más cercanos tiene un mejor ajuste dentro de las interpolaciones planteadas.

9.2 Validación cruzada

Para evaluar la precisión de los métodos de interpolación usados en este notebook, se utiliza la Validación Cruzada Leave-One-Out (LOOCV). Consiste en excluir un punto del conjunto de datos de calibración para luego predecir el valor de ese punto usando los datos restantes y repetir el procedimiento para los demás datos.

Además, se calcula el Error Cuadrático Medio RMSE para evaluar la precisión de las predicciones usando la siguiente fórmula:

Valores cercanos a cero en RMSE indican una alta precisión en sus predicciones, mientras que valores más alejados de este sugieren mayores discrepancias entre los valores observados y los predichos por los modelos.

### Este código se corre desde la consola, linea por línea.
#cv1 = gstat.cv(g1)
#cv2 = gstat.cv(g2)
#cv1 = na.omit(cv1)
#cv1

Conversión de cv1 a un objeto de tipo simple features.

#cv1 = st_as_sf(cv1)

Para tener una representación efectiva se utiliza un gráfico de burbujas, en este, por ejemplo, observamos que las burbujas asemejan el área del departamento del Chocó,el tamaño de cada burbuja indica que tan lejos se encuentra de los valores reales, es decir, a mayor tamaño, mayor error en el cálculo, y de la misma forma, burbujas pequeñas refieren a errores pequeños respecto al cálculo de los datos con los métodos de interpolación. Se observa que los errores más grandes se encuentran en el centro del departamento, con valores positivos hasta 182,176 y negativos hasta -136,789. Así, se explica porqué al calcular el RMSE de OK y IDW, los valores llegan hasta 30, es decir, las predicciones realizadas difieren bastante de los datos reales dentro del departamento sobre el contenido de Carbono Orgánico.

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

# RMSE para IDW interpolation
#sqrt(sum((cv1$var1.pred - cv1$observed)^2) / nrow(cv1))

[1] 36.9009

#cv2 = st_as_sf(cv2)
# RMSE para OK interpolation
#sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))

[1] 30.77824

10.Conclusiones

Los métodos de interpolación como IDW y OK son herramientas esenciales, ya que permiten inferir matemáticamente, los posibles resultados de una propiedad en específico en un terreno, pero, aunque es un método muy efectivo, no siempre posee un acercamiento al valor real. De esta manera, su uso es de gran importancia en la agronomía para la toma decisiones informadas, acompañado con los estudios pertinentes y las metodologías adecuadas, se convierte en una herramienta fundamental, con un uso amplio en diferentes campos de las ciencias agrarias.

11. Referencias

sessionInfo()
## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## 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] curl_6.4.0     dplyr_1.1.4    ggplot2_3.5.2  leafem_0.2.4   leaflet_2.2.2 
##  [6] automap_1.1-16 gstat_2.1-3    stars_0.6-8    abind_1.4-8    sf_1.0-21     
## [11] terra_1.8-54   sp_2.2-0      
## 
## loaded via a namespace (and not attached):
##  [1] generics_0.1.4     sass_0.4.10        class_7.3-23       KernSmooth_2.23-26
##  [5] lattice_0.22-6     digest_0.6.37      magrittr_2.0.3     evaluate_1.0.3    
##  [9] grid_4.5.0         RColorBrewer_1.1-3 fastmap_1.2.0      plyr_1.8.9        
## [13] jsonlite_2.0.0     e1071_1.7-16       reshape_0.8.10     DBI_1.2.3         
## [17] crosstalk_1.2.1    scales_1.4.0       codetools_0.2-20   jquerylib_0.1.4   
## [21] cli_3.6.5          rlang_1.1.6        units_0.8-7        intervals_0.15.5  
## [25] withr_3.0.2        base64enc_0.1-3    cachem_1.1.0       yaml_2.3.10       
## [29] FNN_1.1.4.1        raster_3.6-32      tools_4.5.0        parallel_4.5.0    
## [33] spacetime_1.3-3    png_0.1-8          vctrs_0.6.5        R6_2.6.1          
## [37] zoo_1.8-14         proxy_0.4-27       lifecycle_1.0.4    classInt_0.4-11   
## [41] htmlwidgets_1.6.4  pkgconfig_2.0.3    pillar_1.10.2      bslib_0.9.0       
## [45] gtable_0.3.6       glue_1.8.0         Rcpp_1.0.14        tidyselect_1.2.1  
## [49] tibble_3.2.1       xfun_0.52          rstudioapi_0.17.1  knitr_1.50        
## [53] farver_2.1.2       htmltools_0.5.8.1  rmarkdown_2.29     xts_0.14.1        
## [57] compiler_4.5.0