1. Introduccion

Este cuaderno ilustra dos interpolaciones espaciales. IDW es una tecnica arbitraria y Ok es una probabilistica. Ambas van a ser utilizadas para obtener visualizaciones del SOC.

2. Configuracion

Limpiamos el espacio de trabajo:

rm(list=ls())

Para la lectura y lectura de datos, vamos a usar terra para datos raster y sf para vectoriales.

Primero, vamos a instalar en las consolas las librerias necesarias, esto desde la consola y solo instalando las librerias que hagan falta.

En el siguiente cuadro subimos las librerias pero primero comprobamos la existencia y correcta instalacion de cada una.

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

Lectura de datos de entrada

Para estepaso vamos a tener en cuenta los valores SOC que extrajimos de cuaderno de soilgrids, ahora vamos aplicar un cuadro de codigo que filtre los archivos que tenemos guardados en gpkg.

getwd()
## [1] "C:/Users/DELL/Documents"
list.files(path="Datos", pattern = "*.gpkg")
## [1] "Amazonas Munic.gpkg" "soc_amaz.gpkg"

Ahora vamos a leer los documentos usando sf, tanto el soc anteriormente guardado del cuaderno soilgrids como los datos de nuestro departamente. En este caso el amazonas.

samples <- sf::st_read("Datos/soc_amaz.gpkg")
## Reading layer `soc_amaz' from data source 
##   `C:\Users\DELL\Documents\Datos\soc_amaz.gpkg' using driver `GPKG'
## Simple feature collection with 1953 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.39182 ymin: -4.22843 xmax: -69.3691 ymax: 0.1538592
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs
munic <- sf::st_read("Datos/Amazonas Munic.gpkg")
## Reading layer `cortado' from data source 
##   `C:\Users\DELL\Documents\Datos\Amazonas Munic.gpkg' using driver `GPKG'
## Simple feature collection with 8 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.38985 ymin: -4.228429 xmax: -69.36835 ymax: 0.1603
## Geodetic CRS:  WGS 84

4. Exploracion de los datos de entrada.

Vamos a pedir a R un resumen de los datos de entrada que agregamos.

summary(samples)
##       soc                    geom     
##  Min.   : 10.29   POINT        :1953  
##  1st Qu.: 20.95   epsg:NA      :   0  
##  Median : 30.95   +proj=long...:   0  
##  Mean   : 54.27                       
##  3rd Qu.: 59.16                       
##  Max.   :477.26

A razon de estos datos vamos a proyectar un histograma:

terra::hist(samples$soc,
  col = topo.colors(5),  # Como paleta elegimos esta multicolor
  main = "Distribución estándar del carbono orgánico del suelo",
  xlab = "(SOC)",
  ylab = "Frecuencia",
  border = "white"
)

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

Ahora vamos a plotear nuestros datos para entender como se distribuyen de forma geoespacial:

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")

5. Interpolacion

5.1 Creacion del objeto gstat

Para hacer la interpolacion necesitamos un objeto gstat, este objeto que obtiene la informacion necesaria para hacer la grafica espacial.

5.2 Interpolacion

Vamos a crear del objeto gstat

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

Vamos a utilizar la libreria terra para crear un objeto raster que cubra el area de interes.

(rrr <- terra::rast(xmin=-74.55524, xmax=-72.38985, ymin=5.708308, ymax=8.145474, nrows=370, ncols = 329, vals = 1, crs = "epsg:4326"))            
## class       : SpatRaster 
## size        : 370, 329, 1  (nrow, ncol, nlyr)
## resolution  : 0.006581733, 0.006586935  (x, y)
## extent      : -74.55524, -72.38985, 5.708308, 8.145474  (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 codigo cambia el spatraster a un objeto stars

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  45.35854 46.72108 47.23172 47.17207 47.65948 48.57154      0
## var1.var         NA       NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -74.56  0.006582 WGS 84 [x]
## y    1 370  8.145 -0.006587 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  45.35854 46.72108 47.23172 47.17207 47.65948 48.57154      0
## var1.var         NA       NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -74.56  0.006582 WGS 84 [x]
## y    1 370  8.145 -0.006587 WGS 84 [y]
names(z1) = "soc"
names(z1)
## [1] "soc" NA
# Quitar el atributo NA
z1 <- z1["soc"]
print(names(z1))
## [1] "soc"

Este cuadro corrige la ubicacion del bounding box para que cubra el amazonas ya que en un intento anterior este recuadro aparecia hacia el norte del pais.

bb <- st_bbox(samples)


dims <- st_dimensions(z1)
nx <- dims$x$to
ny <- dims$y$to


dims$x$offset <- bb["xmin"]
dims$x$delta  <- (bb["xmax"] - bb["xmin"]) / nx

dims$y$offset <- bb["ymax"]
dims$y$delta  <- (bb["ymin"] - bb["ymax"]) / ny   

# aplicar
st_dimensions(z1) <- dims

# verificar
print(st_dimensions(z1))
##   from  to offset.xmin delta.xmax refsys x/y
## x    1 329      -74.39    0.01527 WGS 84 [x]
## y    1 370      0.1539   -0.01184 WGS 84 [y]
plot(z1, reset = FALSE)
plot(st_geometry(samples), add = TRUE, pch = 20, col = "red")

Este ploteo rojo cumple la funcion de confirmar si esto confuerda con la informacion que necesitamos clasificar, y ahora si viene el ploteo real dle IDW.

paleta <- colorNumeric(
  palette = c("orange","yellow","cyan","green"),
  domain = z1$soc,
  na.color = "transparent"
)

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,
      colorOptions = colorOptions(
        palette = c("orange","yellow","cyan","green"),
        domain = range(z1$soc, na.rm = TRUE)
      )
  ) %>%
   addCircleMarkers(
    data = samples,
    radius = 1.5,
    label = ~soc,
    color = ~paleta(soc),
    fillOpacity = 1,
    stroke = FALSE
  ) %>%
  addLegend(
    "bottomright",
    pal = paleta,
    values = z1$soc,
    title = "IDW 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
m

5.3 Interpolacion OK

Para realizar este ploteo vamos a realizar unos variogramas previos que den la informacion que se quiere graficar en datos cuantitativos.

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

Ploteamos el variograma:

plot(v_emp_ok)

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

Ahora estos resultados son una lista de objetos.

v_mod_ok$var_model
##   model      psill    range kappa
## 1   Nug   468.0428    0.000   0.0
## 2   Ste 71431.6737 5020.363   0.7
g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2= predict(g2, stars.rrr)
## [using ordinary kriging]

Ahora se renombran los atributos

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

Las predicciones de estos datos se proyectan en el siguiente ploteo:

bb <- st_bbox(samples)
st_dimensions(z2)$x$offset <- bb["xmin"]
st_dimensions(z2)$x$delta  <- (bb["xmax"] - bb["xmin"]) / dim(z2)[1]

st_dimensions(z2)$y$offset <- bb["ymax"]
st_dimensions(z2)$y$delta  <- (bb["ymin"] - bb["ymax"]) / dim(z2)[2]
st_crs(z2) <- 4326
m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 0.7,                
      colorOptions = colorOptions(
        palette = c("orange", "yellow", "cyan", "green"), 
        domain = samples$soc)
    ) %>%
   addCircleMarkers(
    data = samples,
    radius = 1.5, 
    label = ~soc,
    color = ~paleta(soc),
    fillOpacity = 1,
    stroke = FALSE
  ) %>%
    addLegend(
      "bottomright", 
      pal = paleta, 
      values = samples$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
## 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

6. Asignacion de resultados

6.1 Asignacion cualitativa

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 [%]"
)
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m

6.2 Validacion cruzada

Esta validacion nos permite con los dos metodos IDW y OK obtener examinar los datos que nos permiten una vision del SOC

7. Guardar la interpolacion.

write_stars(
  z1, dsn = "Datos/IDW_soc_amaz.tif", layer = 1
)
write_stars(
  z2, dsn = "Datos/OK_soc_amaz.tif", layer = 1
)

8. Conclusiones

En este cuaderno de R se hace una interpolacion en dos modelos, en cuanto aprendizaje se comprende porque estos datos juntos dan una visualiacion mas precisa y en cuanto a aplicaciones porque han sido importante los anteriores cuadernos, sin estos seria mas dificil la solucion de problemas coo sucede en el ploteo de IDW de este cuaderno, sin la comprension del error el plot no hubiese sido correcto. de igual forma se interioriza la lectura de graficos y la profundizacion sobre el departamento.

9. Refencia

sessionInfo()
## R version 4.5.1 (2025-06-13 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] leafem_0.2.5       leaflet_2.2.3      RColorBrewer_1.1-3 automap_1.1-20    
##  [5] gstat_2.1-4        stars_0.6-8        abind_1.4-8        sp_2.2-0          
##  [9] sf_1.0-22          terra_1.8-80      
## 
## loaded via a namespace (and not attached):
##  [1] generics_0.1.4          sass_0.4.10             class_7.3-23           
##  [4] KernSmooth_2.23-26      lattice_0.22-7          digest_0.6.37          
##  [7] magrittr_2.0.4          evaluate_1.0.5          grid_4.5.1             
## [10] fastmap_1.2.0           plyr_1.8.9              jsonlite_2.0.0         
## [13] e1071_1.7-16            reshape_0.8.10          DBI_1.2.3              
## [16] crosstalk_1.2.2         scales_1.4.0            codetools_0.2-20       
## [19] jquerylib_0.1.4         cli_3.6.5               rlang_1.1.6            
## [22] units_1.0-0             intervals_0.15.5        base64enc_0.1-3        
## [25] cachem_1.1.0            yaml_2.3.10             FNN_1.1.4.1            
## [28] raster_3.6-32           tools_4.5.1             parallel_4.5.1         
## [31] dplyr_1.1.4             ggplot2_4.0.0           spacetime_1.3-3        
## [34] 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        
## [40] classInt_0.4-11         leaflet.providers_2.0.0 htmlwidgets_1.6.4      
## [43] pkgconfig_2.0.3         pillar_1.11.1           bslib_0.9.0            
## [46] gtable_0.3.6            glue_1.8.0              Rcpp_1.1.0             
## [49] tidyselect_1.2.1        tibble_3.3.0            xfun_0.53              
## [52] rstudioapi_0.17.1       knitr_1.50              farver_2.1.2           
## [55] htmltools_0.5.8.1       rmarkdown_2.30          xts_0.14.1             
## [58] compiler_4.5.1          S7_0.2.0