1.Introducción

Este código realiza un análisis espacial del Carbono Orgánico del Suelo (SOC) en el departamento de Casanare, Colombia, utilizando interpolación espacial con IDW (Distancia Inversa Ponderada) y Kriging Ordinario.

Este análisis espacial examina la distribución del Carbono Orgánico del Suelo (SOC) en el departamento de Casanare, Colombia, utilizando interpolación espacial con los métodos de IDW (Distancia Inversa Ponderada) y Kriging Ordinario. Casanare es una región caracterizada por su diversidad geográfica, que abarca desde las planicies de la Orinoquía hasta el piedemonte llanero, lo que influye en la distribución de SOC debido a variaciones en el uso del suelo, la vegetación y las prácticas agrícolas y ganaderas.

La distribución espacial del SOC es fundamental para comprender la fertilidad del suelo, la capacidad de almacenamiento de carbono y el impacto de actividades humanas en los ecosistemas del departamento.

2. Configuración del entorno

Antes de ejecutar cualquier código, es importante limpiar el espacio de trabajo para evitar conflictos con datos anteriores:

rm(list=ls())

A continuación, se importan las bibliotecas necesarias para la lectura de datos geoespaciales, interpolación y visualización:

library(terra)
## terra 1.8.15
library(sf)
## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
library(sp)
library(stars)
## Cargando paquete requerido: abind
library(gstat)
library(automap)
library(RColorBrewer)
library(leaflet)
library(leafem)

Antes de usar estas librerías, debes instalarlas en la consola de R con install.packages() si no las tienes.

3. Carga de datos

Se verifica la existencia de archivos en la carpeta “datos”:

list.files(path="datos", pattern = "*.gpkg")
## [1] "agua_areas_2.gpkg"     "agua_lineas_2.gpkg"    "Casanare_ctities.gpkg"
## [4] "depto_2.gpkg"          "mun_casanare.gpkg"     "roads_casanare.gpkg"  
## [7] "soc_csnar.gpkg"

Esto devuelve una lista de archivos en formato GeoPackage (.gpkg) disponibles en la carpeta de datos.

Luego, se cargan los datos de SOC (carbono orgánico del suelo) desde el archivo correspondiente:

samples <- sf::st_read("datos/soc_csnar.gpkg")
## Reading layer `soc_csnar' from data source 
##   `C:\Users\Lilian\OneDrive\Documentos\GB2\Rproyect\datos\soc_csnar.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1930 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -73.07829 ymin: 4.290006 xmax: -69.82965 ymax: 6.345495
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

También se cargan los límites administrativos de los municipios de Casanare:

munic <- sf::st_read("datos/mun_casanare.gpkg")
## Reading layer `municipiosdecolombia' from data source 
##   `C:\Users\Lilian\OneDrive\Documentos\GB2\Rproyect\datos\mun_casanare.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 19 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.07777 ymin: 4.287476 xmax: -69.83591 ymax: 6.346111
## Geodetic CRS:  MAGNA-SIRGAS

4. Exploración de los datos

Para entender la distribución de los valores de SOC, se obtiene un resumen estadístico:

summary(samples)
##       soc                     geom     
##  Min.   :  8.741   POINT        :1930  
##  1st Qu.: 13.731   epsg:NA      :   0  
##  Median : 16.299   +proj=long...:   0  
##  Mean   : 20.939                       
##  3rd Qu.: 20.716                       
##  Max.   :105.884

Se genera un histograma para visualizar la distribución de los valores:

hist(samples$soc)

En el histograma, se observa que los valores de SOC en Casanare presentan una distribución asimétrica, con mayor concentración en valores bajos. Esto indica que la mayoría de los suelos de la región tienen niveles relativamente bajos de carbono orgánico, lo que podría estar relacionado con prácticas de uso del suelo o características naturales del ecosistema.

Se redondean los valores de SOC a dos decimales para facilitar su interpretación:

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

5. Visualización de los datos en un mapa

Se utiliza leaflet() para mostrar la distribución de las muestras de SOC sobre un mapa:

pal <- colorNumeric(
  c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),
  # colors depend on the count variable
  domain = samples$soc,
  )
Caracterización del Departamento de Casanare

Casanare está ubicado en la región de la Orinoquía colombiana y presenta una diversidad de paisajes que influyen directamente en sus suelos y contenido de SOC. Se pueden distinguir tres grandes unidades fisiográficas:

  • Piedemonte llanero: zona de transición entre la cordillera Oriental y los Llanos, con suelos fértiles y alta actividad agrícola.

  • Llanuras aluviales: regiones con suelos arenosos y humedales que pueden influir en la acumulación de carbono orgánico.

  • Sabanas y morichales: áreas de vegetación natural y suelos con baja capacidad de almacenamiento de carbono debido a su uso para la ganadería extensiva.

  • Las prácticas de uso del suelo en Casanare incluyen ganadería, cultivos de arroz y palma de aceite, y explotación petrolera, lo que ha generado impactos en la calidad del suelo y su contenido de carbono. Este código define una escala de colores según los valores de SOC. Luego, se genera el mapa:

library(leaflet)

# Define a color palette for the 'soc' values (you can adjust this based on your data)
pal <- colorNumeric(c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50") , domain = samples$soc)

# This is your visualization code with the corrected `pal` definition
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 = FALSE
  ) %>%
  addLegend(
    data = samples,
    pal = pal,
    values = ~soc,
    position = "bottomleft",
    title = "SOC:",
    opacity = 0.9
  ) %>%
  addProviderTiles("OpenStreetMap")
## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'

Esto genera un mapa interactivo con marcadores de colores según la concentración de SOC en cada punto.

  • En el mapa, se pueden ver variaciones espaciales en los valores de SOC en Casanare. Algunas áreas presentan concentraciones más altas, posiblemente en suelos con mayor cobertura vegetal o menor degradación. Otras zonas muestran valores bajos, lo que podría indicar erosión o suelos con menor capacidad de retención de carbono.

  • En el piedemonte llanero, se observan valores más altos de SOC debido a la presencia de suelos más ricos en materia orgánica y menor nivel de degradación.

  • En las llanuras aluviales, los valores de SOC varían ampliamente, con zonas de alta concentración cerca de cuerpos de agua y otras con valores más bajos debido a la erosión.

  • En las sabanas y morichales, el SOC tiende a ser más bajo, ya que la actividad ganadera intensiva ha reducido el contenido de materia orgánica en el suelo.

Estos patrones reflejan la influencia del clima, el relieve y el uso del suelo en la distribución del SOC en la región.

6. Interpolación de datos

6.1 Creación del objeto gstat

Antes de interpolar, se define el objeto gstat, que almacena la información del modelo:

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

6.2 Interpolación IDW

Se crea un raster vacío que cubre el área de estudio:

(rrr <- terra::rast(xmin=-73.125, xmax=-69.609, ymin=4.214, ymax=6.664, nrows=357, ncols = 513, vals = 1, crs = "epsg:4326"))
## class       : SpatRaster 
## dimensions  : 357, 513, 1  (nrow, ncol, nlyr)
## resolution  : 0.006853801, 0.006862745  (x, y)
## extent      : -73.125, -69.609, 4.214, 6.664  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : lyr.1 
## min value   :     1 
## max value   :     1

Se convierte el raster a formato stars, que es compatible con la interpolación:

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

Se ejecuta la interpolación IDW:

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  8.922556 15.7071 17.02395 20.94176 20.10285 101.9296      0
## var1.var         NA      NA       NA      NaN       NA       NA 183141
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 513 -73.12  0.006854 WGS 84 [x]
## y    1 357  6.664 -0.006863 WGS 84 [y]

Se manejan los valores NA y se renombran los datos:

# 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  8.922556 15.7071 17.02395 20.94176 20.10285 101.9296      0
## var1.var         NA      NA       NA      NaN       NA       NA 183141
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 513 -73.12  0.006854 WGS 84 [x]
## y    1 357  6.664 -0.006863 WGS 84 [y]
names(z1) = "soc"

Se define una paleta de colores para visualizar los resultados:

# change this chunk to meet your preferences
paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), domain = 10:120, na.color = "transparent")

Se visualiza el resultado en un mapa: El método IDW asigna valores a los puntos desconocidos basándose en la proximidad de los puntos de muestreo. En la siguiente imagen se observa la interpolación obtenida con este método

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  # Print the map

Se nota una fuerte influencia de los valores más cercanos, lo que genera un patrón de zonas bien definidas pero con transiciones bruscas. Esto puede ser útil en áreas con datos densos, pero en regiones con menos estaciones puede producir resultados poco precisos.

7. Interpolación Kriging (OK)

7.1 Cálculo del variograma

Se obtiene el variograma empírico, que mide la relación espacial de los datos:

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

Se ajusta un modelo de variograma de manera automática:

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

9En el variograma, se puede observar cómo la variabilidad de los valores de SOC cambia con la distancia. Si el rango es corto, indica que los valores de SOC varían significativamente en distancias pequeñas, lo que puede estar relacionado con cambios en el uso del suelo o condiciones edafológicas locales en Casanare. Si el rango es mayor, sugiere que las diferencias en SOC son más homogéneas a lo largo de grandes extensiones de terreno.

v_mod_ok$var_model
##   model      psill    range kappa
## 1   Nug   5.203208  0.00000   0.0
## 2   Ste 165.633951 98.85009   0.4

7.2 Interpolación Kriging

Se crea el modelo de Kriging:

g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2= predict(g2, stars.rrr)
## [using ordinary kriging]
z2
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean  3rd Qu.      Max.
## var1.pred  9.123210 14.75111 16.28033 21.50055 20.19954  91.46287
## var1.var   8.646039 19.50555 23.58124 33.25783 32.16132 133.59457
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 513 -73.12  0.006854 WGS 84 [x]
## y    1 357  6.664 -0.006863 WGS 84 [y]

Se maneja los valores NA y se renombran los datos:

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

Se visualiza el resultado:

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

El método Kriging es un estimador geoestadístico que considera la estructura espacial de los datos. A continuación, se presenta el mapa generado con este método:

m  # Print the map

Se observa una interpolación más suavizada y continua, lo que da mayor realismo a la distribución espacial. Sin embargo, este método requiere un mayor conocimiento estadístico y puede ser más sensible a la calidad de los datos de entrada.

8. Comparación de resultados

Se comparan los mapas de interpolación IDW y OK en un mismo mapa interactivo:

colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), 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 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
m  # Print the map

La interpolación IDW presenta transiciones más abruptas entre valores de SOC, mientras que el método Kriging (OK) genera una superficie más suavizada. En el contexto de Casanare, el método Kriging podría ser más apropiado si se espera una variación gradual en los suelos, mientras que IDW podría ser útil para resaltar diferencias locales más marcadas. Comparando ambos mapas, se puede ver que en algunas áreas existen discrepancias, lo que sugiere la necesidad de más datos de campo para mejorar la precisión de la interpolación

8.1 Interpolación Kriging Ordinario (OK) y IDW

knitr::include_graphics("datos/residual.png")

knitr::include_graphics("datos/residualOK.png")

Comparación de Interpolaciones: OK vs. IDW

Distribución de residuos:

  • En la interpolación OK (primera imagen), los residuos parecen más equilibrados, con menos valores extremos en comparación con la interpolación IDW (segunda imagen).

  • En IDW, los residuos presentan un rango mayor de valores, con puntos de sobrestimación y subestimación más pronunciados.

  • En zonas agrícolas de Casanare, como los municipios de Yopal y Aguazul, donde la variabilidad del suelo es alta, una interpolación con menos errores es crucial para la planificación de cultivos sostenibles.

Magnitud de los errores:

  • La interpolación OK tiene un rango de residuos entre aproximadamente -16.321 y 58.982.

  • En la interpolación IDW, los residuos van de -24.907 a 44.342.

  • Aunque IDW tiene un menor valor máximo positivo, presenta una mayor subestimación (-24.907 frente a -16.321), lo que sugiere que en ciertos puntos está generando errores más altos en comparación con OK.

  • En áreas como Paz de Ariporo, donde la ganadería extensiva impacta la calidad del suelo, un método de interpolación más preciso ayudaría a evaluar mejor la compactación y degradación del suelo.

Distribución espacial de los errores:

  • En Kriging Ordinario, los residuos parecen más suavizados y menos concentrados en ciertas áreas.

  • En IDW, los errores son más pronunciados en ciertas zonas, especialmente en el cuadrante superior izquierdo, donde hay una acumulación de puntos de gran tamaño, lo que indica errores más elevados en esa región.

  • En Casanare, la heterogeneidad del suelo entre la zona de piedemonte y la llanura aluvial hace que un modelo como Kriging, que considera la estructura espacial de los datos, sea más adecuado.

Precisión de interpolación:

  • Si se toma el Error Cuadrático Medio (RMSE) como referencia, el método con el menor RMSE sería el más preciso.

  • Basándonos en la observación gráfica, OK parece ser más preciso, ya que tiene menos valores extremos en los residuos.

  • IDW, al basarse en la ponderación de la distancia sin considerar la estructura espacial de los datos, puede generar mayores errores en zonas donde hay menos puntos de muestreo.

  • En la región del río Meta y sus afluentes, donde la sedimentación y la erosión pueden alterar la distribución del SOC, el uso de un método más robusto como Kriging permitiría un análisis más detallado.

write_stars(
  z1, dsn = "datos/IDW_soc_csnar.tif", layer = 1
)
write_stars(
  z2, dsn = "datos/OK_soc_csnar.tif", layer = 1
)

9. Bibliografía

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

sessionInfo()
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## 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] leafem_0.2.3       leaflet_2.2.2      RColorBrewer_1.1-3 automap_1.1-16    
##  [5] gstat_2.1-3        stars_0.6-7        abind_1.4-8        sp_2.1-4          
##  [9] sf_1.0-19          terra_1.8-15      
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6            xfun_0.49               bslib_0.8.0            
##  [4] ggplot2_3.5.1           raster_3.6-31           htmlwidgets_1.6.4      
##  [7] lattice_0.22-6          leaflet.providers_2.0.0 vctrs_0.6.5            
## [10] tools_4.4.2             crosstalk_1.2.1         generics_0.1.3         
## [13] parallel_4.4.2          tibble_3.2.1            proxy_0.4-27           
## [16] spacetime_1.3-3         fansi_1.0.6             xts_0.14.1             
## [19] pkgconfig_2.0.3         KernSmooth_2.23-24      lifecycle_1.0.4        
## [22] compiler_4.4.2          farver_2.1.2            FNN_1.1.4.1            
## [25] munsell_0.5.1           codetools_0.2-20        htmltools_0.5.8.1      
## [28] class_7.3-22            sass_0.4.9              yaml_2.3.10            
## [31] pillar_1.9.0            jquerylib_0.1.4         classInt_0.4-10        
## [34] cachem_1.1.0            tidyselect_1.2.1        digest_0.6.37          
## [37] dplyr_1.1.4             fastmap_1.2.0           grid_4.4.2             
## [40] colorspace_2.1-1        cli_3.6.3               magrittr_2.0.3         
## [43] base64enc_0.1-3         utf8_1.2.4              e1071_1.7-16           
## [46] scales_1.3.0            rmarkdown_2.29          zoo_1.8-12             
## [49] png_0.1-8               evaluate_1.0.1          knitr_1.49             
## [52] rlang_1.1.4             Rcpp_1.0.13-1           glue_1.8.0             
## [55] DBI_1.2.3               rstudioapi_0.17.1       reshape_0.8.9          
## [58] jsonlite_1.8.9          R6_2.5.1                plyr_1.8.9             
## [61] intervals_0.15.5        units_0.8-5