1. Introduccion

Este cuaderno ilustra dos técnicas de interpolación espacial: Ponderación por Distancia Inversa (IDW) y Kriging Ordinario (OK).

Ambas técnicas se utilizan en este cuaderno para obtener una superficie continua de Carbono Orgánico del Suelo (SOC) en la profundidad de 15 a 30 cm a partir de muestras puntuales obtenidas del modelo SoilGrids 250m para el departamento de Nariño.

2. Setup

Limpiamos el espacio de trabajo

rm(list=ls())

Para la lectura y escritura de datos geoespaciales utilizamos terra (para datos ráster) y sf (para datos vectoriales). Para la interpolación espacial utilizaremos las librerías gstat y automap. Para la visualización de mapas utilizaremos ggplot2, leaflet y leafem.

Llamamos las librerias que acabmos de instalar:

library(terra)
library(sf)
library(sp)
library(stars)
library(gstat)
library(automap)
library(RColorBrewer)
library(leaflet)
library(leafem)

3. Lectura de datos de entrada

Como se indicó en la introducción, nuestras muestras corresponden a valores de Carbono Orgánico del Suelo (SOC) a una profundidad de 15 a 30 cm, obtenidos desde la plataforma SoilGrids 250 m. Guardamos estas muestras en formato GeoPackage en el directorio datos bajo el nombre soc_narino.gpkg.

Primero, verificamos que los datos estén disponibles para su uso en la carpeta de trabajo:

list.files(path="datos", pattern = "*.gpkg")
## [1] "MUNICIPIOS_9377.gpkg" "soc_narino.gpkg"

Luego, leemos las muestras puntuales de Carbono Orgánico del Suelo para el departamento de Nariño utilizando el paquete sf:

samples <- sf::st_read("datos/soc_narino.gpkg")
## Reading layer `soc_narino' from data source 
##   `C:\Users\Acer\Desktop\GB2R\InterpolacionR\datos\soc_narino.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 2000 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -77.72822 ymin: 0.7230375 xmax: -75.73364 ymax: 2.719939
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

Como nuestra área de estudio corresponde al departamento de Nariño, leemos también el archivo GeoPackage que contiene la delimitación espacial de sus municipios:

munic   <- sf::st_read("datos/MUNICIPIOS_9377.gpkg")
## Reading layer `Municipios' from data source 
##   `C:\Users\Acer\Desktop\GB2R\InterpolacionR\datos\MUNICIPIOS_9377.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 64 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 4330516 ymin: 1598084 xmax: 4573425 ymax: 1855793
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional

REPROYECTAR municipios a WGS84 (EPSG:4326) para que coincida con Leaflet y las muestras

munic <- sf::st_transform(munic, 4326)

4. Exploración de los datos de entrada

A continuación, realizamos un análisis exploratorio preliminar de los valores de Carbono Orgánico del Suelo (SOC) muestreados en el departamento de Nariño para comprender su distribución estadística.

Primero, obtenemos un resumen estadístico general de las muestras:

summary(samples)
##       soc                    geom     
##  Min.   :  0.00   POINT        :2000  
##  1st Qu.: 36.05   epsg:NA      :   0  
##  Median : 50.66   +proj=long...:   0  
##  Mean   : 56.76                       
##  3rd Qu.: 67.83                       
##  Max.   :228.91

Posteriormente, visualizamos la distribución de frecuencias de la variable soc mediante un histograma:

hist(samples$soc, 
     main = "Histograma de SOC en Nariño", 
     xlab = "Carbono Orgánico del Suelo (g/kg)", 
     ylab = "Frecuencia", 
     col = "lightblue", 
     border = "white")

Redondeamos los valores de SOC a dos decimales para facilitar su visualización y lectura:

pal <- colorNumeric(
  palette = c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),
  domain = samples$soc,
  na.color = "transparent"
)

Ahora generamos el mapa:

# Creamos el objeto del mapa
m <- leaflet() %>% 
  addTiles() %>% 
  addPolygons(data = munic, color = "gray", opacity = 1, weight = 1, fillOpacity = 0.2) %>% 
  addCircleMarkers(data = samples, radius = 1.5, label = ~soc, color = ~pal(soc), fillOpacity = 1, stroke = FALSE) %>% 
  addLegend(pal = pal, values = samples$soc, position = "bottomleft", title = "soc (g/kg)", opacity = 0.9)

m

5. Interpolación Espacial

5.1 Creación del objeto gstat

Para realizar la interpolación espacial, primero necesitamos crear un objeto de clase gstat utilizando la función del mismo nombre: gstat().

Un objeto gstat contiene toda la información necesaria para conducir la interpolación espacial, a saber: * La definición del modelo * Los datos de calibración

A partir de sus argumentos, la función gstat “comprende” qué tipo de modelo de interpolación deseamos utilizar: * Sin modelo de variograma \(\rightarrow\) Ponderación por Distancia Inversa (IDW). * Con modelo de variograma, sin covariables \(\rightarrow\) Kriging Ordinario (OK).

Utilizaremos tres parámetros principales de la función gstat: * formula: La “fórmula” de predicción que especifica las variables dependiente e independiente (o covariables). * data: Los datos de calibración (datos de entrenamiento/muestra). * model: El modelo de variograma (aplicable para Kriging).

5.2 Interpolación IDW

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

g1 <- gstat(formula = soc ~ 1, data = samples)

Ahora que nuestro modelo de interpolación g1 está definido, podemos usar la función predict para interpolar formalmente y estimar los valores de SOC en puntos no muestreados.

La función predict acepta:

Un objeto ráster o stars (como un DEM o una grilla objetivo).

Un objeto modelo de gstat (como nuestro g1).

El ráster objetivo cumple dos propósitos:

Especificar las ubicaciones donde queremos realizar las predicciones en el espacio.

Especificar valores de covariables (en el caso de Kriging Universal).

Creamos un objeto SpatRaster con la extensión de Nariño como grilla base para la interpolación:

# Definir la grilla optimizada (más liviana)
bb <- sf::st_bbox(munic)

rrr <- terra::rast(
  xmin = bb$xmin, xmax = bb$xmax, 
  ymin = bb$ymin, ymax = bb$ymax, 
  nrows = 100, ncols = 100, 
  vals = 1, crs = "epsg:4326"
)

rrr
## class       : SpatRaster
## size        : 100, 100, 1  (nrow, ncol, nlyr)
## resolution  : 0.0217653, 0.0232255  (x, y)
## extent      : -79.01021, -76.83368, 0.3613481, 2.683898  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s)   : memory
## name        : lyr.1
## min value   :     1
## max value   :     1

Convertimos el SpatRaster a un objeto de clase stars:

stars.rrr <- stars::st_as_stars(rrr)

Ejecutamos la interpolación IDW sobre la grilla utilizando el modelo g1:

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

z1 es un objeto de clase stars que contiene los resultados del proceso de interpolación espacial IDW. Posee dos atributos principales: var1.pred (los valores predichos de SOC en \(g/kg\)) y var1.var (la varianza de la predicción, que en el caso de IDW es NA debido a que es un método determinístico).

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean  3rd Qu.     Max.  NA's
## var1.pred  10.21678 57.46141 59.22222 61.89519 64.07025 183.7912     0
## var1.var         NA       NA       NA      NaN       NA       NA 10000
## dimension(s):
##   from  to offset    delta refsys x/y
## x    1 100 -79.01  0.02177 WGS 84 [x]
## y    1 100  2.684 -0.02323 WGS 84 [y]

Manejamos los valores NA asignando \(0.0\) e inspeccionamos el objeto resultante:

a <- which(is.na(z1[[1]]))
z1[[1]][a] = 0.0

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean  3rd Qu.     Max.  NA's
## var1.pred  10.21678 57.46141 59.22222 61.89519 64.07025 183.7912     0
## var1.var         NA       NA       NA      NaN       NA       NA 10000
## dimension(s):
##   from  to offset    delta refsys x/y
## x    1 100 -79.01  0.02177 WGS 84 [x]
## y    1 100  2.684 -0.02323 WGS 84 [y]

Renombramos el atributo principal a soc:

names(z1) = "soc"

Creamos una paleta de colores para graficar la superficie interpolada:

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

Visualizamos la superficie interpolada por IDW junto con las muestras puntuales:

# Mapa de la interpolación IDW
m <- leaflet() %>% 
  addTiles() %>% 
  leafem::addGeoRaster(
    z1, 
    opacity = 0.7, 
    colorOptions = colorOptions(
      palette = c("orange", "yellow", "cyan", "green"), 
      domain = range(z1$soc, na.rm = TRUE)
    )
  ) %>% 
  addCircleMarkers(
    data = samples, 
    radius = 1.5, 
    label = ~soc, 
    color = ~paleta(soc), 
    fillOpacity = 1, 
    stroke = FALSE
  ) %>% 
  addLegend(
    "bottomright", 
    pal = paleta, 
    values = z1$soc, 
    title = "IDW SOC interpolation [g/kg]"
  )

m

5.3 Interpolación Kriging Ordinario (OK)

Los métodos de Kriging requieren un modelo de variograma. El modelo de variograma cuantifica el patrón de autocorrelación espacial en los datos y asigna pesos en consecuencia para realizar las predicciones.

Como primer paso, calculamos el variograma empírico utilizando la función variogram(), la cual requiere dos argumentos principales: * formula: Especifica la variable dependiente y las covariables (al igual que en gstat). * data: La capa de puntos con los atributos.

# Cálculo del variograma empírico
v_emp_ok = variogram(soc ~ 1, data = samples)
# Graficar el variograma empírico
plot(v_emp_ok)

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

# Ajuste automático del modelo de variograma
v_mod_ok = automap::autofitVariogram(soc ~ 1, as(samples, "Spatial"))

La función selecciona el tipo de modelo con mejor ajuste y optimiza sus parámetros (Sill, Range y Nugget).

Visualizamos el variograma experimental con la curva del modelo ajustado:

# Graficar variograma empírico y modelo ajustado
plot(v_mod_ok)

El objeto resultante es una lista con varios componentes. El componente $var_model contiene los parámetros del modelo de variograma ajustado:

# Inspeccionar el modelo de variograma ajustado
v_mod_ok$var_model

Ahora pasamos el modelo de variograma a la función gstat y ejecutamos la interpolación por Kriging Ordinario:

# Definir modelo gstat con variograma e interpolar
g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2 = predict(g2, stars.rrr)
## [using ordinary kriging]
# Manejo de NA y cambio de nombre
a <- which(is.na(z2[[1]]))
z2[[1]][a] = 0.0
names(z2) = "soc"

Visualizamos

# Definir paleta dinámica adaptada a los valores predichos por Kriging
paleta_ok <- colorNumeric(
  palette = c("orange", "yellow", "cyan", "green"), 
  domain = range(z2$soc, na.rm = TRUE), 
  na.color = "transparent"
)

# Mapa de la interpolación OK
m_ok <- leaflet() %>% 
  addTiles() %>% 
  leafem::addGeoRaster(
    z2, 
    opacity = 0.7, 
    colorOptions = colorOptions(
      palette = c("orange", "yellow", "cyan", "green"), 
      domain = range(z2$soc, na.rm = TRUE)
    )
  ) %>% 
  addCircleMarkers(
    data = samples, 
    radius = 1.5, 
    label = ~soc, 
    color = ~paleta_ok(soc), 
    fillOpacity = 1, 
    stroke = FALSE
  ) %>% 
  addLegend(
    "bottomright", 
    pal = paleta_ok, 
    values = z2$soc, 
    title = "OK SOC interpolation [g/kg]"
  )

m_ok

6. Evaluación de resultados

6.1 Evaluación cualitativa

Comparamos visualmente las superficies resultantes de los métodos IDW y Kriging Ordinario (OK) superponiendo ambas capas en un mapa interactivo con selector de capas:

# Definir opciones de color para las capas
colores <- colorOptions(
  palette = c("orange", "yellow", "cyan", "green"), 
  domain = 10:100, 
  na.color = "transparent"
)
# Crear mapa de comparación cualitativa con control de capas
m_comp <- 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 = z1$soc,
    title = "Soil organic carbon [%]"
  )

m_comp

6.2 Validación Cruzada (Cross-validation)

Hemos estimado las superficies de SOC mediante dos métodos: IDW y Kriging Ordinario (OK). Para evaluar la precisión de la interpolación de manera objetiva, utilizamos el método de Validación Cruzada Dejando Uno Fuera (Leave-One-Out Cross-Validation - LOOCV) con la función gstat.cv().

En LOOCV se realiza el siguiente procedimiento: * Se excluye un punto de los datos de calibración. * Se realiza la predicción para ese punto específico. * Se repite el proceso para todos los puntos muestreados.

Ejecutamos la validación cruzada para el objeto de Kriging (g2):

#cv2 <- gstat.cv(g2, nmax = 10)

El resultado de gstat.cv() incluye atributos clave:var1.pred: Valor predicho.var1.var: Varianza de la predicción (solo en Kriging).observed: Valor observado/real.residual: Residuo (\(Observado - Predicho\)).zscore: Z-score (solo en Kriging).fold: ID de la validación cruzada.Convertimos el resultado cv2 a un objeto espacial sf:

# 1. Validación Cruzada para Kriging (OK)
#cv2 <- gstat.cv(g2, nmax = 10)
#rmse_ok <- sqrt(sum((cv2$var1.pred - cv2$observed)^2, na.rm = TRUE) / nrow(cv2))

# 2. Validación Cruzada para IDW
#cv1 <- gstat.cv(g1, nmax = 10)
#rmse_idw <- sqrt(sum((cv1$var1.pred - cv1$observed)^2, na.rm = TRUE) / nrow(cv1))

# 3. Ver los resultados
#cat("RMSE Kriging:", rmse_ok, "\nRMSE IDW:", rmse_idw)

Tras correr los códigos en la consola tenemos que: RMSE Kriging: 12.2708 g/kg RMSE IDW: 15.66138 g/kg

Se calcularon los errores cuadráticos medios (RMSE) para ambos métodos de interpolación mediante Leave-One-Out Cross Validation (LOOCV):

Análisis de Resultados y Conclusión

Al comparar los dos métodos de interpolación espacial para la estimación SOC en el departamento de Nariño, se observa una diferencia significativa en la precisión predictiva:

  1. Comparación de precisión: El método de Kriging Ordinario (OK) obtuvo un error cuadrático medio menor (RMSE = 12.27) en comparación con el método IDW (RMSE = 15.66), representando una reducción clara en el margen de error de estimación.

  2. ¿Cuál es el método más confiable y por qué? Kriging Ordinario es el método más confiable. Esto porque IDW es un modelo determinístico inflexible que da pesos basándose en la distancia geográfica predefinida. En cambio, Kriging Ordinario es un método geoestadístico óptimo e insesgado que evalúa la autocorrelación espacial real de los datos a través del modelado del variograma empírico. Esto le permite ajustar ponderaciones acordes a la variabilidad espacial propia de la variable ambiental (SOC), minimizando estadísticamente la varianza del error en las predicciones.

7. Guardado de los resultados de interpolación

Exportamos los rásters resultantes de las interpolaciones en formato GeoTIFF dentro de la carpeta data (o datos según tu estructura):

write_stars(z1, dsn = "datos/IDW_soc_stder.tif", layer = 1)
# Exportar ráster de Kriging Ordinario
write_stars(z2, dsn = "datos/OK_soc_stder.tif", layer = 1)

8. Conclusiones

En este laboratorio se aplicaron y compararon dos técnicas de interpolación espacial deterministicas y geoestadísticas (IDW y Kriging Ordinario) para mapear la distribución de Carbono Orgánico del Suelo (SOC) en el departamento de Nariño. A través de la exploración visual y cuantitativa mediante validación cruzada (LOOCV), se evidenció cómo el Kriging Ordinario logra modelar de forma continua y suave la variabilidad espacial aprovechando la autocorrelación modelada en el variograma empírico. El correcto manejo y alineación del Sistema de Referencia de Coordenadas (CRS) entre las capas vectoriales y ráster resultó fundamental para asegurar la superposición y precisión de las predicciones en el territorio.

9. Referencias

Lizarazo, I. 2025. Spatial interpolation. Available at: https://rpubs.com/ials2un/Spatial_Interpolation

sessionInfo()
## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Spanish_Mexico.utf8  LC_CTYPE=Spanish_Mexico.utf8   
## [3] LC_MONETARY=Spanish_Mexico.utf8 LC_NUMERIC=C                   
## [5] LC_TIME=Spanish_Mexico.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-2        abind_1.4-8        sp_2.2-1          
##  [9] sf_1.1-1           terra_1.9-34      
## 
## 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-7     magrittr_2.0.5     digest_0.6.39      evaluate_1.0.5    
##  [9] grid_4.5.2         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          crosstalk_1.2.2   
## [17] scales_1.4.0       codetools_0.2-20   jquerylib_0.1.4    cli_3.6.6         
## [21] rlang_1.2.0        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         FNN_1.1.4.1       
## [29] raster_3.6-32      tools_4.5.2        parallel_4.5.2     dplyr_1.2.1       
## [33] ggplot2_4.0.3      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       lifecycle_1.0.5   
## [41] classInt_0.4-11    htmlwidgets_1.6.4  pkgconfig_2.0.3    pillar_1.11.1     
## [45] bslib_0.10.0       gtable_0.3.6       glue_1.8.1         Rcpp_1.1.1-1.1    
## [49] tidyselect_1.2.1   tibble_3.3.1       xfun_0.57          rstudioapi_0.18.0 
## [53] knitr_1.51         farver_2.1.2       htmltools_0.5.9    rmarkdown_2.31    
## [57] xts_0.14.2         compiler_4.5.2     S7_0.2.2