El carbono orgánico del suelo (COS) es un indicador clave de la calidad y fertilidad del suelo, ya que influye en su capacidad para retener agua, mejorar la estructura y proporcionar nutrientes esenciales para los ecosistemas agrícolas y naturales. En la región de Magdalena, Colombia, la distribución del carbono orgánico en los suelos es un factor crucial para comprender los procesos de degradación, conservación y manejo sostenible de los recursos naturales.
Este estudio tiene como objetivo analizar la distribución espacial del carbono orgánico en los suelos del Magdalena mediante técnicas de interpolación geoestadística. A través del método de Kriging y la ponderación por distancia inversa (IDW), se busca estimar la variabilidad del COS en la región y evaluar la precisión de los modelos aplicados. Los resultados obtenidos permitirán identificar patrones espaciales del carbono orgánico, lo que contribuirá a la toma de decisiones en la planificación agrícola y ambiental.
Además, se realizará una validación cruzada de los modelos utilizados para estimar la exactitud de las predicciones y evaluar la magnitud de los errores residuales. De este modo, este estudio no solo aportará información relevante sobre la distribución del carbono orgánico en el suelo del Magdalena, sino que también servirá como base para futuras investigaciones en el manejo sostenible de los suelos en la regió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)
library(knitr)
##
## Adjuntando el paquete: 'knitr'
## The following object is masked from 'package:terra':
##
## spin
list.files(path = "Data", pattern = "*.gpkg")
## [1] "areas agua_Magdalena.gpkg" "Ciudade_ Magdalena.gpkg"
## [3] "Departamento_,Magdalena.gpkg" "Mun_Magdalena.gpkg"
## [5] "rios_Magdalena.gpkg" "soc_magdalena.gpkg"
## [7] "soc_Sucre.gpkg" "Vias_Magdalena.gpkg"
samples <- sf::st_read("data/soc_magdalena.gpkg")
## Reading layer `soc_magdalena' from data source
## `C:\Users\Usuario\Documents\Proyecto 4\Magdalena\Data\soc_magdalena.gpkg'
## using driver `GPKG'
## Simple feature collection with 1442 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -75.01995 ymin: 8.940989 xmax: -73.51486 ymax: 11.33822
## Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs
munic <- sf::st_read("data/Mun_Magdalena.gpkg")
## Reading layer `mgn_adm_mpio_grafico' from data source
## `C:\Users\Usuario\Documents\Proyecto 4\Magdalena\Data\Mun_Magdalena.gpkg'
## using driver `GPKG'
## Simple feature collection with 30 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.9466 ymin: 8.936489 xmax: -73.54184 ymax: 11.34891
## Geodetic CRS: MAGNA-SIRGAS
summary(samples)
## soc geom
## Min. : 8.514 POINT :1442
## 1st Qu.: 14.549 epsg:NA : 0
## Median : 18.465 +proj=long...: 0
## Mean : 24.026
## 3rd Qu.: 27.932
## Max. :110.200
El histograma titulado “Histograma of Samples” representa la distribución del carbono orgánico del suelo (SOC) en función de su frecuencia. La mayoría de las muestras presentan valores bajos de SOC, ya que las primeras barras del histograma tienen la mayor frecuencia, superando las 800 observaciones en la categoría más representativa. A medida que el contenido de SOC aumenta, la frecuencia disminuye progresivamente, lo que indica que los valores altos de carbono orgánico en el suelo son menos comunes. La distribución de los datos muestra un sesgo a la derecha, lo que significa que la mayoría de los valores están concentrados en rangos bajos, pero existen algunas muestras con valores más altos. Este patrón es común en suelos donde la acumulación de carbono orgánico depende de factores como la cobertura vegetal, el tipo de suelo y las prácticas agrícolas. Además, la presencia de algunas barras con baja frecuencia en valores más altos sugiere posibles valores atípicos o suelos con características especiales que requieren un análisis más detallado. En general, este análisis indica que la mayoría de los suelos en el conjunto de datos tienen contenidos relativamente bajos de carbono orgánico, con pocas muestras alcanzando valores significativamente altos. Esta información es útil para comprender la variabilidad en el contenido de SOC en la región estudiada y su relación con la fertilidad del suelo y su manejo.
samples$soc <- as.numeric(as.character(samples$soc))
hist(samples$soc)
samples$soc = round(samples$soc,2)
pal <- colorNumeric(
c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),
domain = samples$soc,
)
El mapa resultante permite visualizar la distribución espacial del carbono orgánico en distintos municipios. Las áreas delineadas en cian representan los municipios, mientras que los puntos de muestreo muestran los valores específicos de SOC. La escala de colores ayuda a identificar patrones, como zonas con alta o baja concentración de carbono orgánico en el suelo. Si los puntos con valores altos están agrupados en ciertas áreas, podría indicar regiones con mayor materia orgánica, posiblemente debido a condiciones climáticas, tipo de vegetación o manejo del suelo. En cambio, si hay una dispersión amplia de valores, podría reflejar una variabilidad en los suelos de la región. En general, este mapa es una herramienta útil para analizar la distribución espacial del SOC y facilitar la toma de decisiones en estudios de suelos y conservación ambiental.
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.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'
g1 = gstat(formula = soc ~ 1, data = samples)
(rrr <- terra::rast(xmin=-75.01995, xmax=-73.54184, ymin=8.936489, ymax=11.34891, nrows=370, ncols = 329, vals = 1, crs = "epsg:4326"))
## class : SpatRaster
## dimensions : 370, 329, 1 (nrow, ncol, nlyr)
## resolution : 0.004492736, 0.006520057 (x, y)
## extent : -75.01995, -73.54184, 8.936489, 11.34891 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : lyr.1
## min value : 1
## max value : 1
stars.rrr <- stars::st_as_stars(rrr)
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 9.066499 17.43669 21.53381 24.52109 29.25125 106.7481 0
## var1.var NA NA NA NaN NA NA 121730
## dimension(s):
## from to offset delta refsys x/y
## x 1 329 -75.02 0.004493 WGS 84 [x]
## y 1 370 11.35 -0.00652 WGS 84 [y]
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 9.066499 17.43669 21.53381 24.52109 29.25125 106.7481 0
## var1.var NA NA NA NaN NA NA 121730
## dimension(s):
## from to offset delta refsys x/y
## x 1 329 -75.02 0.004493 WGS 84 [x]
## y 1 370 11.35 -0.00652 WGS 84 [y]
names(z1) = "soc"
paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), domain = 10:120, na.color = "transparent")
El mapa generado con leaflet y la interpolación IDW permite visualizar la variabilidad del carbono orgánico en el suelo a nivel espacial. La combinación de puntos de muestreo con una capa interpolada ayuda a identificar zonas con valores altos o bajos de SOC y facilita la toma de decisiones en estudios ambientales y de conservación del suelo.
Las zonas con tonos verdes indican valores más altos de carbono orgánico en el suelo. Las áreas con tonos naranja o amarillo representan valores más bajos de SOC. La interpolación IDW estima valores en áreas no muestreadas en función de la cercanía a los puntos medidos.
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
v_emp_ok = variogram(soc ~ 1, data=samples)
La curva muestra un crecimiento de la semivarianza con la distancia hasta alcanzar el umbral en torno a los 34 unidades de distancia. Esto indica que la variabilidad espacial del SOC tiene una influencia hasta esa distancia, y más allá de ese punto, los valores de SOC son independientes entre sí.
En términos prácticos, este variograma es útil para modelar la estructura espacial del SOC y puede utilizarse en técnicas de interpolación como el kriging, que ayuda a estimar valores en áreas no muestreadas con base en la dependencia espacial observada en los datos.
plot(v_emp_ok)
v_mod_ok = automap::autofitVariogram(soc ~ 1, as(samples, "Spatial"))
plot(v_mod_ok)
v_mod_ok$var_model
## model psill range kappa
## 1 Nug 26.40112 0.00000 0.0
## 2 Ste 173.58404 49.28774 0.6
g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2= predict(g2, stars.rrr)
## [using ordinary kriging]
a <- which(is.na(z2[[1]]))
z2[[1]][a] = 0.0
names(z2) = "soc"
Este mapa muestra la distribución espacial del SOC en la región de interés, interpolada mediante kriging ordinario.
Las áreas en naranja y amarillo indican niveles más bajos de SOC. Las zonas en cian y verde representan concentraciones más altas de SOC. Los valores de SOC en puntos de muestreo ayudan a validar la interpolación. El uso del kriging en este contexto permite generar una superficie continua de estimaciones basadas en los datos muestreados, optimizando la predicción en áreas donde no se tienen mediciones directas.
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 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")
Este mapa permite comparar visualmente los resultados de la interpolación del SOC utilizando dos métodos diferentes:
IDW (Inverse Distance Weighting): Método basado en la proximidad, donde los puntos más cercanos influyen más en la estimación. Puede generar patrones más suaves y con menos variabilidad. OK (Ordinary Kriging): Método basado en modelos geoestadísticos que considera la correlación espacial entre puntos. Puede producir estimaciones más precisas en áreas con datos dispersos. Naranja y amarillo → Indican valores bajos de SOC. Cian y verde → Representan valores altos de SOC.
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
El gráfico presentado muestra los residuales del modelo de interpolación utilizado para estimar los valores de carbono orgánico en el suelo. En este gráfico, cada círculo representa un punto de muestreo donde se ha calculado la diferencia entre el valor predicho por el modelo y el valor observado en campo.
El tamaño de los círculos indica la magnitud del error: los círculos más grandes representan errores más altos, mientras que los más pequeños corresponden a diferencias menores entre la predicción y el valor real. Por otro lado, el color de los círculos refleja la dirección del error: los círculos verdes indican que el modelo subestimó los valores observados, mientras que los círculos rosados representan puntos donde el modelo sobreestimó los valores reales.
La escala de la derecha muestra los valores numéricos de los residuales. Por ejemplo, un valor de -39.383 indica que en ese punto el modelo predijo un valor mucho más alto que el observado, mientras que un residual de 54.628 señala que el modelo subestimó significativamente el valor real. Valores cercanos a cero (como -0.52) sugieren que en esas ubicaciones la predicción del modelo fue bastante precisa.
include_graphics("C:/Users/Usuario/Documents/Proyecto 4/Magdalena/Data/Residual.png.")
un RMSE de 6.95653 significa que, en promedio, las predicciones del modelo presentan una desviación de aproximadamente 6.96 unidades en la misma escala de los datos. La magnitud de este error debe ser analizada en función del rango de la variable estudiada. Si los valores observados oscilan dentro de un rango amplio, como de 0 a 100, este error podría considerarse aceptable. Sin embargo, si los valores se encuentran dentro de un rango reducido, por ejemplo, de 0 a 10, un error de 6.96 indicaría una baja precisión del modelo.
write_stars(
z1, dsn = "data/IDW_soc_magdalena.tif", layer = 1
)
write_stars(
z2, dsn = "data/OK_soc_magdalena.tif", layer = 1
)
sessionInfo()
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## 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
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] knitr_1.49 leafem_0.2.3 leaflet_2.2.2 RColorBrewer_1.1-3
## [5] automap_1.1-16 gstat_2.1-3 stars_0.6-8 abind_1.4-8
## [9] sp_2.1-4 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-13
## [49] png_0.1-8 evaluate_1.0.1 rlang_1.1.4
## [52] Rcpp_1.0.13-1 glue_1.8.0 DBI_1.2.3
## [55] rstudioapi_0.17.1 reshape_0.8.9 jsonlite_1.8.9
## [58] R6_2.5.1 plyr_1.8.9 intervals_0.15.5
## [61] units_0.8-5