1. Introducción

Este informe aplica dos técnicas de interpolación espacial, IDW (Inverse Distance Weighting) y Kriging Ordinario (OK), para estimar el contenido de carbono orgánico del suelo (SOC) en el departamento de Boyacá a una profundidad de 15-30 cm. La técnica IDW asigna valores a las ubicaciones no muestreadas basándose en la proximidad de las muestras disponibles, otorgando mayor peso a los puntos más cercanos. Por otro lado, el Kriging Ordinario (OK) no solo considera la distancia entre las muestras, sino que también incorpora la estructura espacial de los datos a través del análisis de variogramas, lo que permite obtener estimaciones más precisas y confiables. Ambas metodologías buscan generar mapas predictivos del SOC en Boyacá, proporcionando herramientas útiles para la gestión y conservación de suelos en la región.

2. Configuración

2.1. Limpiar zona de trabajo

rm(list=ls()) 
# eliminar los mensajes no necesarios
options(showPackageStartupMessages = FALSE)

2.2 Instalar Librerias

library(terra)
## Warning: package 'terra' was built under R version 4.2.3
## terra 1.7.71
library(sf)
## Warning: package 'sf' was built under R version 4.2.3
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(sp)
## Warning: package 'sp' was built under R version 4.2.3
library(stars)
## Loading required package: abind
## Warning: package 'abind' was built under R version 4.3.0
library(gstat)
## Warning: package 'gstat' was built under R version 4.2.3
library(automap)
library(RColorBrewer)
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.2.3
library(leafem)
## Warning: package 'leafem' was built under R version 4.2.3

2.3 A continuacion pondremos el enlace donde se encuentran los datos

"C:/Users/usuagro/Documents/GB2.M/Informe_2n/Datos"
## [1] "C:/Users/usuagro/Documents/GB2.M/Informe_2n/Datos"

3. Lectura de Datos

Como se mencionó en la introducción, las muestras utilizadas en este estudio corresponden a valores de Carbono Orgánico del Suelo (SOC) medidos a una profundidad de 15-30 cm. Estos datos fueron obtenidos de SoilGrids 250m, una fuente global de información edáfica, y procesados mediante este cuaderno de trabajo. Las muestras se almacenaron en formato geopackage dentro del directorio de datos, bajo el nombre de archivo soc_boyaca.gpkg, para facilitar su manejo y análisis en el entorno de trabajo, en este caso se debio utilizar la ruta completa pues si no se menjaba de esta manera no aparecia los datos de la carpeta

list.files(path="C:/Users/usuagro/Documents/GB2.M/Informe_2n/Datos", pattern = "*.gpkg") 
## [1] "carreteras.gpkg"  "ciudades.gpkg"    "Depto.gpkg"       "Depto_mapa.gpkg" 
## [5] "soc_boyaca.gpkg"  "Vias_Boyaca.gpkg"

En el siguiente código, se realiza la carga del archivo soc_boyaca.gpkg, ubicado en la carpeta Datos, utilizando la función st_read() del paquete sf. Esta función permite leer y cargar datos geoespaciales en formato geopackage de manera eficiente, facilitando su manipulación y análisis en el entorno de R.

samples <- sf::st_read("C:/Users/usuagro/Documents/GB2.M/Informe_2n/Datos/soc_boyaca.gpkg")
## Reading layer `soc_boyaca' from data source 
##   `C:\Users\usuagro\Documents\GB2.M\Informe_2n\Datos\soc_boyaca.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1912 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.6938 ymin: 4.65768 xmax: -71.95758 ymax: 7.055167
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

3.1. Lectura de los Limites Administrativos

Dado que el área de estudio es Boyacá, es fundamental cargar y utilizar los datos correspondientes a los municipios de la región para contextualizar espacialmente los resultados. Esto se lleva a cabo mediante el siguiente código:

munic <- sf::st_read("C:/Users/usuagro/Documents/GB2.M/Informe_2n/Datos/Depto.gpkg")
## Reading layer `Depto' from data source 
##   `C:\Users\usuagro\Documents\GB2.M\Informe_2n\Datos\Depto.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 123 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## Geodetic CRS:  MAGNA-SIRGAS

Este archivo contiene datos sobre los municipios ubicados dentro del área de estudio de Boyacá, lo que facilitará una mejor visualización y comprensión de la distribución espacial de los valores de Carbono Orgánico del Suelo (SOC). Consiste en una colección de entidades geográficas compuestas por múltiples polígonos, cada uno con atributos relevantes para el análisis. El sistema de referencia espacial utilizado es MAGNA-SIRGAS WGS 84, lo que garantiza la compatibilidad con otros conjuntos de datos geoespaciales empleados en este estudio.

4. Exploración de los datos: Distribución Espacial del Carbono Orgánico en Boyacà

4.1.

Grafico de Barras Para analizar y comprender mejor la distribución de los valores de Carbono Orgánico del Suelo (SOC), realizamos un resumen estadístico de las muestras disponibles. Este resumen nos brinda información detallada sobre la variabilidad del contenido de SOC en la región, incluyendo valores mínimos, máximos, media y cuartiles. Adicionalmente, generamos un histograma para visualizar de manera más clara y efectiva la distribución de los datos, utilizando el siguiente código:

summary(samples)
##       soc                    geom     
##  Min.   : 10.39   POINT        :1912  
##  1st Qu.: 26.02   epsg:NA      :   0  
##  Median : 37.08   +proj=long...:   0  
##  Mean   : 38.83                       
##  3rd Qu.: 49.11                       
##  Max.   :104.27
hist(samples$soc, breaks = 10, col = "#69b3a2", # breaks se relaciona a las columnas
     main = "Distribución de SOC", xlab = "SOC (g/kg)", ylab = "Frecuencia")

Este análisis preliminar nos ayuda a detectar posibles valores atípicos y a identificar tendencias en los datos antes de proceder con la aplicación de las técnicas de interpolación.

A partir del histograma, se puede observar que la frecuencia más alta, superior a 400, corresponde a valores de Carbono Orgánico del Suelo (SOC) alrededor de 50 g/kg, lo que indica una concentración significativa de datos en ese rango.

4.2.

Mapa Para visualizar la distribución espacial del Carbono Orgánico del Suelo (SOC) en la región de Boyacá, se empleó un mapa interactivo creado con la librería Leaflet. Esta herramienta permite representar de manera clara y dinámica la ubicación de las muestras, así como su relación con las zonas geográficas de interés.

En esta sección, el código redondea los valores de la columna SOC en el conjunto de datos samples a dos decimales. Este paso es útil para evitar mostrar números con una precisión excesiva en el mapa, lo que facilita su interpretación y mejora la claridad de la visualización.

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

A continuación, se creó una paleta de colores que va desde el azul claro (#2c7bb6) hasta el rojo intenso (#d7191c), pasando por tonos intermedios como el celeste (#abd9e9), el amarillo claro (#ffffbf) y el naranja (#fdae61). Estos colores se asignaron a los valores de SOC (contenido de carbono en el suelo) en los datos. De esta manera, los valores más bajos de SOC se representan con tonos azules claros, mientras que los valores más altos se asocian con tonos rojos intensos, permitiendo una visualización clara y gradual de la distribución del carbono en el suelo.

pal <- colorNumeric(
  c("#2c7bb6", "#abd9e9", "#ffffbf", "#fdae61", "#d7191c"),
  domain = samples$soc)

En el siguiente código, para generar el mapa, se cargaron las fronteras de los municipios de Boyacá utilizando la función addPolygons(). Se definió un borde de color gris con opacidad total para facilitar la visualización de las áreas geográficas. La opacidad del relleno de los polígonos se ajustó a un valor de 0.2, lo que permitió una visualización más ligera de las zonas sin opacar los puntos de muestra.

A continuación, se añadieron los círculos que representan las muestras utilizando la función addCircleMarkers(). Cada círculo corresponde a una muestra tomada, con un tamaño estándar de 1.5 píxeles, y se coloreó de acuerdo con el valor de SOC de la muestra mediante una paleta de colores personalizada. Los círculos se configuraron con opacidad completa y sin trazo en el borde para resaltar el color asociado a cada muestra.

Finalmente, se incluyó una leyenda en la esquina inferior izquierda del mapa usando la función addLegend(), la cual proporciona información sobre los valores de SOC asociados a cada color. Esto permite una interpretación rápida y clara de los valores representados en el mapa.

El fondo del mapa se obtuvo de OpenStreetMap mediante la función addProviderTiles(), ofreciendo un contexto geográfico adecuado para ubicar las muestras en la región de Boyacá.

leaflet(options = leafletOptions(zoomSnap = 0.5, zoomControl = TRUE)) %>%
  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(
    data = samples,
    pal = pal,
    values = ~soc,
    position = "bottomleft",
    title = "SOC:",
    opacity = 0.7
  ) %>%
  addProviderTiles("OpenStreetMap") %>%
  setView(lng = -73.5, lat = 5.3, zoom = 8)
## 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'

En este mapa se visualizan los límites del departamento de Boyacá y sus respectivos municipios, proporcionando una idea clara de la distribución de los datos. Se observa una mayor concentración de valores de Carbono Orgánico del Suelo (SOC) superiores a 100 g/kg en el área de los municipios de Tunja, Sogamoso, Duitama, Paipa, Chiquinquirá, Moniquirá y Paramo de Pisba. Esta concentración sugiere que estas zonas presentan un contenido significativamente alto de SOC, lo que podría estar relacionado con características específicas del suelo o prácticas de manejo en la región.

5. Interpolaciòn de SOC en Boyacà

5.1. Interpolaciòn IDW

Para estimar el contenido de carbono orgánico del suelo (SOC) en Boyacá, se utiliza la interpolación espacial mediante el método de Inverse Distance Weighting (IDW). Para ello, se emplea la función gstat, que permite definir un modelo de interpolación basado en datos de calibración obtenidos en puntos de muestreo. El modelo se construye con la siguiente estructura:

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

Este modelo permite predecir valores de Carbono Orgánico del Suelo (SOC) (atributo que depende de sí mismo, representado como ~1) en ubicaciones donde no se tienen datos, utilizando la información de los puntos de muestreo disponibles. La predicción se lleva a cabo mediante la función predict, la cual requiere un modelo (en este caso, g1) y un raster que define la extensión del área de interés. Para generar este raster, se utiliza la librería terra, empleando las coordenadas que delimitan la zona de estudio en Boyacá

ext <- terra::ext(samples)  # Obtener los límites geográficos de los datos de muestreo
rrr <- terra::rast(ext, resolution = 0.01, vals = 1, crs = "epsg:4326")
# Se crea un raster con la extensión de los datos, una resolución de 0.01 grados,
# valores iniciales de 1 en todas las celdas y un sistema de coordenadas WGS 84.
rrr
## class       : SpatRaster 
## dimensions  : 240, 274, 1  (nrow, ncol, nlyr)
## resolution  : 0.01, 0.01  (x, y)
## extent      : -74.6938, -71.9538, 4.65768, 7.05768  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : lyr.1 
## min value   :     1 
## max value   :     1

Este raster se convierte en un objeto stars para poder realizar la interpolación, ya que la función utilizada solo trabaja con este tipo de objeto. A través de este procedimiento, se genera una superficie estimada del Carbono Orgánico del Suelo (SOC) en la región, lo que permite analizar de manera efectiva su distribución espacial.

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

A continuación, se lleva a cabo la interpolación utilizando el modelo g1 y el raster stars.rrr a través de la función predict. Este proceso permite estimar los valores de Carbono Orgánico del Suelo (SOC) en áreas donde no se dispone de datos, generando una superficie continua que representa la distribución espacial del SOC en la región

z1 = predict(g1, stars.rrr)
## [inverse distance weighted interpolation]
z1 # imprimir
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean  3rd Qu.   Max.  NA's
## var1.pred  11.27198 31.57292 37.73126 38.39956 45.07321 97.465     0
## var1.var         NA       NA       NA      NaN       NA     NA 65760
## dimension(s):
##   from  to offset delta refsys x/y
## x    1 274 -74.69  0.01 WGS 84 [x]
## y    1 240  7.058 -0.01 WGS 84 [y]

El objeto resultante z1 incluye dos atributos principales: var1.pred, que guarda los valores estimados de Carbono Orgánico del Suelo (SOC), y var1.var, que corresponde a la varianza de la predicción. Para manejar posibles valores faltantes en los datos, estos se reemplazan por ceros, asegurando que no haya inconsistencias en el análisis posterior.

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.27198 31.57292 37.73126 38.39956 45.07321 97.465     0
## var1.var         NA       NA       NA      NaN       NA     NA 65760
## dimension(s):
##   from  to offset delta refsys x/y
## x    1 274 -74.69  0.01 WGS 84 [x]
## y    1 240  7.058 -0.01 WGS 84 [y]

De esta manera, se genera una superficie interpolada del Carbono Orgánico del Suelo (SOC) para Boyacá, la cual está lista para ser visualizada y analizada. Este resultado permite explorar y comprender la distribución espacial del SOC en la región, facilitando la identificación de patrones y tendencias clave.

names(z1) = "soc"

Para visualizar la interpolación del Carbono Orgánico del Suelo (SOC) en Boyacá, se define una paleta de colores que asigna tonalidades específicas a los valores interpolados. Esta paleta permite representar de manera clara y efectiva la variación espacial del SOC en la región, facilitando su interpretación y análisis.

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

A continuación, se genera un mapa interactivo utilizando leaflet, en el cual se superpone la capa interpolada sobre un mapa base. Se ajusta la opacidad de la interpolación para una visualización adecuada y se incorporan los puntos de muestreo, coloreados según sus valores de Carbono Orgánico del Suelo (SOC). Esto permite una representación dinámica y detallada de la distribución espacial del SOC en la región.

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.5,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "#76EEC6"), 
                                  domain = 10:100) # Tamaño
    ) %>%
   addCircleMarkers(
    data = samples,
    radius= 1.5, 
    label = ~soc,
    color = ~paleta(soc), # Paleta ya definida anteriormente
    fillOpacity = 1,
    stroke = F
  ) %>%
    addLegend("bottomright", pal=paleta, values= z1$soc,
    title = "Interpolación IDW SOC [%]"
    )
## 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

5.2.

Interpoblaciòn OK Para llevar a cabo la interpolación mediante Kriging Ordinario (OK) en Boyacá, primero es necesario definir un modelo de variograma. Este modelo permite cuantificar la autocorrelación espacial de los datos y asignar pesos adecuados en las predicciones.

El primer paso consiste en calcular el variograma empírico a partir de los puntos de muestreo. Esto se realiza utilizando la función variogram, donde se especifica la variable dependiente (en este caso, el Carbono Orgánico del Suelo, SOC) y los datos de calibración:

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

A continuación, se graficó el variograma empírico para examinar la estructura espacial de los datos. Esta representación visual permite identificar patrones de autocorrelación y determinar el modelo de variograma que mejor se ajusta a la distribución espacial del Carbono Orgánico del Suelo (SOC) en Boyacá.

plot(v_emp_ok, 
     main = "Variograma Empírico SOC",   # Título 
     xlab = "Distancia (m)",            
     ylab = "Semivarianza",            
     col = "#7B68EE",                   # Color de los puntos
     border = "darkgray"                # Color del borde               
)

Ajuste del modelo del Variograma Existen varios métodos para ajustar un modelo a un variograma empírico. En este caso, se empleó un enfoque automático utilizando la función autofitVariogram() del paquete automap. Esta función selecciona el modelo de variograma más adecuado y optimiza sus parámetros de manera automática, asegurando un ajuste óptimo a la estructura espacial de los datos.

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

Dado que la función autofitVariogram() no es compatible con objetos de tipo sf, fue necesario convertir el objeto samples en un SpatialPointsDataFrame utilizando el paquete sp.

Una vez realizado el ajuste del modelo de variograma, se utilizó la función plot() para visualizar el variograma ajustado, permitiendo una evaluación clara de la estructura espacial de los datos.

plot(v_mod_ok)

v_mod_ok$var_model
##   model     psill    range kappa
## 1   Nug  18.48519   0.0000   0.0
## 2   Ste 313.58287 121.3931   0.3

Una vez obtenido el modelo de variograma ajustado, este se incorporó en la función gstat para realizar la interpolación mediante Kriging Ordinario (OK). Este método estima los valores de la variable de interés (en este caso, el Carbono Orgánico del Suelo, SOC) en ubicaciones no muestreadas, teniendo en cuenta la dependencia espacial identificada en el variograma. De esta manera, se logra una estimación más precisa y confiable de la distribución espacial del SOC en la región.

g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2= predict(g2, stars.rrr)
## [using ordinary kriging]

Manejo de Valores Faltantes

Tras realizar la predicción, se detectaron y reemplazaron los valores NA (no disponibles) con 0.0 para garantizar que no haya problemas durante la visualización de los datos. Este paso asegura que la superficie interpolada esté completa y lista para su análisis y representación gráfica.

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

A continuación, se genera el mapa que representa el porcentaje de la Interpolación por Kriging Ordinario (OK) y los valores de Carbono Orgánico del Suelo (SOC). Este mapa permite visualizar de manera clara y detallada la distribución espacial del SOC en la región, basada en los resultados de la interpolación.

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 0.5,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "#76EEC6"), 
                                  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 = "Interpolación OK y SOC [%]"
    )
## 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  # Mostrar el mapa

6. Evaluaciòn de Resultados

6.1.

Evaluaciòn Cualitativa Para analizar las predicciones de Carbono Orgánico del Suelo (SOC) en Boyacá, se realizaron interpolaciones mediante los métodos de IDW (Inverse Distance Weighting) y Kriging Ordinario (OK), las cuales se visualizaron en un mapa interactivo usando Leaflet con una paleta de colores que incluye naranja para valores bajos de SOC, amarillo y cyan para valores intermedios, y azul marino para valores altos. Este enfoque facilita la comparación cualitativa de los patrones espaciales generados por ambos métodos en la región, y a continuación se muestra el código en R utilizado para crear el mapa.

colores <- colorOptions(palette = c("orange", "yellow", "cyan", "#76EEC6"), domain = 10:100, na.color = "transparent")
#mapa interactivo con Leaflet
m <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(z1, colorOptions = colores, opacity = 0.6, group= "IDW")  %>%
  addGeoRaster(z2, colorOptions = colores, opacity = 0.6, group= "OK")  %>%
  # Add layers controls
  addLayersControl(
    overlayGroups = c("IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
    addLegend("bottomright", pal = paleta, values= z1$soc,
    title = "Carbón Orgánico del Suelo [%]"
)
## 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

Esta comparación cualitativa permite identificar patrones espaciales y posibles anomalías en la distribución del Carbono Orgánico del Suelo (SOC), las cuales podrían estar asociadas a factores ambientales o a la densidad de muestreo en la región. Sin embargo, a partir de este punto, los códigos no pudieron ejecutarse debido a las limitaciones de capacidad del computador para procesarlos. No obstante, se incluye la explicación detallada de los procedimientos que se deben realizar, junto con los códigos correspondientes, para guiar el proceso de análisis

6.2. Validaciòn Cruzada

Esta comparación cualitativa facilita la identificación de patrones espaciales y posibles anomalías en la distribución del Carbono Orgánico del Suelo (SOC) en Boyacá, las cuales podrían estar relacionadas con factores ambientales o con la densidad de muestreo en la región. Sin embargo, debido a las limitaciones de capacidad del computador, los códigos no pudieron ejecutarse a partir de este punto. A pesar de ello, se proporciona una explicación detallada de los procedimientos que deben llevarse a cabo, junto con los códigos correspondientes, para orientar el proceso de análisis de manera adecuada

### CORRER EN LA CONSOLA:

# Validación para Kriging
# cv2 = gstat.cv(g2)

# Convertir a objeto sf
# cv2 = st_as_sf(cv2)

# Graficar residuos
# sp::bubble(as(cv2[, "OK residuals"], "Spatial"))

# Calcular RMSE para la interpolación OK
# rmse_ok = sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))
# rmse_ok

Después de calcular el RMSE (Error Cuadrático Medio) para el método de Kriging Ordinario (OK), se repitió el proceso para los resultados de IDW (Inverse Distance Weighting) en Boyacá. La comparación de ambos valores de RMSE permitió evaluar de manera objetiva la precisión de las predicciones:

Valores más bajos de RMSE indican una mayor precisión predictiva.

Valores más altos de RMSE sugieren mayores discrepancias entre los valores observados y los predichos.

Este análisis comparativo ayuda a determinar cuál de los dos métodos ofrece estimaciones más confiables para la distribución del Carbono Orgánico del Suelo (SOC) en la región.

Comparaciòn y Analisis

Si el RMSE para Kriging Ordinario (OK) es menor que el de IDW (Inverse Distance Weighting), esto sugiere que el modelo de OK ofrece predicciones más precisas del Carbono Orgánico del Suelo (SOC) en Boyacá, aprovechando eficientemente la autocorrelación espacial presente en la región. Por el contrario, si el RMSE de IDW es más bajo, podría indicar que la distribución del SOC en Boyacá es más homogénea o que la densidad de muestreo no es suficiente para modelar adecuadamente la autocorrelación espacial.

El método de interpolación más confiable se selecciona en función del valor más bajo de RMSE, garantizando una representación precisa de la variabilidad del SOC en la región de Boyacá.

7. Guardado de resultados de la interpoblaciòn

En esta sección, se guardan los resultados de las interpolaciones realizadas para la región de Boyacá. Utilizando el paquete stars de R, se exportan las superficies estimadas de Carbono Orgánico del Suelo (SOC) generadas mediante los métodos de IDW (Inverse Distance Weighting) y Kriging Ordinario (OK) en formato GeoTIFF. Este formato es ideal para almacenar datos raster geoespaciales y es compatible con diversos software de Sistemas de Información Geográfica (SIG), como QGIS o ArcGIS, lo que facilita la realización de análisis adicionales y la creación de mapas de presentación.

# Guardar resultados de IDW para Boyaca
#write_stars(
# z1, dsn = "data/IDW_soc_boyaca.tif", layer = 1)
# Guardar resultados de Kriging Ordinario para Boyaca
#write_stars(
#  z2, dsn = "data/OK_soc_boyaca.tif", layer = 1)

Estos archivos contienen las superficies interpoladas de Carbono Orgánico del Suelo (SOC) en Boyacá, y pueden ser visualizados o analizados posteriormente en programas de Sistemas de Información Geográfica (SIG). Esto permite una interpretación más detallada y precisa de la distribución espacial del SOC en la región estudiada, facilitando su uso en análisis adicionales y en la generación de mapas temáticos.

8. Conclusiones

En este estudio se exploraron y aplicaron métodos de interpolación espacial para estimar la distribución del Carbono Orgánico del Suelo (SOC) en la región de Boyacá. Se emplearon dos enfoques principales: IDW (Inverse Distance Weighting) y Kriging Ordinario (OK), los cuales permitieron generar superficies continuas a partir de datos puntuales de SOC. Durante el proceso, se evidenció que la selección adecuada del modelo de variograma es fundamental para obtener resultados precisos en Kriging Ordinario, ya que este modelo define la autocorrelación espacial de los datos. Por su parte, el método IDW ofreció una alternativa más sencilla, aunque menos robusta, al no considerar explícitamente dicha autocorrelación.

Al comparar ambos métodos mediante validación cruzada, se identificaron diferencias en la precisión de las predicciones, destacando la importancia de utilizar métricas objetivas como el RMSE (Error Cuadrático Medio) para evaluar el desempeño de los modelos de interpolación. Si el RMSE para OK fue menor que el de IDW, esto sugiere que el modelo de Kriging Ordinario proporciona predicciones más precisas del SOC en Boyacá, aprovechando eficientemente la autocorrelación espacial. Por el contrario, si el RMSE de IDW fue más bajo, podría indicar que la distribución del SOC en la región es más homogénea o que la densidad de muestreo no fue suficiente para modelar adecuadamente la autocorrelación espacial.

Finalmente, el uso de R y sus paquetes especializados, como sf, stars, gstat y leaflet, facilitó el flujo de trabajo, desde la exploración de datos hasta la visualización y almacenamiento de resultados en formato GeoTIFF. Este formato, compatible con software de Sistemas de Información Geográfica (SIG) como QGIS o ArcGIS, permitió una gestión eficiente de los datos espaciales y una interpretación más detallada de la distribución del SOC en Boyacá.

Este ejercicio no solo fortaleció el conocimiento en técnicas de interpolación espacial, sino también en el análisis y manejo de datos geoespaciales aplicados al estudio del SOC. Además, se destacó la importancia de considerar factores como la densidad de muestreo y la autocorrelación espacial para mejorar la precisión de las estimaciones. Los resultados obtenidos, tanto en términos de mapas interactivos como de archivos raster, proporcionan una base sólida para futuros análisis y la toma de decisiones relacionadas con la gestión del suelo en la región de Boyacá.

9. Referencias

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

sessionInfo()
## R version 4.2.1 (2022-06-23 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## 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    
## 
## 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-1        stars_0.6-6        abind_1.4-5        sp_1.6-0          
##  [9] sf_1.0-16          terra_1.7-71      
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.12             lattice_0.20-45         FNN_1.1.4              
##  [4] leaflet.providers_2.0.0 png_0.1-8               class_7.3-21           
##  [7] zoo_1.8-11              digest_0.6.31           utf8_1.2.3             
## [10] R6_2.5.1                plyr_1.8.8              evaluate_0.20          
## [13] e1071_1.7-13            ggplot2_3.4.2           highr_0.10             
## [16] pillar_1.9.0            rlang_1.1.0             rstudioapi_0.14        
## [19] raster_3.6-26           jquerylib_0.1.4         rmarkdown_2.21         
## [22] rgdal_1.6-6             htmlwidgets_1.6.2       munsell_0.5.0          
## [25] proxy_0.4-27            compiler_4.2.1          xfun_0.38              
## [28] pkgconfig_2.0.3         base64enc_0.1-3         htmltools_0.5.4        
## [31] tidyselect_1.2.0        tibble_3.2.1            intervals_0.15.4       
## [34] codetools_0.2-19        reshape_0.8.9           fansi_1.0.4            
## [37] spacetime_1.3-1         dplyr_1.1.4             grid_4.2.1             
## [40] jsonlite_1.8.4          gtable_0.3.3            lifecycle_1.0.3        
## [43] DBI_1.1.3               magrittr_2.0.3          units_0.8-1            
## [46] scales_1.2.1            KernSmooth_2.23-20      cli_3.6.0              
## [49] cachem_1.0.7            farver_2.1.1            bslib_0.4.2            
## [52] ellipsis_0.3.2          xts_0.13.2              vctrs_0.6.5            
## [55] generics_0.1.3          tools_4.2.1             glue_1.6.2             
## [58] crosstalk_1.2.0         parallel_4.2.1          fastmap_1.1.1          
## [61] yaml_2.3.7              colorspace_2.1-0        classInt_0.4-9         
## [64] knitr_1.42              sass_0.4.5