1. Introducción

Este cuaderno ilustra dos técnicas de interpolación espacial: distancia ponderada inversa (IDW) y Kriging ordinario (OK). IDW es una técnica determinista. OK es probabilístico. Ambas técnicas se utilizan aquí para obtener una superficie continua de SOC a 15-30 cm a partir de muestras obtenidas de SoilGrids 250 m.

2. Setup

Primero debemos limpiar la memoria

rm(list=ls())

Debemos asegurarnos de tener instaladas las librerias necesarias. Luego debemos cargar las librerias

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

3. Leer los datos de entrada

Necesitamos leer un conjunto de datos para imitar datos del mundo real. Por lo tanto, leamos la capa SOC que descargamos de (ISRIC) usando la biblioteca terra:

archivo <- ("soc_igh_15_30.tif")
(soc <- rast(archivo))
## class       : SpatRaster
## size        : 759, 696, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8495500, -8321500, 453750, 643500  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine
## source      : soc_igh_15_30.tif
## name        : soc_igh_15_30

Ahora, vamos a convertir los datos del SOC en porcentaje.

soc.perc <-  soc/10

¿Cuál es el sistema de referencia de coordenadas (CRS) de los datos del mundo real? Necesitamos una transformación de dicho CRS al conocido CRS WGS84:

geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(soc.perc, geog))
## class       : SpatRaster
## size        : 772, 739, 1  (nrow, ncol, nlyr)
## resolution  : 0.002208737, 0.002208737  (x, y)
## extent      : -76.25632, -74.62406, 4.075514, 5.780659  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## source(s)   : memory
## name        : soc_igh_15_30
## min value   :      11.30435
## max value   :    257.617828

Vamos a convertir la capa SpatRaster en un objeto de estrellas:

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

4. Muestreo del mundo

Consigamos una muestra de aprox. 500 sitios a partir de datos del mundo real utilizando una muestra ubicada aleatoriamente:

set.seed(123456)

# Muestreo aleatorio de 500 puntos.
(samples <- spatSample(geog.soc, 500, "random", as.points=TRUE))
## class       : SpatVector
## geometry    : points
## dimensions  : 500, 1  (geometries, attributes)
## extent      : -76.25521, -74.62516, 4.085453, 5.777346  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names       : soc_igh_15_30
## type        :         <num>
## values      :            NA
##                     83.2353
##                     78.1224
##               ...

El objeto samples es un SpatVector que contiene una muestra aleatoria de 500 puntos extraídos del raster geog.soc. Cada punto representa una ubicación geográfica y almacena un valor numérico correspondiente al carbono orgánico del suelo (SOC). El objeto está en formato vectorial de puntos, conserva el sistema de referencia WGS84 y su extensión espacial corresponde al área cubierta por el raster original.

Tabla 1. Características principales del objeto samples

tabla_muestra <- data.frame(
  Característica = c(
    "Clase",
    "Geometría",
    "Dimensiones",
    "Extensión",
    "Sistema de referencia",
    "Atributo",
    "Tipo de dato",
    "Valores"
  ),
  Descripción = c(
    "Objeto vectorial del paquete terra (SpatVector).",
    "Está compuesto por puntos.",
    "Contiene 500 puntos y 1 atributo (SOC).",
    "Corresponde a los límites geográficos del Eje Cafetero donde se distribuyen las muestras.",
    "WGS84 (latitud y longitud).",
    "soc_igh_15_30: representa los valores de carbono orgánico del suelo.",
    "Numérico.",
    "Incluye valores de SOC; algunos registros presentan NA o espacio en blanco por ausencia de datos."
  )
)

kable(tabla_muestra)
Característica Descripción
Clase Objeto vectorial del paquete terra (SpatVector).
Geometría Está compuesto por puntos.
Dimensiones Contiene 500 puntos y 1 atributo (SOC).
Extensión Corresponde a los límites geográficos del Eje Cafetero donde se distribuyen las muestras.
Sistema de referencia WGS84 (latitud y longitud).
Atributo soc_igh_15_30: representa los valores de carbono orgánico del suelo.
Tipo de dato Numérico.
Valores Incluye valores de SOC; algunos registros presentan NA o espacio en blanco por ausencia de datos.

Necesitamos convertir el objeto spatVector en un objeto de característica simple:

muestras <- sf::st_as_sf(samples)

muestras_tabla <- data.frame(
  soc_igh_15_30 = muestras$soc_igh_15_30,
  geometry = "<sf_POINT>"
)

DT::datatable(
  muestras_tabla,
  rownames = TRUE,
  options = list(
    dom = "tp",
    pageLength = 10,
    ordering = FALSE
  )
) %>%
  DT::formatRound(
    columns = "soc_igh_15_30",
    digits = 5
  )
nmuestras <- na.omit(muestras)

Vamos a visualizar las 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)

DT::datatable(
  sitios,
  rownames = FALSE,
  options = list(
    dom = "tp",
    pageLength = 10,
    ordering = FALSE
  )
) %>%
  DT::formatRound(
    columns = c("longit", "latit", "soc"),
    digits = 5
  )

Eliminamos los valores de NA, en este caso espacios en blanco:

Tabla 2. Sitios de muestreo sin valores NA.

sitios <- na.omit(sitios)

rownames(sitios) <- NULL
sitios$id <- seq_len(nrow(sitios))

DT::datatable(
  head(sitios),
  rownames = FALSE,
  options = list(
    dom = "tp",
    pageLength = 10,
    ordering = FALSE
  ),
) %>%
  DT::formatRound(
    columns = c("longit", "latit", "soc"),
    digits = 5
  )

Vamos a visualizar las muestras

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())

Ahora estamos listos para realizar las tareas de interpolación.

5. Interpolación

5.1 Creación del objeto gstat

Para interpolar, primero necesitamos crear un objeto de la clase gstat, utilizando la función del mismo nombre: gstat. Un objeto gstat contiene toda la información necesaria para realizar la interpolación espacial, a saber:

• La definición del modelo

• Los datos de calibración

Según sus argumentos, la función gstat “entiende” qué tipo de modelo de interpolación queremos usar:

• Sin modelo de variograma → IDW

• Modelo de variograma, sin covariables → Kriging ordinario

Vamos a usar tres parámetros de la función gstat:

• Fórmula: la fórmula de predicción que especifica las variables dependientes e independientes (también llamadas covariables)

• Datos: los datos de calibración (también llamados datos de entrenamiento)

• Modelo: el modelo de variograma

5.2 Interpolación ISW

Para interpolar el SOC utilizando el método IDW, creamos el siguiente objeto gstat, especificando únicamente la fórmula y los datos:

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

Ahora que nuestro modelo de interpolación g1 está definido, podemos usar la función predict para interpolar, es decir, para estimar los valores de precipitación.

La función predict acepta:

• Un objeto rasterizado (estrellas)

• Un modelo: objeto gstat, como g1

La imagen rasterizada tiene dos propósitos:

• Especificar las ubicaciones donde queremos realizar predicciones (en todos los métodos).

• Especificación de valores de covariables (solo en Kriging Universal)

Vamos a crear un objeto ráster con valores de celda iguales a 1:

# una simple copia
rrr = aggregate(geog.soc, 4)

¿Qué es rrr?

rrr
## class       : SpatRaster
## size        : 193, 185, 1  (nrow, ncol, nlyr)
## resolution  : 0.008834949, 0.008834949  (x, y)
## extent      : -76.25632, -74.62185, 4.075514, 5.780659  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## source(s)   : memory
## name        : soc_igh_15_30
## min value   :     11.999339
## max value   :     224.85436

Definir nuevos valores:

values(rrr) <-1

Definir nuevos nombres:

names(rrr) <- "valor"

¿Qué es rrr ahora?

rrr
## class       : SpatRaster
## size        : 193, 185, 1  (nrow, ncol, nlyr)
## resolution  : 0.008834949, 0.008834949  (x, y)
## extent      : -76.25632, -74.62185, 4.075514, 5.780659  (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)

Por ejemplo, la siguiente expresión interpola los valores de SOC según el modelo definido en g1 y la plantilla ráster definida en stars.rrr:

## [Interpolación ponderada por distancia inversa]
z1 = predict(g1, stars.rrr)
## [inverse distance weighted interpolation]

¿Qué 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  13.22538 33.03615 45.79884 46.57871 57.35566 216.9049     0
## var1.var         NA       NA       NA      NaN       NA       NA 35705
## dimension(s):
##   from  to offset     delta                       refsys x/y
## x    1 185 -76.26  0.008835 +proj=longlat +datum=WGS8... [x]
## y    1 193  5.781 -0.008835 +proj=longlat +datum=WGS8... [y]

Ahora obtenemos los nombres de los dos atributos incluidos en el objeto z1.

Podemos seleccionar solo el primer atributo y renombrarlo como “soc”:

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

Vamos a necesitar una paleta de colores.

paleta <- colorNumeric(
  palette = c("orange", "yellow", "cyan", "green"),
  domain = z1$soc,
  na.color = "transparent"
)

El ráster SOC interpolado, utilizando IDW, se muestra en la siguiente figura:

soc_idw <- as.vector(z1[[1]])

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

m <- leaflet() %>%
  addTiles() %>%
  leafem::addGeoRaster(
    z1,
    opacity = 0.7,
    colorOptions = leafem::colorOptions(
      palette = c("orange", "yellow", "cyan", "green"),
      domain = range(soc_idw, na.rm = TRUE)
    )
  ) %>%
  addMarkers(
    lng = sitios$longit,
    lat = sitios$latit,
    popup = sitios$soc,
    clusterOptions = markerClusterOptions()
  ) %>%
  addLegend(
    "bottomright",
    pal = paleta,
    values = soc_idw,
    bins = seq(20, 220, by = 20),
    title = "IDW SOC interpolation [%]"
  )

m

5.3 Interpolación

Los métodos de Kriging requieren un modelo de variograma. El modelo de variograma es una forma objetiva de cuantificar el patrón de autocorrelación en los datos y asignar ponderaciones en consecuencia al realizar predicciones.

Como primer paso, podemos calcular y examinar el variograma empírico utilizando la función de variograma.

La función requiere dos argumentos:

• fórmula: especifica la variable dependiente y las covariables, al igual que en gstat

• datos: la capa de puntos con la variable dependiente y las covariables como atributos de punto

Por ejemplo, la siguiente expresión calcula el variograma empírico de muestras, sin covariables:

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

Vamos a trazar el variograma:

plot(v_emp_ok)

Existen varias formas de ajustar un modelo de variograma a un variograma empírico. Utilizaremos la más sencilla: el ajuste automático mediante la función autofitVariogram del paquete automap.

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

La función selecciona el tipo de modelo que mejor se ajusta y también optimiza sus parámetros. Puede usar show.vgms() para visualizar los tipos de modelos de variograma.

Tenga en cuenta que la función “autofitVariogram” no funciona con objetos “sf”, por lo que convertimos el objeto a un “SpatialPointsDataFrame” (paquete “sp”).

El modelo ajustado se puede graficar con plot:

plot(v_mod_ok)

El objeto resultante es en realidad una lista con varios componentes, incluyendo el variograma empírico y el modelo de variograma ajustado. El componente $var_model del objeto resultante contiene el modelo propiamente dicho:

Tabla 3. Parámetros del modelo de variograma ajustado

tabla_variograma <- as.data.frame(v_mod_ok$var_model)

DT::datatable(
  tabla_variograma,
  rownames = FALSE,
  options = list(
    dom = "tp",
    pageLength = 5,
    ordering = FALSE
  ),
) %>%
  DT::formatRound(
    columns = c("psill", "range", "kappa"),
    digits = 3
  )

En tu cuaderno, explica el significado de cada elemento del modelo anterior.

Interpretación de la figura

En el variograma se observa que:

Los círculos representan el variograma experimental calculado a partir de los datos. La línea azul corresponde al modelo Stable ajustado automáticamente. La semivarianza aumenta con la distancia, indicando que las muestras cercanas presentan valores de SOC más parecidos que las lejanas. A partir de una distancia cercana a 44, la curva comienza a estabilizarse alrededor de un valor de 603, lo que indica que ya no existe dependencia espacial importante entre las muestras.

¿Qué significan los números sobre cada punto?

En la figura aparecen números como 41, 97, 151, 421, 499, 3057, 10976, etc.

Esos números no son valores del SOC. Representan el número de pares de muestras utilizados para calcular la semivarianza en cada intervalo de distancia (bin). Por ejemplo:

• 41 → la primera clase de distancia se calculó usando 41 pares de puntos. • 10976 → esa clase de distancia utilizó 10 976 pares de puntos.

Cuantos más pares tenga un punto del variograma, más confiable suele ser la estimación de esa semivarianza.

Ahora, el modelo de variograma se puede pasar a la función gstat y podemos continuar con la interpolación de Kriging ordinario:

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

Nuevamente, seleccionaremos un subconjunto del atributo de valores predichos y lo renombraremos:

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

Las predicciones del Kriging Ordinario se muestran en la siguiente figura:

# Convertir los valores del objeto stars en un vector numérico
soc_ok <- as.vector(z2[[1]])

# Crear la paleta de colores
paleta <- colorNumeric(
  palette = c("orange", "yellow", "cyan", "green"),
  domain = range(soc_ok, na.rm = TRUE),
  na.color = "transparent"
)

# Crear el mapa
m <- leaflet() %>%
  addTiles() %>%
  leafem::addGeoRaster(
    z2,
    opacity = 0.7,
    colorOptions = leafem::colorOptions(
      palette = c("orange", "yellow", "cyan", "green"),
      domain = range(soc_ok, na.rm = TRUE)
    )
  ) %>%
  addMarkers(
    lng = sitios$longit,
    lat = sitios$latit,
    popup = sitios$soc,
    clusterOptions = markerClusterOptions()
  ) %>%
  addLegend(
    position = "bottomright",
    pal = paleta,
    values = soc_ok,
    bins = seq(20, 220, by = 20),
    title = "OK SOC interpolation [%]"
  )

m

6. Evaluación de resultados

6.1 Evaluación cualitativa

Otra perspectiva de los tres resultados de la interpolación:

colores <- leafem::colorOptions(
  palette = c("orange", "yellow", "cyan", "green"),
  domain = range(z1$soc, z2$soc, stars.soc[[1]], na.rm = TRUE),
  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 [%]"
)
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
m

6.2 Validación cruzada

Hemos estimado superficies climáticas utilizando dos métodos diferentes: IDW y Kriging ordinario. Si bien es útil examinar y comparar los resultados gráficamente, también necesitamos una forma objetiva de evaluar la precisión de la interpolación. De esta manera, podemos elegir objetivamente el método más preciso entre los métodos de interpolación disponibles.

En la validación cruzada de exclusión de un elemento, nosotros:

  • Resta un punto de los datos de calibración.

  • Haz una predicción para ese punto.

  • Repetir para todos los puntos. Al final, lo que obtenemos es una tabla con un valor observado y un valor predicho para todos los puntos.

Podemos ejecutar la validación cruzada de exclusión de un elemento utilizando la función gstat.cv, que acepta un objeto gstat.

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

El resultado de gstat.cv tiene los siguientes atributos:

• Var1.pred — Valor predicho

• Var1.var — Varianza (solo para Kriging)

• Observed — Valor observado

• Residual — Observado - Predicho

• Zscore — Puntuación Z (solo para Kriging)

• Fold — ID de validación cruzada

cv1 = na.omit(cv1)
cv1
## class       : SpatialPointsDataFrame 
## features    : 472 
## extent      : -76.24859, -74.64504, 4.085453, 5.777346  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 6
## names       :        var1.pred, var1.var,         observed,          residual, zscore, fold 
## min values  : 15.4146598586264,       NA, 12.8426065444946, -33.4410159770908,     NA,    1 
## max values  : 94.0013415568955,       NA, 224.055267333984,  136.534663001667,     NA,  472

Vamos a convertir el objeto cv1:

cv1 = st_as_sf(cv1)

Ahora, vamos a graficar los residuos.

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

Ahora, calculamos índices de precisión de predicción, como el error cuadrático medio (RMSE):

# Este es el valor RMSE para la interpolación IDW.
sqrt(sum((cv1$var1.pred - cv1$observed)^2) / nrow(cv1))
## [1] 14.87041

Ahora, repita el proceso con los resultados correctos:

Tiempo de conversión:

cv2 = st_as_sf(cv2)

Calcula el RSME para obtener resultados OK:

# Este es el valor RMSE para la interpolación OK.
sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))
## [1] 12.45652

La validación cruzada permitió evaluar objetivamente la precisión de los métodos de interpolación mediante el cálculo del error cuadrático medio (RMSE). El método IDW obtuvo un RMSE de 14.87, mientras que el Kriging Ordinario presentó un RMSE de 12.46. Debido a que un menor RMSE indica una mayor precisión en las predicciones, se concluye que el Kriging Ordinario fue el método más preciso para estimar el carbono orgánico del suelo en el área de estudio. Esto sugiere que, al incorporar la estructura de dependencia espacial mediante el variograma, el Kriging logró representar mejor la variabilidad espacial de los datos que el método IDW.

7. Referencias

Lizarazo, I., 2023. Interpolación espacial del carbono orgánico del suelo. Disponible en https://rpubs.com/ials2un/soc_interp.

sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Tahoe 26.5.2
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] DT_0.34.0      knitr_1.51     curl_7.1.0     dplyr_1.2.1    ggplot2_4.0.3 
##  [6] leafem_0.2.5   leaflet_2.2.3  automap_1.1-20 gstat_2.1-6    stars_0.7-2   
## [11] abind_1.4-8    sf_1.1-1       terra_1.9-34   sp_2.2-3      
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6       xfun_0.58          bslib_0.11.0       raster_3.6-32     
##  [5] htmlwidgets_1.6.4  lattice_0.22-7     vctrs_0.7.3        tools_4.5.2       
##  [9] crosstalk_1.2.2    generics_0.1.4     parallel_4.5.2     tibble_3.3.1      
## [13] proxy_0.4-29       spacetime_1.3-3    xts_0.14.2         pkgconfig_2.0.3   
## [17] KernSmooth_2.23-26 RColorBrewer_1.1-3 S7_0.2.2           lifecycle_1.0.5   
## [21] compiler_4.5.2     farver_2.1.2       FNN_1.1.4.1        codetools_0.2-20  
## [25] htmltools_0.5.9    class_7.3-23       sass_0.4.10        yaml_2.3.12       
## [29] pillar_1.11.1      jquerylib_0.1.4    classInt_0.4-11    cachem_1.1.0      
## [33] tidyselect_1.2.1   digest_0.6.39      fastmap_1.2.0      grid_4.5.2        
## [37] cli_3.6.6          magrittr_2.0.5     base64enc_0.1-6    e1071_1.7-17      
## [41] withr_3.0.2        scales_1.4.0       rmarkdown_2.31     otel_0.2.0        
## [45] zoo_1.8-15         png_0.1-9          evaluate_1.0.5     rlang_1.2.0       
## [49] Rcpp_1.1.2         glue_1.8.1         DBI_1.3.0          rstudioapi_0.18.0 
## [53] reshape_0.8.10     jsonlite_2.0.0     R6_2.6.1           plyr_1.8.9        
## [57] intervals_0.15.5   units_1.0-1