1 1. Introducción

Este cuaderno aplica dos tecnicas de interpolacion espacial: Ponderación en funcion inversa a la distancia (IDW) Kriging ordinaria (OK).

El proposito es realizar interpolacion del Carbono organico del suelo (COS) a profunidad de 15-30cm en el area correspondiente al area del departamento del cordoba.

2 2. Configuración

# Uncomment at first-time run
rm(list=ls())
# DO NOT INSTALL  PACKAGES THAT YOU ALREADY INSTALLED
#paquetes = c('terra', 'sf', 'sp', 'stars', 'gstat', 'automap', 'leaflet', 'leafem')
#install.packages(paquetes)
#
library(terra)
## terra 1.8.70
library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(sp)
library(stars)
## Loading required package: abind
library(gstat)
library(automap)
library(RColorBrewer)
library(leaflet)
library(leafem)

3 2. Lectura de los datos de entrada.

Describir datos de entrada.

list.files(path="data", pattern = "*.gpkg")
## [1] "mun_cordoba.gpkg" "soc_cordoba.gpkg"
samples <- sf::st_read("data/soc_cordoba.gpkg")
## Reading layer `soc_cordoba' from data source 
##   `C:\Users\jsmal\OneDrive\Documents\G2\Project6\data\soc_cordoba.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1682 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -76.52913 ymin: 7.35109 xmax: -74.70232 ymax: 9.449183
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

En otra oportunidad leer los municios con el archivo gpkg

munic <- sf::st_read("data/mun_cordoba.gpkg")
## Reading layer `col_adm2' from data source 
##   `C:\Users\jsmal\OneDrive\Documents\G2\Project6\data\mun_cordoba.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 26 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.5189 ymin: 7.3516 xmax: -74.7731 ymax: 9.448473
## Geodetic CRS:  WGS 84

4 4. Exploración de datos de entrada

summary(samples)
##       soc                    geom     
##  Min.   : 10.25   POINT        :1682  
##  1st Qu.: 19.33   epsg:NA      :   0  
##  Median : 26.10   +proj=long...:   0  
##  Mean   : 30.40                       
##  3rd Qu.: 36.22                       
##  Max.   :138.03

Ahora se realiza valida el histograma

hist(samples$soc)

Redondeemos los valores de COS a dos dígitos:

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

Ahora, visualicemos las muestras para conocer su distribución espacial. Puedes visitar la documentación de leaflet para entender cómo funciona esta librería.

pal <- colorNumeric(
  c("#FF9999", "#E84A4A", "#AB3A3A", "#660000", "#1C0000"),
  # colors depend on the count variable
  domain = samples$soc,
  )
# This a simple visualization
# Change the code to suit your preferences
leaflet() %>%
  addPolygons(
    data = munic,
    color = "gray",
    # set the opacity of the outline
    opacity = 1,
    # set the stroke width in pixels
    weight = 1,
    # set the fill opacity
    fillOpacity = 0.2) %>%
 addCircleMarkers(
    data = samples,
    radius= 1.5, 
    label = ~soc,
    color = ~pal(soc),
    fillOpacity = 1,
    stroke = F
  ) %>%
  addLegend(
    data = samples,
    pal = pal,
    values = ~soc,
    position = "bottomleft",
    title = "COS:",
    opacity = 0.9) %>%
  addProviderTiles("OpenStreetMap")

5 5. Interpolacion.

5.1 5.1 Creacion del objeto gstat

Para interpolar primero necesitamos crear un objeto que tenga la clase gstat, usando una función con el mismo nombre gstat.

El objeto gstat contiene toda la información necesaria para la interpolación espacial:

  • El modelo de la función
  • La calibración de la información

Basado en estos argumentos, la función gstat determina qué tipo de interpolación usar:

  • Sin modelo de variograma → Inverso a la distancia (IDW)
  • Modelo de variograma → Kriging Ordinario (OK)

5.2 5.2Empezamos con la interpolacion IDW

En función de la muestra vamos a usar los datos de las variables.

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

Cuando vamos a interpolar usaremos la función llamada predict, que requiere un raster donde se organizan los datos en una cuadrícula. Sin esta cuadrícula, no se puede determinar dónde ubicar los valores interpolados.

El modelo es el objeto gstat que creamos anteriormente.

El raster define dónde se deben colocar los datos interpolados. Sin una referencia espacial dentro del departamento, necesitamos proporcionar el DEM (Modelo Digital de Elevación) del departamento como base.

DEM <- rast("data/DEMCordoba.tif")
DEM
## class       : SpatRaster 
## size        : 250, 209, 1  (nrow, ncol, nlyr)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -76.51667, -74.775, 7.358333, 9.441667  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : DEMCordoba.tif 
## name        : DEMCordoba 
## min value   :         -1 
## max value   :       2069
(rrr <- terra::rast(xmin=-76.51667, xmax=-74.775, ymin=7.358333, ymax=9.441667, nrows=150, ncols=125, vals=1, crs="epsg:4326"))
## class       : SpatRaster 
## size        : 150, 125, 1  (nrow, ncol, nlyr)
## resolution  : 0.01393336, 0.01388889  (x, y)
## extent      : -76.51667, -74.775, 7.358333, 9.441667  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : lyr.1 
## min value   :     1 
## max value   :     1

Nota: Si tu computadora tiene recursos limitados (por ejemplo, poca memoria RAM), usa una resolución espacial más baja (por ejemplo, 0.008 o 0.01 grados).

Ahora necesitamos convertir este SpatRaster en un objeto stars:

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

Ahora realizamos la interpolación con IDW (Ponderación Inversa de la Distancia):

## [inverse distance weighted interpolation]
## running this code takes several minutes
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  11.89306 23.40692 28.61628 30.57707 34.91851 131.3214     0
## var1.var         NA       NA       NA      NaN       NA       NA 18750
## dimension(s):
##   from  to offset    delta refsys x/y
## x    1 125 -76.52  0.01393 WGS 84 [x]
## y    1 150  9.442 -0.01389 WGS 84 [y]

La variable Var1.var representa la varianza. En IDW es un método determinista y estima que no tiene ningún error asociado.

# dealing with NA values
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  11.89306 23.40692 28.61628 30.57707 34.91851 131.3214     0
## var1.var         NA       NA       NA      NaN       NA       NA 18750
## dimension(s):
##   from  to offset    delta refsys x/y
## x    1 125 -76.52  0.01393 WGS 84 [x]
## y    1 150  9.442 -0.01389 WGS 84 [y]
names(z1) = "soc"

Para trazar la superficie interpolada, necesitamos una paleta de colores:

# change this chunk to meet your preferences
paleta <- colorNumeric(palette = c("#FF9999", "#E84A4A", "#AB3A3A", "#660000", "#1C0000"), domain = 10:140, na.color = "transparent")

Ahora, visualicemos el resultado de la interpolación:

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 1,                
      colorOptions = colorOptions(palette = c("#FF9999", "#E84A4A", "#AB3A3A", "#660000", "#1C0000"), 
                                  domain = 10:140)
    ) %>%
   addCircleMarkers(
    data = samples,
    radius= 1.5, 
    label = ~soc,
    color = ~paleta(soc),
    fillOpacity = 1,
    stroke = F
  ) %>%
    addLegend("bottomright", pal=paleta, values= z1$soc,
    title = "IDW SOC interpolation [%]"
    )
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m  # Print the map

6 5.3 Interpolación OK (Kriging Ordinario)

Los métodos de kriging requieren un modelo de variograma. Este modelo 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 variogram.

Esta función requiere dos argumentos:

formula: especifica la variable dependiente y las covariables, al igual que en gstat.

data: la capa de puntos con la variable dependiente y las covariables como atributos de punto.

La siguiente expresión calcula el variograma empírico de muestras sin covariables:

v_emp_ok = variogram(soc ~ 1, data=samples)

Grafiquemos el variograma empírico:

plot(v_emp_ok)

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

#make sure you understand the parameters needed to run this fitting function
v_mod_ok = automap::autofitVariogram(soc ~ 1, as(samples, "Spatial"))

La función selecciona el tipo de modelo más adecuado y ajusta sus parámetros. Puede usar show.vgms()

para visualizar varios tipos de modelos de variograma.

Tenga en cuenta que la función autofitVariogram no funciona con objetos sf, por lo que convertimos el objeto de muestras en 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:

v_mod_ok$var_model
##   model     psill    range kappa
## 1   Nug  34.34964   0.0000   0.0
## 2   Ste 381.81423 189.4159   0.3

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

## [using ordinary kriging]
## This takes several minutes 
## (or even ages depending on your raster cell size and your computer resources)
g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2= predict(g2, stars.rrr)
## [using ordinary kriging]

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

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

Las predicciones correctas se muestran en el siguiente gráfico:

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 1,                
      colorOptions = colorOptions(palette = c("#FF9999", "#E84A4A", "#AB3A3A", "#660000", "#1C0000"), domain = 10:140)
    ) %>%
   addCircleMarkers(
    data = samples,
    radius= 1.5, 
    label = ~soc,
    color = ~paleta(soc),
    fillOpacity = 1,
    stroke = F
  ) %>%
    addLegend("bottomright", pal = paleta, values= z2$soc,
    title = "OK SOC interpolation [%]"
    )
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m  # Print the map

7 6.Assessment of results

8 6.1 Qualitative assessment

colores <- colorOptions(palette = c("#FF9999", "#E84A4A", "#AB3A3A", "#660000", "#1C0000"), domain = 10:100, na.color = "transparent")
m <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(z1, colorOptions = colores, opacity = 0.8, group= "IDW")  %>%
  addGeoRaster(z2, colorOptions = colores, opacity = 0.8, group= "OK")  %>%
  # Add layers controls
  addLayersControl(
    overlayGroups = c("IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
    addLegend("bottomright", pal = paleta, values= z1$soc,
    title = "Soil organic carbon [%]"
)
## 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  # Print the map

9 6.2 Validación cruzada

Hemos estimado las superficies de SOC mediante 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 dejar uno fuera (Leave-One-Out Cross Validation), hacemos lo siguiente:

Extraemos un punto de los datos de calibración.

Realizamos una predicción para ese punto.

Repetimos el proceso para todos los puntos. Al final, obtenemos una tabla con un valor observado y un valor predicho para todos los puntos.

Podemos ejecutar la validación cruzada de dejar uno fuera (Leave-One-Out Cross Validation) utilizando la función gstat.cv, que acepta un objeto gstat.

Al escribir el siguiente fragmento de código, oculte el mensaje y los resultados.

El resultado de gstat.cv tiene los siguientes atributos:

var1.pred—Valor predicho

var1.var—Varianza (solo para Kriging)

observed—Valor observado

residual—Valor observado - Valor predicho

zscore—Puntuación Z (solo para Kriging)

fold—ID de validación cruzada

Convirtamos el objeto cv2:

Grafiquemos los residuales:

#Análisis de las diferencias visualizadas en el Carbono Orgánico del Suelo (COS)

##1. Análisis del gráfico de residuales (OK.residuals - Imagen 1) En el gráfico de residuales del Kriging Ordinario observo:

Residuales negativos (puntos magenta/morados): Indican que el modelo sobreestimó el COS en esas ubicaciones (el valor predicho fue mayor que el observado). El residual más bajo es -39.924. Residuales positivos (puntos verdes): Indican que el modelo subestimó el COS (el valor predicho fue menor que el observado). El residual más alto es 70.94. Distribución espacial: Los residuales parecen estar distribuidos de manera relativamente aleatoria en el espacio, lo cual es deseable. Sin embargo, hay algunos clusters de residuales positivos grandes (verdes grandes) en ciertas zonas, lo que sugiere que el modelo tiene dificultades para capturar la variabilidad en esas áreas específicas. Magnitud: La presencia de residuales grandes (especialmente el 70.94) indica que hay puntos donde el modelo OK tiene errores considerables.

##2. Análisis de la superficie interpolada (Imagen 2) En el mapa interpolado observo:

Zonas con bajo COS (tonos verdes claros/amarillos): Predominan en la parte central del departamento de Córdoba. Zonas con alto COS (tonos naranjas/rojizos): Se concentran principalmente en:

Esquina inferior izquierda (suroeste) Algunas áreas dispersas en el norte Zonas específicas en el este

Patrón espacial suave: El Kriging produce una superficie más suavizada en comparación con IDW, debido a que considera la estructura de autocorrelación espacial.

##3. Principales diferencias entre IDW y Kriging Ordinario

#IDW (Ponderación Inversa de la Distancia):

Produce superficies con efecto “ojo de buey” alrededor de los puntos muestreados No considera la estructura espacial de autocorrelación Es un método determinista (no proporciona medidas de incertidumbre) Tiende a crear transiciones más bruscas entre valores altos y bajos

#Kriging Ordinario (OK):

Produce superficies más suaves y realistas Considera el variograma (estructura de autocorrelación espacial) Es un método estocástico (proporciona varianza de predicción) Mejor interpolador lineal insesgado (BLUP) Respeta mejor los patrones naturales de distribución del COS

10 7. Guardar los resultados de la interpolación

write_stars(
  z1, dsn = "data/IDW_soc_stder.tif", layer = 1
)
write_stars(
  z2, dsn = "data/OK_soc_stder.tif", layer = 1
)

11 8. Conclusiones

Este trabajo representa un ejercicio aplicado de geoestadística para el mapeo del Carbono Orgánico del Suelo (COS) en el departamento de Córdoba, Colombia, a una profundidad de 15-30 cm. A través de la implementación comparativa de dos metodologías de interpolación espacial —Ponderación Inversa de la Distancia (IDW) y Kriging Ordinario (OK)—, se lograron generar superficies predictivas que revelan patrones espaciales significativos en la distribución del COS en la región. La aplicación del Kriging Ordinario demostró ser técnicamente superior a IDW por varias razones fundamentales. Primero, el análisis variográfico permitió cuantificar objetivamente la estructura de autocorrelación espacial del COS, identificando el rango de influencia y el grado de dependencia espacial entre las observaciones. Esta información, capturada mediante el modelo de variograma ajustado, proporciona una base teóricamente robusta para la asignación óptima de pesos durante la interpolación. Segundo, el Kriging no solo genera predicciones puntuales sino que también cuantifica la incertidumbre asociada a cada estimación mediante la varianza de kriging, lo que resulta invaluable para la toma de decisiones informadas en manejo de suelos y planificación territorial. La validación cruzada Leave-One-Out (LOOCV) reveló que el Kriging Ordinario presentó un Error Cuadrático Medio (RMSE) menor en comparación con IDW, confirmando su mayor precisión predictiva. Los estudios comparativos entre estos métodos han demostrado consistentemente que, aunque ambos interpoladores pueden tener desempeños similares en ciertos contextos, el Kriging generalmente es superior al predecir la variación espacial de propiedades del suelo RStudio

12 9. Reference

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

sessionInfo()
## R version 4.5.1 (2025-06-13 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=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.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-4        stars_0.6-8        abind_1.4-8        sp_2.2-0          
##  [9] sf_1.0-21          terra_1.8-70      
## 
## 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-7          digest_0.6.37          
##  [7] magrittr_2.0.4          evaluate_1.0.5          grid_4.5.1             
## [10] fastmap_1.2.0           plyr_1.8.9              jsonlite_2.0.0         
## [13] e1071_1.7-16            reshape_0.8.10          DBI_1.2.3              
## [16] crosstalk_1.2.2         scales_1.4.0            codetools_0.2-20       
## [19] jquerylib_0.1.4         cli_3.6.5               rlang_1.1.6            
## [22] units_1.0-0             intervals_0.15.5        base64enc_0.1-3        
## [25] cachem_1.1.0            yaml_2.3.10             FNN_1.1.4.1            
## [28] raster_3.6-32           tools_4.5.1             parallel_4.5.1         
## [31] dplyr_1.1.4             ggplot2_4.0.0           spacetime_1.3-3        
## [34] 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        
## [40] classInt_0.4-11         leaflet.providers_2.0.0 htmlwidgets_1.6.4      
## [43] pkgconfig_2.0.3         pillar_1.11.1           bslib_0.9.0            
## [46] gtable_0.3.6            glue_1.8.0              Rcpp_1.1.0             
## [49] tidyselect_1.2.1        tibble_3.3.0            xfun_0.53              
## [52] rstudioapi_0.17.1       knitr_1.50              farver_2.1.2           
## [55] htmltools_0.5.8.1       rmarkdown_2.30          xts_0.14.1             
## [58] compiler_4.5.1          S7_0.2.0