Este informe aplica dos técnicas de interpolación espacial, IDW y Kriging Ordinario (OK), para estimar el contenido de carbono orgánico del suelo (SOC) en Nariño a 15-30 cm de profundidad. IDW asigna valores según la proximidad de las muestras, mientras que OK incorpora la estructura espacial de los datos para mayor precisión.
rm(list=ls())
# Borrar los mensajes no necesarios
options(showPackageStartupMessages = FALSE)
library(terra)
## terra 1.8.5
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)
Como se indicó en la introducción, las muestras corresponden a valores de SOC a 15-30 cm de profundidad y fueron obtenidas de SoilGrids 250m mediante este notebook. Estas muestras se guardaron en formato geopackage dentro del directorio de datos, utilizando el nombre de archivo soc_narino.gpkg.
Para verificar la disponibilidad de los datos se realiza lo siguiente:
list.files(path="datos", pattern = "*.gpkg")
## [1] "NARINO_.gpkg" "Nariño_Cities.gpkg" "soc_narino.gpkg"
# path (la carpeta), pattern (para que solo aparezcan los archivos gpkg)
En el siguiente código, se carga correctamente el archivo
soc_narino.gpkg de la carpeta Datos utilizando la
función st_read() del paquete sf.
samples <- sf::st_read("datos/soc_narino.gpkg")
## Reading layer `soc_narino' from data source
## `C:\Users\danna\Desktop\Geomática Básica 2024-2\Proyecto5\GB_sem_12\Datos\soc_narino.gpkg'
## using driver `GPKG'
## Simple feature collection with 1758 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -79.01206 ymin: 0.365009 xmax: -76.82632 ymax: 2.684845
## Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs
Dado que la área de estudio es Nariño, es importante leer también los datos correspondientes a los municipios de la región para contextualizar espacialmente los resultados. Esto se raliza a partir del siguiente código:
munic <- sf::st_read("datos/NARINO_.gpkg")
## Reading layer `municipios_colombia' from data source
## `C:\Users\danna\Desktop\Geomática Básica 2024-2\Proyecto5\GB_sem_12\Datos\NARINO_.gpkg'
## using driver `GPKG'
## Simple feature collection with 64 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -79.01021 ymin: 0.3613481 xmax: -76.83368 ymax: 2.683898
## Geodetic CRS: MAGNA-SIRGAS
Este archivo contiene información sobre los municipios dentro del área de estudio, lo que permitirá visualizar mejor la distribución espacial de los valores de SOC. Se trata de una colección de entidades geográficas con múltiples polígonos y atributos relevantes para el análisis. Su sistema de referencia es MAGNA-SIRGAS WGS 84, lo que asegura compatibilidad con otros datos geoespaciales utilizados en el estudio.
Para comprender mejor la distribución de los valores de SOC, realizamos un resumen estadístico de las muestras. Este resumen nos proporciona información sobre la variabilidad del contenido de SOC en la región, incluyendo valores mínimos, máximos, media y cuartiles. Además, graficamos un histograma para visualizar la distribución de los datos de una mejor manera con el siguiente código:
summary(samples)
## soc geom
## Min. : 18.24 POINT :1758
## 1st Qu.: 41.87 epsg:NA : 0
## Median : 61.58 +proj=long...: 0
## Mean : 78.12
## 3rd Qu.:103.89
## Max. :215.61
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 inicial nos permite identificar posibles valores atípicos
y tendencias en los datos antes de aplicar las técnicas de
interpolación.
Del histograma, se observa una frecuencia superior a 400 en los valores de SOC en el suelo, específicamente alrededor de 50 g/kg de SOC.
Para llevar a cabo la visualización de la distribución espacial del Carbono Orgánico del Suelo (SOC) en la región de Nariño, se utilizó un mapa interactivo generado con la librería Leaflet. Esta herramienta permite representar de forma clara y dinámica la ubicación de las muestras y su relación con las zonas geográficas de interés.
En esta parte, el código redondea los valores de la columna soc en el conjunto de datos samples a dos decimales. Esto es útil para evitar mostrar números con demasiada precisión en el mapa, facilitando su interpretación.
samples$soc = round(samples$soc, 2)
A continuación, se crea una paleta de colores que va desde el verde claro hasta el rojo intenso. Estos colores se asignan a los valores de soc (el contenido de carbono en el suelo) en los datos. Es decir, los valores más bajos de SOC estarán representados por un verde más claro, y los más altos por un rojo más intenso.
pal <- colorNumeric(
c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),
domain = samples$soc)
En el siguiente código, para generar el mapa, se cargaron las fronteras de los municipios utilizando la función addPolygons(), donde se definió un borde 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 restar visibilidad a los puntos de muestra.
A continuación, se añadieron los círculos de muestra utilizando la función addCircleMarkers().Cada círculo representa una muestra tomada, con un tamaño estándar de 1.5 píxeles, y se colorea de acuerdo con el valor de SOC de la muestra mediante una paleta de colores personalizada. Los círculos fueron configurados con opacidad completa, y se evitó el trazo del borde para resaltar el color de cada muestra.
Finalmente, se incluyó una leyenda en la esquina izquierda del mapa, usando la función addLegend(), que proporciona información sobre los valores de SOC asociados a cada color. Esto permite al usuario interpretar rápidamente los valores representados en el mapa.
El fondo del mapa fue provisto por OpenStreetMap mediante la función addProviderTiles(), ofreciendo un contexto geográfico adecuado para la ubicación de las muestras en la región.
leaflet() %>%
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 = F
) %>%
addLegend(
data = samples,
pal = pal,
values = ~soc,
position = "bottomleft",
title = "SOC:",
opacity = 0.7) %>%
addProviderTiles("OpenStreetMap")
En este mapa se observan los limites del departamento de Nariño y sus respectivos municipios, dandonos una idea de sus datos. Se puede ver que hay una mayor concentración de más de 150 de SOC en el área de los municipios de Olaya Herra, La Tota, El Charco, Barbacoas, Roberto Payán y Magui Payán.
Para estimar el contenido de carbono orgánico del suelo (SOC) en Nariño, se emplea la interpolación espacial mediante el método de Inverse Distance Weighting (IDW). Para ello, se utiliza la función gstat, que permite definir un modelo de interpolación con base 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 SOC (atributo que depende de si mismo (~1)) en ubicaciones desconocidas a partir de los datos de muestreo. La predicción se realiza 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 emplea la librería terra, utilizando las coordenadas que delimitan la zona de estudio en Nariño.
ext <- terra::ext(samples) # Para obtener la extensión real de los datos
rrr <- terra::rast(ext, resolution = 0.01, vals = 1, crs = "epsg:4326")
# vals: valores (le dice: "pongale 1 a todas las celdas")
rrr
## class : SpatRaster
## dimensions : 232, 219, 1 (nrow, ncol, nlyr)
## resolution : 0.009980555, 0.009999294 (x, y)
## extent : -79.01206, -76.82632, 0.365009, 2.684845 (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 es convertido en un objeto stars para realizar la interpolación, ya que solo lee este tipo de objeto. Con este procedimiento, se obtiene una superficie estimada de SOC en la región, permitiendo analizar su distribución espacial.
stars.rrr <- stars::st_as_stars(rrr)
Posteriormente, se ejecuta la interpolación con el modelo g1 y el raster stars.rrr mediante la función predict.
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 21.87028 53.45083 69.33817 81.59401 108.5943 207.9112 0
## var1.var NA NA NA NaN NA NA 50808
## dimension(s):
## from to offset delta refsys x/y
## x 1 219 -79.01 0.009981 WGS 84 [x]
## y 1 232 2.685 -0.009999 WGS 84 [y]
El objeto resultante z1 contiene dos atributos: var1.pred, que almacena los valores estimados de SOC, y var1.var, correspondiente a la varianza de predicción. Para manejar valores faltantes, estos se reemplazan por ceros:
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 21.87028 53.45083 69.33817 81.59401 108.5943 207.9112 0
## var1.var NA NA NA NaN NA NA 50808
## dimension(s):
## from to offset delta refsys x/y
## x 1 219 -79.01 0.009981 WGS 84 [x]
## y 1 232 2.685 -0.009999 WGS 84 [y]
De esta forma, se obtiene una superficie de SOC interpolada para Nariño, lista para su visualización y análisis.
names(z1) = "soc"
Para visualizar la interpolación de SOC en Nariño, se define una paleta de colores que representa los valores interpolados:
paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "#76EEC6"),
domain = range(z1$soc, na.rm = TRUE),
na.color = "transparent")
A continuación, se crea un mapa interactivo con leaflet, donde se superpone la capa interpolada sobre un mapa base. Se configura la opacidad de la interpolación y se añaden los puntos de muestreo, coloreados de acuerdo con sus valores de SOC:
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
Para realizar interpolación mediante Kriging Ordinario (OK) en Nariño, 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 es calcular el variograma empírico a partir de los puntos de muestreo. Esto se hace con la función variogram, donde se especifica la variable dependiente y los datos de calibración:
v_emp_ok = variogram(soc ~ 1, data=samples)
Posteriormente, se representó gráficamente el variograma empírico para analizar la estructura espacial de los datos:
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
)
Existen diversas maneras de ajustar un modelo a un variograma empírico. En este caso, se utilizó el método automático mediante la función autofitVariogram() del paquete automap, que selecciona el modelo más adecuado y ajusta sus parámetros de forma óptima.
v_mod_ok = automap::autofitVariogram(soc ~ 1, as(samples, "Spatial"))
Dado que autofitVariogram() no es compatible con objetos de tipo sf, fue necesario convertir samples en un SpatialPointsDataFrame utilizando el paquete sp.
Para visualizar el modelo de variograma ajustado, se utilizó la función plot():
plot(v_mod_ok)
v_mod_ok$var_model
Una vez obtenido el modelo de variograma ajustado, este se utilizó en la función gstat para llevar a cabo la interpolación mediante Kriging Ordinario. Este método estima los valores de la variable de interés en ubicaciones no muestreadas, considerando la dependencia espacial identificada en el variograma.
g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2= predict(g2, stars.rrr)
## [using ordinary kriging]
Después de realizar la predicción, se identificaron y reemplazaron los valores NA con 0.0 para evitar problemas en la visualización:
a <- which(is.na(z2[[1]]))
z2[[1]][a] = 0.0
names(z2) = "soc"
Ahora se genera el mapa correspondiente a el porcentaje de Interpolación OK y SOC. El cual es el siguiente:
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 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 # Mostrar el mapa
Para analizar las predicciones de carbono orgánico del suelo (SOC) en Nariño, se realizaron interpolaciones utilizando los métodos de IDW (Inverse Distance Weighting) y Kriging Ordinario (OK). Las predicciones se visualizaron en un mapa interactivo usando Leaflet con la siguiente paleta de colores:
Este enfoque facilita la comparación cualitativa de los patrones espaciales generados por ambos métodos de interpolación en la región de estudio.
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 ayuda a identificar patrones espaciales y posibles anomalías en la distribución del SOC, que podrían estar relacionadas con factores ambientales o la densidad de muestreo.
Para evaluar objetivamente la precisión de los métodos de interpolación, se utilizó la Validación Cruzada Leave-One-Out (LOOCV). Este enfoque consiste en:
El Error Cuadrático Medio (RMSE) se calculó para evaluar la precisión de las predicciones utilizando la fórmula:
\[RMSE = \sqrt{\frac{\sum{(predicho - observado)^2}}{n}}\]
### 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 para el método OK, se repitió el proceso para los resultados de IDW en Nariño. La comparación de ambos valores de RMSE permitió evaluar objetivamente la precisión de las predicciones:
Si el RMSE para OK es menor que el de IDW, esto sugiere que el modelo de Kriging Ordinario proporciona predicciones de SOC más precisas en Nariño, aprovechando la autocorrelación espacial en la región. Por otro lado, si IDW obtiene un RMSE más bajo, podría indicar que la distribución de SOC en Nariño es más homogénea o que la densidad de muestreo no es suficiente para modelar la autocorrelación espacial de manera efectiva.
El método de interpolación más confiable se selecciona en función del valor más bajo de RMSE, asegurando una representación precisa de la variabilidad del SOC en la región de Nariño.
En este apartado, se guardan los resultados de las interpolaciones realizadas para la región de Nariño. Utilizando el paquete stars de R, se exportan las superficies estimadas de Carbono Orgánico del Suelo (SOC) obtenidas mediante los métodos de IDW (Inverse Distance Weighting) y Kriging Ordinario (OK) en formato GeoTIFF.
Este formato permite almacenar datos raster geoespaciales y es compatible con diversos software de SIG (Sistemas de Información Geográfica) como QGIS o ArcGIS, lo cual facilita análisis adicionales y la generación de mapas de presentación.
# Guardar resultados de IDW para Nariño
#write_stars(
# z1, dsn = "data/IDW_soc_stder.tif", layer = 1)
# Guardar resultados de Kriging Ordinario para Nariño
#write_stars(
# z2, dsn = "data/OK_soc_stder.tif", layer = 1)
Estos archivos almacenan las superficies interpoladas de SOC en Nariño, y pueden ser visualizados o analizados posteriormente en programas de SIG, permitiendo una mejor interpretación de la distribución espacial del Carbono Orgánico del Suelo en la región estudiada.
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 Nariño. Se utilizaron dos enfoques: 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 comprobó que la selección adecuada del modelo de variograma es fundamental para obtener resultados precisos en OK, ya que define la autocorrelación espacial de los datos y, por ende, influye en la precisión de las predicciones. Por otro lado, el método IDW ofreció una alternativa más sencilla y directa, aunque menos robusta al no considerar dicha autocorrelación.
Sin embargo, no se logró ejecutar la validación cruzada para comparar ambos métodos debido a las limitaciones de capacidad del computador utilizado. Esta restricción impidió obtener resultados concluyentes sobre la precisión de las predicciones y, por lo tanto, no fue posible determinar cuál método ofrece un mejor desempeño en la estimación del SOC en la región. A pesar de esta limitación, el enfoque metodológico empleado sienta las bases para futuros análisis. Se recomienda utilizar un equipo con mayores recursos computacionales para poder realizar la validación cruzada y evaluar objetivamente la precisión de los modelos de interpolación.
Finalmente, el uso de R y sus paquetes especializados facilitó el flujo de trabajo, desde la exploración de datos hasta la visualización y almacenamiento de resultados en GeoTIFF, promoviendo una gestión eficiente de datos espaciales. Este ejercicio no solo fortaleció el conocimiento en interpolación espacial, sino también en el análisis y manejo de datos geoespaciales aplicados al estudio del SOC en Nariño.
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-8 abind_1.4-8 sp_2.1-4
## [9] sf_1.0-19 terra_1.8-5
##
## 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-30 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