Interpolacion Espacial del Carbón Organico del Suelo en el Departamento de Casanare

Introducción

La interpolación espacial es una rama de la geoestadística que se basa en el cálculo de los valores desconocidos de una variable spacial a partir de otros datos con valor conocido, basándose en la suposición de que los fenomenos geograficos cercanos se encuentran mas relacionados entre si que los lejanos. [figura 1] QGIS Proyec.(2017). Esta herramienta permite modelar la distribucion de variables, por ejemplo, el soc (soil carbo organic), que funciona como indicador para conocer la calidad y fertilidad en el suelo de estudio.
En este cuaderno se aplican metodos de análissi de interpolacion espacial usando el lenguaje de programación R para analisar y visualizar la distribucion de soc en el departamento de Casanare. En el se desarrolló y comparo tres modelos espaciales; El primero fue basado en valores reales muestreados (Real World), el segundo se realizo con la interpolación por distancia inversa ponderada (IDW), y por ultimo se usó Kriging Ordinary (OK).

[Figura 1]

rm(list=ls())

Entorno de trabajo y herramientas

Durante el desarrollo de este análisis se usó diferentes herramientas provenientes de los paquetes de trabajo de R, en el desarrollo se manipularón datos vectoriales y ráster con ayuda de las librerias (sf, sp, terra, stars), para la interpolación geoestadistica (gstat y automap), y finalmente para la visualizacion de mapas interactivos se usó las librerias de (ggplot2, leaflet y leafem).

library(sp)
library(terra)
library(sf)
library(stars)
library(gstat)
library(automap)
library(leaflet)
library(leafem)
library(ggplot2)
library(dplyr)

curl permite la descarga de datos desde fuentes externas. Con la integración de estos paquetes se puede desarrollar un análisis de mejor calidad y exascitud.

library("curl")

## Using libcurl 8.14.1 with Schannel

Preparación de datos

Aquí se usa el paquete curl, se configura un “handle” personalizado que es compatible con los servidores que operan bajo el protoclo HTTP/2.

h <- new_handle()
handle_setopt(h, http_version = 2)

archivo <- ("C:/Users/Janus/Desktop/GB2/p7/CASANARE/soc_igh_15_30.tif")
(soc <- rast(archivo))

## class       : SpatRaster 
## size        : 917, 1413, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8145750, -7792500, 477250, 706500  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source      : soc_igh_15_30.tif 
## name        : soc_igh_15_30

Los valores del ráster son reescalados para expresar el SOC en porcentaje, y la capa que es reproycta al sistema de coordenadas WGS84, para que se acompatible con las herramientas de análisis y visualización espacial.

soc.perc  <- soc/10

geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(soc.perc, geog))

## class       : SpatRaster 
## size        : 934, 1489, 1  (nrow, ncol, nlyr)
## resolution  : 0.002205488, 0.002205488  (x, y)
## extent      : -73.09924, -69.81527, 4.286671, 6.346597  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_igh_15_30 
## min value   :      6.961041 
## max value   :    116.768845

Visualisación

Se convierte el archivo raster descargado de soc, a un objeto stars, para que sea compatible con leafem.

stars.soc = st_as_stars(geog.soc)

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

Extracción de puntos para la interpolación

El conjunto de datos espaciales que se visualizan estan representados por un SpatVctor de tipo puntos, estan dimensionado a 1000 geometrias o puntos, asociados a un atributo ; Que en este caso es a soc_igh_15_30, la extension geografica indinca la delimitacion de los puntos que se estan localizando

set.seed(123456)
(samples <- spatSample(geog.soc, 1000, "random", as.points=TRUE))

##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 1000, 1  (geometries, attributes)
##  extent      : -73.09814, -69.81638, 4.287774, 6.345495  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : soc_igh_15_30
##  type        :         <num>
##  values      :         10.59
##                        18.32
##                        11.98

El objeto samples fue generado como un conjunto de puntos espaciales con spatSample, con este comado lo convertimos a un objeto sf, que es un formato mas compatible para trabajar con datos vectoriles.

(muestras <- sf::st_as_sf(samples))

Es necesario verificar que todos los datos no tengan valores faltantes en las variables para que los datos de interpolación sean validos, para eso se uso la variable na.omit, que elimina u omite esos valores.

nmuestras <- na.omit(muestras)

longit <- st_coordinates(muestras)[,1]
latit <- st_coordinates(muestras)[,2]
soc <- muestras$soc_igh_15_30

se creo un identificador numerico para el conjunto de puntos creando una secuencia que va desde el 1 hasta el 1000, adoptando un punto cada identificador numerico.

id <- seq(1,1000,1)

El data.frame agrupa 4 vectores que se usaron como molde para diferentes modelos de interpolación

(sitios <- data.frame(id, longit, latit, soc))

Por último se realiza un chequeo final de todos los datos, asegurando que no existan valores nulos que puedan afectar al modelo de interpolación

sitios <- na.omit(sitios)

head(sitios)

Visualización de puntos de muestreo

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

El mapa interativo permite observar el raster del soc con los puntos seleccionados, se representan en las etiquetas flotantes.

Definición del modelo previo a la interpolacción

La función gstat define un modelo geoestadistico en el que la variable respuesta son los valores muestreados de soc, la modelacion de los datos establece una relacion espacial necesaria para la interpolación.

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

aggregate reduce la resolución del raster agrupando celdas del raster en bloques mas gandes, para efectos de este análisis se agruparon en factor de agregación 4, es decir cada grupo de 4*4 son iguales a 16 celdas del rater original

rrr = aggregate(geog.soc, 4)

rrr es el nuevo raster con la resolucion reducida, esto facilita el análisis computacional y lo hace mas eficiente, pero las visualizaciones son menos detalladas. Hijmans, R. J., & van Etten, J. (2023).

rrr

## class       : SpatRaster 
## size        : 234, 373, 1  (nrow, ncol, nlyr)
## resolution  : 0.008821954, 0.008821954  (x, y)
## extent      : -73.09924, -69.80866, 4.28226, 6.346597  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_igh_15_30 
## min value   :      8.061533 
## max value   :    103.183873

Los valores en la celdas del raster se modifican a 1, esto para que puedan usarse como base en la interpolación.

values(rrr) <-1

se asigna el nombre valor a la unica capa del raster.

names(rrr) <- "valor"

rrr ahora es un raster donde cada celda tiene un valor de 1,que cubre la misma areas espacial que el original pero con menor resolución

rrr

## class       : SpatRaster 
## size        : 234, 373, 1  (nrow, ncol, nlyr)
## resolution  : 0.008821954, 0.008821954  (x, y)
## extent      : -73.09924, -69.80866, 4.28226, 6.346597  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : valor 
## min value   :     1 
## max value   :     1

Kriging Ordinario

El kriging es una procedimiento geoestadistico avanzado que genera una superficie estimada a partir de un conjunto de punstos dispersos con valores z. (esri. 2025). Este metodo asume que los datos presentan autocorrelación espacial,los puntos cercanos son mas parecidos que los lejanos; se basa tanto en la distancia entre cada punto como en el grado de variacion espacial.

[Figura 2. (esri. 2025)]

Se paso el raster de entrada rrr a una estructura start, se aplico la interpolación espacial OK definida en g1, generando una predicción espacial

stars.rrr = st_as_stars(rrr)
z1 = predict(g1, stars.rrr)

## [inverse distance weighted interpolation]

z1 es un raster con los nuevos valores estimados de la interpolación.

z1

## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean 3rd Qu.     Max.  NA's
## var1.pred  8.828741 15.67117 17.25569 20.98964 19.7485 87.51642     0
## var1.var         NA       NA       NA      NaN      NA       NA 87282
## dimension(s):
##   from  to offset     delta                       refsys x/y
## x    1 373  -73.1  0.008822 +proj=longlat +datum=WGS8... [x]
## y    1 234  6.347 -0.008822 +proj=longlat +datum=WGS8... [y]

este codigo extrae la capa var1.pred donde etan los valores calculados por el kriging, y lo renombra a soc

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

se usa una paleta que permita visualizar en este caso el carbono orgánico en el suelo, de manera intuitiva y fácil de interpretar por los lectores.

paleta <- colorNumeric(palette = c("#ffffcc", "#a1dab4", "#41b6c4", "#2c7fb8", "#253494"), domain = 10:100, na.color = "transparent")

visualización IDW

La herramienta de distancia inversa ponderada (IDW) se conoce como metodo deterministico de interpolación porque se basan directamente en los valores medidos.

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("#ffffcc", "#a1dab4", "#41b6c4", "#2c7fb8", "#253494"), 
                                  domain = 11:55)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
    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

m 

El variograma empírico es una herramienta que permite medir la correlación de una variable según la distancia entre los puntos

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

plot(v_emp_ok)

La función autofitVariogram prueba diferentes modelos y ajusta automáticamente a los datos del variograma, esto se utiliza para realizar las interpolaciones en el modelo kriging.

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

plot(v_mod_ok)

v_mod_ok$var_model

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

## [using ordinary kriging]

Z2 tiene las capas extraidas de la prediccion evaluada por kriging, y se seleciona unicamente var1.pred y ,, mantiene las dimensiones espaciales

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

Visualización Kriging

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("#ffffcc", "#a1dab4", "#41b6c4", "#2c7fb8", "#253494"), 
                                  domain = 11:55)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
    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

m 

Se pueden ver resultados más precisos de la variable que se está estudiando

colores <- colorOptions(palette = c("#ffffcc", "#a1dab4", "#41b6c4", "#2c7fb8", "#253494"), domain = 10:100, na.color = "transparent")

Interpolacion de soc usando modelos de : Datos Reales vs IDW vs Kriging

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 

Para fines de este trabajo, este mapa nos permite evaluar la interpolación con los métodos implementados y compararlos con valores reales.

cv1 = na.omit(cv1)

cv1

## class       : SpatialPointsDataFrame 
## features    : 954 
## extent      : -73.07609, -69.86049, 4.287774, 6.345495  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 6
## names       :        var1.pred, var1.var,         observed,          residual, zscore, fold 
## min values  : 9.76613431796397,       NA, 8.74386882781982, -22.4448776265792,     NA,    1 
## max values  : 76.1968498250655,       NA, 88.4289169311523,  37.4075628973918,     NA,  954

cv1 = st_as_sf(cv1)

Mapa de burbujas residual del modelo IDW

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

cv2 = na.omit(cv2)

cv2 = st_as_sf(cv2)

mapa de burbujas residual del modelo Kriging

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

sqrt(sum((cv1$var1.pred - cv1$observed)^2) / nrow(cv1))

## [1] 5.578804

cv2 = st_as_sf(cv2)

sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))

## [1] 4.818877

conclusión

Al comparar los modelos se puede observar que Kriging tiene una mejor distribución de los datos espaciales, es mas homogénea y presenta ayor similitud con la visualización real. Esto demuestra que el modelo de Kriging realiza mejor la prediccion de los datos en los puntos muestreados, su modelo que interpreta la cercania y tambien la estructura espacial y la correlación de los dato lo hace mas confiable.

Bibliografia

Esri. (consultado el 20/07/2025). Cómo funciona Kriging. ArcGIS Pro. https://pro.arcgis.com/es/pro-app/latest/tool-reference/3d-analyst/how-kriging-works.htm

Hijmans, R. J., & van Etten, J. (2023). Aggregate raster cells or SpatialPolygons/Lines. In raster: Geographic Data Analysis and Modeling (Version 3.6-26) [R package documentation]. CRAN. https://search.r-project.org/CRAN/refmans/raster/html/aggregate.html

sessionInfo()

## R version 4.5.0 (2025-04-11 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_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] curl_6.4.0     dplyr_1.1.4    ggplot2_3.5.2  leafem_0.2.4   leaflet_2.2.2 
##  [6] automap_1.1-16 gstat_2.1-3    stars_0.6-8    abind_1.4-8    sf_1.0-21     
## [11] terra_1.8-54   sp_2.2-0      
## 
## 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-6     digest_0.6.37      magrittr_2.0.3     evaluate_1.0.3    
##  [9] grid_4.5.0         RColorBrewer_1.1-3 fastmap_1.2.0      plyr_1.8.9        
## [13] jsonlite_2.0.0     e1071_1.7-16       reshape_0.8.10     DBI_1.2.3         
## [17] crosstalk_1.2.1    scales_1.4.0       codetools_0.2-20   jquerylib_0.1.4   
## [21] cli_3.6.5          rlang_1.1.6        units_0.8-7        intervals_0.15.5  
## [25] withr_3.0.2        base64enc_0.1-3    cachem_1.1.0       yaml_2.3.10       
## [29] FNN_1.1.4.1        raster_3.6-32      tools_4.5.0        parallel_4.5.0    
## [33] spacetime_1.3-3    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    classInt_0.4-11   
## [41] htmlwidgets_1.6.4  pkgconfig_2.0.3    pillar_1.10.2      bslib_0.9.0       
## [45] gtable_0.3.6       glue_1.8.0         Rcpp_1.0.14        tidyselect_1.2.1  
## [49] tibble_3.2.1       xfun_0.52          rstudioapi_0.17.1  knitr_1.50        
## [53] farver_2.1.2       htmltools_0.5.8.1  rmarkdown_2.29     xts_0.14.1        
## [57] compiler_4.5.0

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.