Para el presente documento se usaron los datos del libro correspondiente a soil grid https://rpubs.com/churtadoj/Soilgrid

Descarga de libreria

library(terra)
## Warning: package 'terra' was built under R version 4.3.3
## terra 1.8.15
library(sf)
## Warning: package 'sf' was built under R version 4.3.3
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(sp)
## Warning: package 'sp' was built under R version 4.3.3
library(stars)
## Warning: package 'stars' was built under R version 4.3.3
## Loading required package: abind
## Warning: package 'abind' was built under R version 4.3.3
library(gstat)
## Warning: package 'gstat' was built under R version 4.3.3
library(automap)
## Warning: package 'automap' was built under R version 4.3.3
library(RColorBrewer)
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.3.3
library(leafem)
## Warning: package 'leafem' was built under R version 4.3.3
warning=FALSE
list.files(path="Datos", pattern = "*.gpkg")
## [1] "soc_Vlle.gpkg"
samples <- sf::st_read("Datos/soc_Vlle.gpkg")
## Reading layer `soc_Vlle' from data source 
##   `C:\Users\lacj7\OneDrive\Escritorio\GB2\Proyecto_fin\Datos\soc_Vlle.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1744 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -77.54262 ymin: 3.092656 xmax: -75.70974 ymax: 5.047424
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs
munic <- sf::st_read("C:/Users/lacj7/OneDrive/Escritorio/GB2/PROYECTO_1/DATOS/mun_Valle.gpkg")
## Reading layer `municipios_colombia' from data source 
##   `C:\Users\lacj7\OneDrive\Escritorio\GB2\PROYECTO_1\DATOS\mun_Valle.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 42 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -77.54977 ymin: 3.091239 xmax: -75.70724 ymax: 5.047394
## Geodetic CRS:  MAGNA-SIRGAS
summary(samples)
##       soc                    geom     
##  Min.   : 15.57   POINT        :1744  
##  1st Qu.: 37.75   epsg:NA      :   0  
##  Median : 61.89   +proj=long...:   0  
##  Mean   : 93.78                       
##  3rd Qu.:163.15                       
##  Max.   :269.87
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'
g1 = gstat(formula = soc ~ 1, data = samples)
(rrr <- terra::rast(
  xmin = -77.6, xmax = -75.5,   # Longitud (de oeste a este)
  ymin = 3.0, ymax = 5.1,       # Latitud (de sur a norte)
  nrows = 185, ncols = 165, 
  vals = 1, 
  crs = "epsg:4326"
))
## class       : SpatRaster 
## dimensions  : 185, 165, 1  (nrow, ncol, nlyr)
## resolution  : 0.01272727, 0.01135135  (x, y)
## extent      : -77.6, -75.5, 3, 5.1  (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  17.04321 53.49453 65.85637 94.97385 143.5604 255.9833     0
## var1.var         NA       NA       NA      NaN       NA       NA 30525
## dimension(s):
##   from  to offset    delta refsys x/y
## x    1 165  -77.6  0.01273 WGS 84 [x]
## y    1 185    5.1 -0.01135 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  17.04321 53.49453 65.85637 94.97385 143.5604 255.9833     0
## var1.var         NA       NA       NA      NaN       NA       NA 30525
## dimension(s):
##   from  to offset    delta refsys x/y
## x    1 165  -77.6  0.01273 WGS 84 [x]
## y    1 185    5.1 -0.01135 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 [%]"
    )
## 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)
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   308.6201   0.0000   0.0
## 2   Ste 44295.3691 782.6086   0.7
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 [%]"
    )
## 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
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 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

Realizamos validacion con la funcion gstat.cv

###Al realizar este cuaderno de interpolación espacial en R, aprendí a manejar datos geoespaciales desde su carga y exploración hasta su análisis e interpolación con gstat. Primero, cargué y revisé los datos, entendiendo su distribución mediante histogramas y resúmenes estadísticos. Luego, construí un modelo de variograma para estudiar la estructura espacial de los datos y ajusté un modelo de interpolación basado en Kriging. También descubrí cómo optimizar la resolución de los rásteres usando terra, algo clave para trabajar en equipos con poca memoria RAM. Sin embargo, tuve varios problemas porque mi computador tiene recursos limitados y una conexión a internet lenta. La carga de datos grandes y la creación de interpolaciones espaciales a veces tardaban demasiado o hacían que R se bloqueara. Para solucionarlo, reduje la resolución de los rásteres, usé menos puntos de muestreo y evité procesos que consumieran mucha memoria, como gráficos muy detallados o mapas interactivos pesados. Además, aprendí que, con una conexión lenta, es mejor descargar bibliotecas y datos con anticipación y trabajar con archivos locales para evitar demoras por el uso de memorias USB