---

title: “Interpolacion Espacial del carbono organico del suelo en Meta”

author: “Nawel Sebastian Ballesteros Sierra”

date: 13/11/2025

##1.Introduccion

Este cuaderno muestra cómo dos métodos de interpolación—Inverse Distance Weighting (IDW) y Kriging Ordinario (OK)—permiten transformar datos puntuales de carbono orgánico del suelo (SOC) en superficies continuas. Usando información de SoilGrids (15–30 cm), se ilustra la diferencia entre un método basado únicamente en la distancia (IDW) y uno que incorpora la estructura espacial del fenómeno (OK).

##2.Configuracones

Limpiar el entorno

Elimina todos los objetos cargados en la sesión para evitar errores con variables viejas.

rm(list=ls())

Carga de librerías

Estas librerías permiten: terra / sf / sp / stars → trabajar con datos espaciales (raster y vector). gstat / automap → realizar interpolaciones (IDW, Kriging). leaflet / leafem → crear mapas interactivos. RColorBrewer → generar paletas de colores.

library(terra)

## Warning: package 'terra' was built under R version 4.5.2

## terra 1.8.86

library(sf)

## Warning: package 'sf' was built under R version 4.5.2

## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE

library(sp)

## Warning: package 'sp' was built under R version 4.5.2

library(stars)

## Warning: package 'stars' was built under R version 4.5.2

## Cargando paquete requerido: abind

## Warning: package 'abind' was built under R version 4.5.2

library(gstat)

## Warning: package 'gstat' was built under R version 4.5.2

library(automap)

## Warning: package 'automap' was built under R version 4.5.2

library(RColorBrewer)

## Warning: package 'RColorBrewer' was built under R version 4.5.2

library(leaflet)

## Warning: package 'leaflet' was built under R version 4.5.2

library(leafem)

## Warning: package 'leafem' was built under R version 4.5.2

##3.Cargar archivos

Revisar si existe el archivo de municipios en la carpeta data.

list.files(path="data", pattern = "metamunicipios.gpkg")

## [1] "metamunicipios.gpkg"

Se carga los puntos con el valor de SOC (carbono orgánico del suelo).

samples <- sf::st_read("data/soc_meta.gpkg")

## Reading layer `soc_meta' from data source 
##   `C:\Users\nawel\Downloads\GB2\Proyecto 6\data\soc_meta.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1947 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.93382 ymin: 1.570154 xmax: -71.03937 ymax: 4.861022
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

Se carga el mapa de municipios del Meta.

munic <- sf::st_read("data/metamunicipios.gpkg")

## Reading layer `MetaMunicipios' from data source 
##   `C:\Users\nawel\Downloads\GB2\Proyecto 6\data\metamunicipios.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 29 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.9206 ymin: 1.570401 xmax: -71.0749 ymax: 4.860201
## Geodetic CRS:  WGS 84

##4.Revisar y explorar los datos

summary(): estadísticas descriptivas del SOC.

summary(samples)

##       soc                    geom     
##  Min.   : 8.525   POINT        :1947  
##  1st Qu.:12.215   epsg:NA      :   0  
##  Median :15.183   +proj=long...:   0  
##  Mean   :20.770                       
##  3rd Qu.:20.957                       
##  Max.   :82.989

hist(): histograma de la distribución del SOC.

hist(samples$soc)

round(): redondea los valores a 2 decimales.

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

Se crea una paleta de colores para Leaflet donde se define los colores que representarán los valores de SOC.

pal <- colorNumeric(
  c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),
  # colors depend on the count variable
  domain = samples$soc,
  )

Mapa interactivo de puntos en los que se muestra: polígonos de municipios, puntos con valor SOC,leyenda,fondo de OpenStreetMap

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 = "SOC:",
    opacity = 0.9) %>%
  addProviderTiles("OpenStreetMap")

##5.Interpolaciones

Se le indica a gstat que queremos interpolar la variable soc sin variables auxiliares.

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

Se carga el raster de altitud del departamento del meta para generar una malla donde se puedan interpolar los datos

dem <- terra::rast("rastermeta.tif")

Crear el raster vacío donde se interpolará

Este raster define: área de interpolación (coordenadas de meta),dimensiones(resolución),sistema de referencia El valor “1” es solo de relleno.

(rrr <- terra::rast(xmin=-74.9, xmax=-71.1, ymin=1.62, ymax=4.91, nrows=370, ncols = 329, vals = 1, crs = "epsg:4326"))

## class       : SpatRaster 
## size        : 370, 329, 1  (nrow, ncol, nlyr)
## resolution  : 0.01155015, 0.008891892  (x, y)
## extent      : -74.9, -71.1, 1.62, 4.91  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : lyr.1 
## min value   :     1 
## max value   :     1

Convertir raster a formato STARS es necesario porque gstat trabaja con STARS para predicción.

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

Se hace una Interpolación IDW al generar un raster interpolado

z1 = predict(g1, stars.rrr)

## [inverse distance weighted interpolation]

se verifica que 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  8.605225 13.93298 16.14701 20.74495 21.65034 82.4907      0
## var1.var         NA       NA       NA      NaN       NA      NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329  -74.9   0.01155 WGS 84 [x]
## y    1 370   4.91 -0.008892 WGS 84 [y]

Rellenando valores NA ya que algunos no tendran valores por lo cual Los huecos sin datos los reemplaza por 0.

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

Nuevamente se verifica si estamos trabajando con Z1

z1

## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean  3rd Qu.    Max.   NA's
## var1.pred  8.605225 13.93298 16.14701 20.74495 21.65034 82.4907      0
## var1.var         NA       NA       NA      NaN       NA      NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329  -74.9   0.01155 WGS 84 [x]
## y    1 370   4.91 -0.008892 WGS 84 [y]

se asigna el nombre a z1 por soc

names(z1) = "soc"

se genera la paleta de colores y

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

Se agrega el raster interpolado al mapa Leaflet con raster IDW

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 10:100)
    ) %>%
   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 pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA

## Warning in paleta(soc): Some values were outside the color scale and will be
## treated as NA
## Warning in paleta(soc): 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 

Variograma empírico

Analiza la autocorrelación espacial del SOC.

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

generar el diagrama

plot(v_emp_ok)

Ajuste automático del variograma

Encuentra el mejor modelo teórico (Spherical, Exponential, etc.).

v_mod_ok = automap::autofitVariogram(soc ~ 1, as(samples, "Spatial"))

Genera el variograma teorico

plot(v_mod_ok)

v_mod_ok$var_model

##   model     psill    range kappa
## 1   Nug  10.69485   0.0000   0.0
## 2   Ste 281.31428 264.5914   0.4

Crea el modelo gstat con variograma

g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2= predict(g2, stars.rrr)

## [using ordinary kriging]

Rellena valores NA los huecos sin datos los reemplaza por 0.

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

Mapa Leaflet con Kriging

Igual que con IDW, pero usando z2.

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 10:100)
    ) %>%
   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 pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA

## Warning in paleta(soc): Some values were outside the color scale and will be
## treated as NA
## Warning in paleta(soc): 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 

colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), domain = 10:100, na.color = "transparent")

Comparación de IDW y Kriging con controles de capasy permite activar y desactivar cada raster.

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 pal(c(r[1], cuts, r[2])): 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
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored

Generamos el mapa comparativo entre IDW Y Kriking

m

Exportar los resultados

Guardamos los raster interpolados en archivos TIFF.

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

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

Una vez calculados ambos valores RMSE, compárelos y analice sus diferencias. ¿Cuál es el método de interpolación más fiable? Explique.

En interpolación espacial, el RMSE (Root Mean Square Error) mide la diferencia entre los valores observados y los estimados por cada modelo. Es decir:

RMSE alto → peor precisión RMSE bajo → mejor precisión

El método más fiable es el que presente menor RMSE, y en la mayoría de casos es Kriging, porque modela explícitamente la autocorrelación espacial y produce estimaciones óptimas.

##8.Conclusiones

En este cuaderno se realizó un flujo completo de análisis geoestadístico para estimar el SOC en el Meta, desde la carga y exploración de datos hasta la creación de mapas interactivos y la generación de interpolaciones. Se aplicaron y compararon dos métodos, IDW y Kriging Ordinario, entendiendo sus diferencias y la importancia del variograma en la modelación espacial. La comparación permitió identificar que Kriging suele ofrecer resultados más consistentes al incorporar la autocorrelación espacial. Finalmente, se exportaron los productos interpolados, logrando una comprensión general del proceso de interpolación aplicado a datos ambientales.

##9.Referencias

Tomado desde

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 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Spanish_Latin America.utf8 
## [2] LC_CTYPE=Spanish_Latin America.utf8   
## [3] LC_MONETARY=Spanish_Latin America.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=Spanish_Latin America.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-86      
## 
## 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            FNN_1.1.4.1             raster_3.6-32          
## [28] tools_4.5.1             parallel_4.5.1          dplyr_1.1.4            
## [31] ggplot2_4.0.0           spacetime_1.3-3         png_0.1-8              
## [34] vctrs_0.6.5             R6_2.6.1                zoo_1.8-14             
## [37] proxy_0.4-27            lifecycle_1.0.4         classInt_0.4-11        
## [40] leaflet.providers_2.0.0 htmlwidgets_1.6.4       pkgconfig_2.0.3        
## [43] pillar_1.11.1           bslib_0.9.0             gtable_0.3.6           
## [46] glue_1.8.0              Rcpp_1.1.0              tidyselect_1.2.1       
## [49] tibble_3.3.0            xfun_0.54               rstudioapi_0.17.1      
## [52] knitr_1.50              farver_2.1.2            htmltools_0.5.8.1      
## [55] rmarkdown_2.30          xts_0.14.1              compiler_4.5.1         
## [58] S7_0.2.0