Interpolación Espacial de Bolívar

1.INTRODUCCION

En este notebook se realizará la ilustración de dos técnicas de interpolación espacial, una de ellas es la distancia ponderada inversa (IDW) y Kriging ordinario (OK) y la otra es la de IDW es una técnica determinista. El objetivo de este notebook es poder emplear las técnicas mencionadas anteriormente y lograr una obtención de una superficie continua de carbono orgánico del suelo (COS) a 15-30 cm a partir de muestras obtenidas de SoilGrids 250 m. El análisis de la distribución espacial del carbono orgánico del suelo (SOC) es fundamental para que podamos comprender la variabilidad en la calidad del suelo y su relación con la productividad agrícola y la conservación ambiental.

2.CONFIGURACIÓN

rm(list=ls())

library(terra)
library(sf)
library(sp)
library(stars)
library(gstat)
library(automap)
library(RColorBrewer)
library(leaflet)
library(leafem)

3.LECTURA DE DATOS

list.files(path="DATA", pattern = "*gpk*")

## [1] "AREAS.gpkg"         "cities2.gpkg"       "DEPTO_BOLIVAR.gpkg"
## [4] "LINEAS.gpkg"        "mpios_bolivar.gpkg" "Mun_Bolivar.gpkg"  
## [7] "soc_Bolivar.gpkg"   "soc_stder.gpkg"     "VIAS.gpkg"

samples <- sf::st_read("data/soc_Bolivar.gpkg")

## Reading layer `soc_Bolivar' from data source 
##   `C:\Users\Tecnologia\Desktop\GB2\PROYECTO_4\Bolivar\DATA\soc_Bolivar.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1533 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -76.31233 ymin: 7.001934 xmax: -73.57113 ymax: 10.79678
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

Este código carga un conjunto de datos geoespaciales en GeoPackage con 1533 puntos y con un solo campo de datos adicionales. Donde cada uno tiene coordenadas en longitud y latitud, y los puntos están distribuidos en un área delimitada por la caja de coordenadas que va desde la longitud -76.31 hasta -73.57 y desde la latitud 7.00 hasta 10.80.

munic <- sf::st_read("data/Mun_Bolivar.gpkg")

## Reading layer `mapa_colombia' from data source 
##   `C:\Users\Tecnologia\Desktop\GB2\PROYECTO_4\Bolivar\DATA\Mun_Bolivar.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 46 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.19063 ymin: 6.99916 xmax: -73.74578 ymax: 10.80147
## Geodetic CRS:  MAGNA-SIRGAS

4.LECTURA DE DATOS

summary(samples)

##       soc                     geom     
##  Min.   :  8.425   POINT        :1533  
##  1st Qu.: 16.674   epsg:NA      :   0  
##  Median : 24.223   +proj=long...:   0  
##  Mean   : 29.158                       
##  3rd Qu.: 36.613                       
##  Max.   :131.358

samples$soc <- as.numeric(as.character(samples$soc))
hist(samples$soc)

Se puede observar que la mayoría de los valores se encuentran concentrados en un rango bajo, en especial los que están entre 0 y 40, por lo cual, hay una distribución sesgada a la derecha. La frecuencia es alta en los primeros intervalos, permitiendo que alcance su punto máximo alrededor de 20.

A medida que los valores aumentan, la frecuencia disminuye de manera progresiva, sugiriendo la presencia de algunos valores atípicos en el extremo superior (mayores a 100), aunque estos son poco frecuentes.

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"),
  # colors depend on the count variable
  domain = samples$soc,
  )

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 = 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'

Se puede evidenciar que la mayoría de los valores bajos que están en color amarillo, se concentran en la parte norte y central del departamento, en especial en áreas cercanas a Cartagena y los municipios de la Depresión Momposina. Por otro lado, observamos un contraste los valores más altos que están en los colores naranja y rojo, estos se encuentran en la zona sur y suroeste.

5.INTERPOLACIÓN

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

(rrr <- terra::rast(xmin=-76.19063, xmax=-73.74578, ymin=6.99916, ymax=10.80147, nrows=370, ncols = 329, vals = 1, crs = "epsg:4326"))

## class       : SpatRaster 
## dimensions  : 370, 329, 1  (nrow, ncol, nlyr)
## resolution  : 0.007431155, 0.01027651  (x, y)
## extent      : -76.19063, -73.74578, 6.99916, 10.80147  (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  8.70542 21.4968 26.66365 28.65421 33.17702 130.2811      0
## var1.var        NA      NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset    delta refsys x/y
## x    1 329 -76.19 0.007431 WGS 84 [x]
## y    1 370   10.8 -0.01028 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  8.70542 21.4968 26.66365 28.65421 33.17702 130.2811      0
## var1.var        NA      NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset    delta refsys x/y
## x    1 329 -76.19 0.007431 WGS 84 [x]
## y    1 370   10.8 -0.01028 WGS 84 [y]

names(z1) = "soc"

paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), domain = 10:120, na.color = "transparent")

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 [%]"
    )

m  # Print the map

Este mapa nos muestra la forma en la que se distribuye el carbono orgánico del suelo en Bolívar. Los valores más bajos que son (amarillo y naranja) están en el norte y centro, cerca de Cartagena, una posibilidad es porque sus suelos son más pobres y tienen mayor actividad humana. Los valores más altos que son (verde y azul) aparecen en el sur y en otras zonas, donde el suelo puede ser más fértil y menos afectado por la agricultura. Esto nos da a entender que el contenido de carbono en el suelo varía según el uso del terreno y las condiciones naturales.

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

plot(v_emp_ok)

En este análisis podemos entender de qué manera se distribuye el carbono orgánico en el suelo. A través del método IDW, identificamos patrones de variabilidad en los valores, mostrando zonas con mayor y menor concentración. Por otro lado, podemos observar que el variograma muestra que la semivarianza aumenta con la distancia, lo que significa que los valores del suelo tienden a ser más diferentes entre más alejados están entre ellos.

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

plot(v_mod_ok)

Allí logramos evidenciar los valores de carbono orgánico en el suelo son parecidos cuando los puntos están cerca uno del otro. Sin embargo, a medida que aumenta la distancia entre ellos, las diferencias también crecen hasta un punto en el que ya no tienen relación.

En este caso, esa distancia es de 79 metros, lo que significa que si tomamos muestras más alejadas, los valores de carbono ya no estarán conectados entre sí, lo cual es útil para que decidamos cómo tomar muestras en el suelo y poder entender mejor su distribución.

v_mod_ok$var_model

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"

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 [%]"
    )

m  # Print the map

En este mapa podemos observar que en su gran mayoría, en especial en la zona norte y centro del departamento, hay una presencia de niveles bajos de SOC, lo cual puede estar relacionado con la actividad agrícola intensiva, la deforestación o la composición del suelo. Mientras que en la parte del sur y se suroeste, logramos evidenciar áreas con mayores niveles de carbono orgánico, lo cual indica que existe una mejor calidad del suelo, en este caso sería muy posible que se deba a la presencia de vegetación un poco densa o suelos menos degradados.

6.LECTURA DE LOS DATOS

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 [%]"
)

m  # Print the map

En este mapa se observa la manera en la que se distribuye el carbono orgánico en el suelo (SOC). En este caso, empleamos dos métodos para analizar y estimar los valores en áreas sin datos, los cuales son: IDW e OK. Allí se evidencian los colores, donde indican la cantidad de SOC: amarillo y naranja, estos representan valores bajos a medios, mientras que los tonos verdes muestran las zonas con más carbono en el suelo. La mayoría del territorio tiene valores medios, pero hay algunas áreas con más carbono, sobre todo en el sur y suroeste.

write_stars(
  z1, dsn = "DATA/IDW_soc_Bolivar.tif", layer = 1
)

write_stars(
  z2, dsn = "DATA/IDW_soc_Bolivar.tif", layer = 1
)

CONCLUSIONES

En pirmer lugar, en este notebook se analizó la distribución del carbono orgánico del suelo (SOC) en este caso, en el departamento de Bolívar, utilizando los métodos de interpolación IDW y Kriging Ordinario. Se encontró que los valores de SOC varían entre aproximadamente 20% y 120%, con las concentraciones más altas ubicadas en el sur y suroeste del área de estudio, mientras que las más bajas se presentan en el centro y norte. La interpolación mostró diferencias en la estimación de valores de acuerdo al método realizado, por un lado, con IDW resaltando transiciones más marcadas y por otro, el Kriging mostrando una distribución más suavizada.

En segundo lugar,los lenguajes que se usaron con paquetes como (sf, gstat y leaflet), permitió cargar datos, transformarlos a coordenadas geográficas correctas y la generación de mapas interactivos. En otra instancia se observó que el método IDW genera transiciones más marcadas, mientras que Kriging ofrece una distribución más suavizada.

Por último, se dedujo que estos resultados nos indican que el SOC no es homogéneo y está influenciado por factores ambientales y de manejo del suelo. Obteniendo como resultado que estos métodos son herramientas útiles para la planificación agrícola y ambiental, los cuales pueden adaptarse para analizar otros tipos de datos ambientales y mejorar la toma de decisiones en gestión del suelo

sessionInfo()

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

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.