1. Introducción

Este cuaderno explica brevemente como descargar datos raster de SoilGrigds usando R. Este proceso se lleva a cabo con datos de carbono organico en el suelo como una variable de interes, claramente el codigo se puede replicar usando una variable distinta.

Posterior a la descarga de los datos procederemos a ilustrar dos metodos de interpolacion espacial: Inverse Dinstance Weighted (IDW) y Ordinary Kriging (OK). IDW es una tecnica deterministica y OK es una tecnica probabilistica. Ambas técnicas se utilizan aquí para obtener una superficie continua de carbono orgánico del suelo (SOC) a una profundidad de 15-30 cm a partir de muestras obtenidas de SoilGrids 250 m.

2. Configuración inicial

limpiar el workspace

rm(list=ls())

Cargamos las librerias necesarias

Fijamos las variables de interes

Definimos el link que corresponde al sitio web de ISRIC

url = "https://files.isric.org/soilgrids/latest/data/"

Entonces, creamos algunos objetos indicando que necesitamos de ISRIC

#strings basicos
voi = "soc" #Soil organic carbon
depth = "15-30cm"
quantile = "mean"
#concatenar los strings
(variable = paste(url, voi, sep = ""))
## [1] "https://files.isric.org/soilgrids/latest/data/soc"

Ahora definimos la propiedad del suelo que queremos descargar

(layer = paste(variable,depth,quantile, sep = "_")) #capa de interes
## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean"

asi, la verificación VRT esta completa

(vrt_layer = paste(layer, '.vrt', sep = "" ))
## [1] "https://files.isric.org/soilgrids/latest/data/soc_15-30cm_mean.vrt"

3. Descargar TIFF de la region de interes

Lee un mapa vector con el area de interes previamente descargado

#Cambia la ruta y el nombre del archivo para que coincidan con tus datos.
(stder = st_read("C:/Users/vocav/Documents/GB2-2025 nicolin-20250723T070234Z-1-001/GB2-2025 nicolin/GB2-2025/GB2/PR4/nariñito.gpkg"))
## Reading layer `municipios' from data source 
##   `C:\Users\vocav\Documents\GB2-2025 nicolin-20250723T070234Z-1-001\GB2-2025 nicolin\GB2-2025\GB2\PR4\nariñito.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 64 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -79.01021 ymin: 0.3613481 xmax: -76.83368 ymax: 2.683898
## Geodetic CRS:  MAGNA-SIRGAS
## Simple feature collection with 64 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -79.01021 ymin: 0.3613481 xmax: -76.83368 ymax: 2.683898
## Geodetic CRS:  MAGNA-SIRGAS
## First 10 features:
##    dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1          52        001      52001     NARIÑO      PASTO
## 2          52        019      52019     NARIÑO      ALBÁN
## 3          52        022      52022     NARIÑO     ALDANA
## 4          52        036      52036     NARIÑO     ANCUYA
## 5          52        051      52051     NARIÑO   ARBOLEDA
## 6          52        079      52079     NARIÑO  BARBACOAS
## 7          52        083      52083     NARIÑO      BELÉN
## 8          52        110      52110     NARIÑO    BUESACO
## 9          52        203      52203     NARIÑO      COLÓN
## 10         52        207      52207     NARIÑO    CONSACÁ
##                           mpio_crslc mpio_tipo mpio_narea mpio_nano shape_Leng
## 1                               1540 MUNICIPIO 1099.37825      2023  1.8748401
## 2               Ordenanza 41 de 1903 MUNICIPIO   38.61009      2023  0.3402693
## 3               Ordenanza 11 de 1911 MUNICIPIO   47.64055      2023  0.4520361
## 4                               1544 MUNICIPIO   68.90636      2023  0.3736985
## 5                               1859 MUNICIPIO   60.18576      2023  0.3754062
## 6                               1916 MUNICIPIO 2732.93328      2023  3.6991487
## 7  Ordenanza 53 Noviembre 29 de 1985 MUNICIPIO   41.84541      2023  0.3732840
## 8                               1899 MUNICIPIO  635.96083      2023  1.2292312
## 9               Ordenanza 37 de 1921 MUNICIPIO   61.75053      2023  0.4592866
## 10                              1900 MUNICIPIO  119.35364      2023  0.4650471
##     shape_Area                           geom
## 1  0.089062266 MULTIPOLYGON (((-77.31099 1...
## 2  0.003129126 MULTIPOLYGON (((-77.06531 1...
## 3  0.003855327 MULTIPOLYGON (((-77.69156 0...
## 4  0.005578855 MULTIPOLYGON (((-77.50924 1...
## 5  0.004877181 MULTIPOLYGON (((-77.13479 1...
## 6  0.220963473 MULTIPOLYGON (((-77.79456 1...
## 7  0.003391678 MULTIPOLYGON (((-77.07227 1...
## 8  0.051533090 MULTIPOLYGON (((-77.23516 1...
## 9  0.005005108 MULTIPOLYGON (((-77.04473 1...
## 10 0.009664908 MULTIPOLYGON (((-77.50158 1...

Nota el CRS de stder. Como las capas de ISRIC usan la proyección Homolosine, necesitamos reproyectar nuestra capa usando la biblioteca sf.

repro = '+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
stder_reproj = st_transform(stder, repro)

Ahora definimos los limites de nuestra area de interes

(box = st_bbox(stder_reproj))
##        xmin        ymin        xmax        ymax 
## -8796326.48    40225.08 -8554103.93   298770.17

Ahora, usemos los datos de bbox para definir los límites de la caja de delimitación tal como los usa la biblioteca GDAL. Por cierto, esta es una de las partes más complicadas de usar GDAL.

ulx = box$xmin
uly = box$ymax
lrx = box$xmax
lry = box$ymin
(bbx = c(ulx, uly, lrx,lry))
##        xmin        ymax        xmax        ymin 
## -8796326.48   298770.17 -8554103.93    40225.08

Ahora, podemos usar la función gdal_translate para descargar la capa VRT. Primero, definamos dónde guardar los datos:

sg_url="/vsicurl/https://files.isric.org/soilgrids/latest/data/"
datos = 'soc/soc_15-30cm_mean.vrt'
file = "C:/Users/vocav/Documents/GB2-2025 nicolin-20250723T070234Z-1-001/GB2-2025 nicolin/GB2-2025/GB2/PR4/carbono_orgánico_suelo.tif"

Ten en cuenta que la siguiente tarea puede tardar varios minutos. ¡Ejecútala solo una vez!

Las unidades de la capa SOC no son las más comunes, para facilitar su uso vamos a usar un factor de conversión del cual obtenemos unidades familiares.

(stder_soc <- terra::rast(file)/10)
## class       : SpatRaster 
## size        : 972, 963, 1  (nrow, ncol, nlyr)
## resolution  : 0.002260748, 0.002389609  (x, y)
## extent      : -79.0102, -76.8331, 0.3613, 2.684  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## varname     : carbono_orgánico_suelo 
## name        : carbono_orgánico_suelo 
## min value   :                    0.0 
## max value   :                  262.3

##4. Exploración de datos

Un histograma puede ayudar con esto

## añade un titulo significativo y las unidades de SOC

(media_soc = mean(values(stder_soc), na.rm = T))
## [1] 71.41871
par(bg = "white")
options(scipen = 10)
hist(values(stder_soc), main= "Contenido de carbono orgánico del suelo en la fracción de tierra fina", 
     xlab= "g/kg", xlim = c(0, 120), col="lightblue", breaks=20)

Un resumen también es útil.

summary(stder_soc)
## Warning: [summary] used a sample
##  carbono_orgánico_suelo
##  Min.   :  0.0         
##  1st Qu.: 35.6         
##  Median : 56.4         
##  Mean   : 71.4         
##  3rd Qu.: 98.3         
##  Max.   :251.8

Ahora cambiamos el nombre del atributo SOC

(names(stder_soc) <-  "soc")
## [1] "soc"

Es hora de graficar. ### Ponemos los valores de SOC en una variable

val = values(stder_soc, na.rm = T)

creamos una paleta de colores

orangecyan = colorNumeric(c("orange","yellow2", "darkseagreen", "cyan" ), val,
                          na.color = "transparent")

Por ultimo la grafica.

#Grafica un mapa interactivo
stder <- st_transform(stder, crs = 4326)

leaflet::leaflet(stder) %>% 
  addTiles() %>% 
  setView(-74, 6, 9) %>% 
  addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
              opacity = 0.5, fillOpacity = 0.10,
              ) %>%
  addRasterImage(stder_soc, colors ="Spectral", opacity = 0.8)  %>%
  addLegend(pal = orangecyan, values = val, title = "Soil Organic Carbon (SOC) [g/kg]")

5. Muestreando el mundo

Supongamos que también estamos interesados en obtener una muestra del conjunto de datos de SOC. En este cuaderno, usaremos la biblioteca terra para este propósito.

set.seed(123446)
#convertimos al SRC de SOC
geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(stder_soc, geog))
## class       : SpatRaster 
## size        : 998, 936, 1  (nrow, ncol, nlyr)
## resolution  : 0.00232667, 0.00232667  (x, y)
## extent      : -79.0102, -76.83244, 0.3619834, 2.684  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        :      soc 
## min value   :   0.0000 
## max value   : 257.1497
(samples = terra::spatSample(geog.soc, 2000, "random", as.points=TRUE))
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 2000, 1  (geometries, attributes)
##  extent      : -79.00904, -76.8336, 0.3631468, 2.682837  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       :   soc
##  type        : <num>
##  values      :     0
##                111.6
##                    0

Nota dos cosas: (i) la clase del objeto en la salida y (ii) el CRS del conjunto de datos de la muestra.

Por lo tanto, necesitamos realizar varios procesos. Explica cada uno de los bloques de código que siguen.

(muestras <- sf::st_as_sf(samples))
## Simple feature collection with 2000 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -79.00904 ymin: 0.3631468 xmax: -76.8336 ymax: 2.682837
## Geodetic CRS:  GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]
## First 10 features:
##          soc                    geometry
## 1    0.00000  POINT (-78.98577 1.298468)
## 2  111.57671  POINT (-78.39712 1.666082)
## 3    0.00000   POINT (-78.7438 1.998796)
## 4   51.50857  POINT (-76.87083 1.801029)
## 5   79.23051 POINT (-77.48507 0.5050736)
## 6   25.33458 POINT (-78.66236 0.6423472)
## 7  113.99451  POINT (-78.48786 2.036022)
## 8    0.00000  POINT (-78.69494 1.880136)
## 9   51.72680 POINT (-77.82011 0.9262009)
## 10 135.99634  POINT (-77.85501 2.419923)
muestras %>% as_tibble()
## # A tibble: 2,000 × 2
##      soc              geometry
##    <dbl>           <POINT [°]>
##  1   0    (-78.98577 1.298468)
##  2 112.   (-78.39712 1.666082)
##  3   0     (-78.7438 1.998796)
##  4  51.5  (-76.87083 1.801029)
##  5  79.2 (-77.48507 0.5050736)
##  6  25.3 (-78.66236 0.6423472)
##  7 114.   (-78.48786 2.036022)
##  8   0    (-78.69494 1.880136)
##  9  51.7 (-77.82011 0.9262009)
## 10 136.   (-77.85501 2.419923)
## # ℹ 1,990 more rows

Omitimos los valores de NA

(nmuestras <- na.omit(muestras))
## Simple feature collection with 2000 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -79.00904 ymin: 0.3631468 xmax: -76.8336 ymax: 2.682837
## Geodetic CRS:  GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]
## First 10 features:
##          soc                    geometry
## 1    0.00000  POINT (-78.98577 1.298468)
## 2  111.57671  POINT (-78.39712 1.666082)
## 3    0.00000   POINT (-78.7438 1.998796)
## 4   51.50857  POINT (-76.87083 1.801029)
## 5   79.23051 POINT (-77.48507 0.5050736)
## 6   25.33458 POINT (-78.66236 0.6423472)
## 7  113.99451  POINT (-78.48786 2.036022)
## 8    0.00000  POINT (-78.69494 1.880136)
## 9   51.72680 POINT (-77.82011 0.9262009)
## 10 135.99634  POINT (-77.85501 2.419923)

seguimos con la visualización de las muestras limpias

longit <- st_coordinates(nmuestras)[,1]
latit <- st_coordinates(nmuestras)[,2]
soc <- nmuestras$soc

¿Qué obtuvimos?

summary(soc)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   35.23   54.41   71.40  100.67  228.83
length(soc)
## [1] 2000

¿Cuántas muestras no contienen datos validos?

id <- seq(1,2000,1)
(sitios <- data.frame(id, longit, latit, soc))

Ahora, graficamos las muestras

m <- leaflet() %>%
  addTiles() %>%  
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions())
m

6. Guardamos las muestras

#Usar la ruta deseada
sf::st_write(nmuestras, "C:/Users/vocav/Documents/GB2-2025 nicolin-20250723T070234Z-1-001/GB2-2025 nicolin/GB2-2025/GB2/PR4/soc_stder.gpkg", driver = "GPKG", append = F)

7. Leer los datos de entrada

list.files(path="C:/Users/vocav/Documents/GB2-2025 nicolin-20250723T070234Z-1-001/GB2-2025 nicolin/GB2-2025/GB2/PR4", pattern = "soc_stder.gpkg")
## [1] "soc_stder.gpkg"

Entonces, leemos las muestras usando sf:

samples <- sf::st_read("C:/Users/vocav/Documents/GB2-2025 nicolin-20250723T070234Z-1-001/GB2-2025 nicolin/GB2-2025/GB2/PR4/soc_stder.gpkg")
## Reading layer `soc_stder' from data source 
##   `C:\Users\vocav\Documents\GB2-2025 nicolin-20250723T070234Z-1-001\GB2-2025 nicolin\GB2-2025\GB2\PR4\soc_stder.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 2000 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -79.00904 ymin: 0.3631468 xmax: -76.8336 ymax: 2.682837
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

Leemos también el geopackage correspondiente a nuestra area de estudio.

munic <- sf::st_read("C:/Users/vocav/Documents/GB2-2025 nicolin-20250723T070234Z-1-001/GB2-2025 nicolin/GB2-2025/GB2/PR4/nariñito.gpkg")
## Reading layer `municipios' from data source 
##   `C:\Users\vocav\Documents\GB2-2025 nicolin-20250723T070234Z-1-001\GB2-2025 nicolin\GB2-2025\GB2\PR4\nariñito.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 64 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -79.01021 ymin: 0.3613481 xmax: -76.83368 ymax: 2.683898
## Geodetic CRS:  MAGNA-SIRGAS

8. Explorar los datos de entrada

Obtenemos un resumen de las muestras

summary(samples)
##       soc                    geom     
##  Min.   :  0.00   POINT        :2000  
##  1st Qu.: 35.23   epsg:NA      :   0  
##  Median : 54.41   +proj=long...:   0  
##  Mean   : 71.40                       
##  3rd Qu.:100.67                       
##  Max.   :228.83

Un histograma puede ayudar de nuevo

hist(samples$soc, main="Histograma SOC", xlab = "g/Kg", ylab = "frecuencia", col= "green") 
abline(h = 0, v = 28.03, col = "firebrick",lty = 2 )

Redondeamos los valores de SOC a dos digitos

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

Ahora podemos visualizar la distribucion espacial de las muestras

pal <- colorOptions(
  c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),
  domain = samples$soc,
)

reproyectamos las capas al mismo SRC

munic <- st_transform(munic, crs = 4326)
samples <- st_transform(samples, crs = 4326)
pal <- colorNumeric("viridis", domain = samples$soc)

leaflet()%>%
  addPolygons(
    data = munic,
    color = "cyan",
    # set the opacity of the outline
    opacity = 0.5,
    # 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 = samples$soc,
    position = "bottomleft",
    title = "SOC:",
    opacity = 0.9) %>%
  addProviderTiles("OpenStreetMap")  
st_geometry_type(munic)
##  [1] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
##  [6] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [11] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [16] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [21] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [26] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [31] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [36] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [41] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [46] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [51] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [56] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [61] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

9. Interpolación

9.1 Creación de objeto gstat

Para interpolar, primero necesitamos crear un objeto de clase gstat utilizando una función con el mismo nombre: gstat.

Un objeto gstat contiene toda la información necesaria para realizar la interpolación espacial, incluyendo:

La definición del modelo Los datos de calibración Según sus argumentos, la función gstat “interpreta” qué tipo de modelo de interpolación queremos usar:

Sin modelo de variograma → Inverse Distance Weighting (IDW) Modelo de variograma, sin covariables → Ordinary Kriging (OK) Vamos a utilizar tres parámetros en la función gstat:

formula: la “fórmula” de predicción que especifica la variable dependiente y las independientes (también llamadas covariables). data: los datos de calibración (también conocidos como datos de entrenamiento). model: el modelo de variograma.”

9.2 Interpolación IDW

Para interpolar el carbono orgánico del suelo (SOC) utilizando el método IDW, creamos el siguiente objeto gstat, especificando solo la fórmula y los datos.

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

Ahora que nuestro modelo de interpolación g1 está definido, podemos utilizar la función predict para realizar la interpolación, es decir, estimar los valores de SOC.

La función predict acepta:

Un raster (stars object), como dem. Un modelo (gstat object), como g1. El raster cumple dos propósitos:

Especificar las ubicaciones donde queremos hacer predicciones (en todos los métodos). Especificar valores de covariables (solo en Kriging Universal). Vamos a utilizar la librería terra para crear un objeto raster con valores de celda iguales a 1. La extensión de este raster debe cubrir nuestra área de interés.

Nota: Puedes usar la metadata de cualquier raster que cubra tu departamento (por ejemplo, un modelo digital de elevación) para definir los parámetros de creación del raster.”

(rrr <- terra::rast(xmin=-79.0102, xmax=-76.8331, ymin= 0.3613 , ymax=2.684, nrows=370, ncols = 329, vals = 1, crs = "epsg:4326"))
## class       : SpatRaster 
## size        : 370, 329, 1  (nrow, ncol, nlyr)
## resolution  : 0.006617325, 0.006277568  (x, y)
## extent      : -79.0102, -76.8331, 0.3613, 2.684  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : lyr.1 
## min value   :     1 
## max value   :     1

Ahora, necesitamos convertir este SpatRaster en un objeto stars:

stars.rrr <- stars::st_as_stars(rrr)

La siguiente expresión interpola los valores de SOC según el modelo definido en g1 y el ráster definido en stars.rrr

z1 = predict(g1, stars.rrr)
## [inverse distance weighted interpolation]

¿Qué es z1?

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                   Min.  1st Qu. Median     Mean  3rd Qu.     Max.   NA's
## var1.pred  0.005282132 47.17663 62.349 72.20531 89.39917 221.9297      0
## var1.var            NA       NA     NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -79.01  0.006617 WGS 84 [x]
## y    1 370  2.684 -0.006278 WGS 84 [y]

Toma nota de los nombres de los dos atributos incluidos dentro del objeto z1.

Podemos seleccionar solo el primer atributo y renombrarlo:

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  0.005282132 47.17663 62.349 72.20531 89.39917 221.9297      0
## var1.var            NA       NA     NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -79.01  0.006617 WGS 84 [x]
## y    1 370  2.684 -0.006278 WGS 84 [y]
names(z1) = "soc"

Para imprimir la superficie interpolada, necesitaremos una paleta de colores

paleta <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), domain = range(samples$soc, na.rm = TRUE))
paleta <- colorOptions(
  palette = c("orange", "yellow", "cyan", "green"),
  domain = range(samples$soc, na.rm = TRUE)
)

Ahora visualicemos el resultado de la interpolación:

z1 <- z1["soc"]

msoc <- leaflet() %>%
  addTiles() %>%
  addStarsImage(
    z1,
    opacity = 0.7,
    colorOptions = list(palette = c("orange", "yellow", "cyan", "green"), domain = range(samples$soc, na.rm = TRUE))
  ) %>%
  addCircleMarkers(
    data = samples,
    radius = 1.5,
    label = ~soc,
    color = ~pal(soc),
    fillOpacity = 1,
    stroke = FALSE
  ) %>%
  addLegend(
    "bottomright",
    pal = pal,
    values = samples$soc,
    title = "IDW SOC interpolation [%]"
  )
msoc #imprimir el mapa

9.3 Interpolación OK

os métodos de Kriging requieren un modelo de variograma. El modelo de variograma es una forma objetiva de cuantificar el patrón de autocorrelación en los datos y asignar pesos en consecuencia al hacer predicciones.

Como primer paso, podemos calcular y examinar el variograma empírico utilizando la función variogram.

Esta función requiere dos argumentos:

formula: especifica la variable dependiente y las covariables, al igual que en gstat. data: la capa de puntos con la variable dependiente y las covariables como atributos de los puntos. La siguiente expresión calcula el variograma empírico de samples, sin covariables:

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

Imprimimos el variograma empirico

plot(v_emp_ok)

Ajustamos el variograma a una función

v_mod_ok = automap::autofitVariogram(soc ~ 1, as(samples, "Spatial"), model = c("Sph", "Exp", "Gau"))

Imprimimos el modelo

plot(plot(v_mod_ok)
)

El objeto resultante es en realidad una lista con varios componentes, incluyendo el variograma empírico y el modelo de variograma ajustado.

El componente $var_model del objeto resultante contiene el modelo real:

v_mod_ok$var_model
##   model     psill    range
## 1   Nug  287.9209   0.0000
## 2   Sph 3606.9052 145.4271

Ahora, el modelo de variograma puede ser pasado a la función gstat, y podemos continuar con la interpolación mediante Kriging Ordinario.

g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2= predict(g2, stars.rrr)
## [using ordinary kriging]

Nuevamente, seleccionaremos el atributo de valores predichos y lo renombraremos.

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

Las predicciones de OK podran verse en la siguiente impresión

paleta <- colorNumeric(
     palette = c("orange", "yellow", "cyan", "green"),
     domain = range(z2[[1]], na.rm = TRUE)  # o z2$soc si existe
 )
mp <- leaflet() %>%
   addTiles() %>%
   addStarsImage(
     z2, 
     opacity = 0.7,
     colorOptions = colorOptions(
       palette = c("orange", "yellow", "cyan", "green"),
       domain = range(z2$soc, na.rm = TRUE)
     )
   ) %>%
   addCircleMarkers(
     data = samples,
     radius = 1.5, 
     label = ~soc,
     color = ~paleta(soc),
     fillOpacity = 1,
     stroke = FALSE
   ) %>%
   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
mp

10. Evaluación de resultados

10.1 Evaluación cualitativa

Comparamos las graficas

# Establecemos la paleta de colores
paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), 
                        domain = range(z1[[1]], z2[[1]], na.rm = TRUE), 
                        na.color = "transparent")

Crear mapa interactivo

convertimos nuestros datos en objetos raster

z1_terra <- rast(z1) 

z2_terra <- rast(z2["soc"])

guardamos las nuevas capas raster

writeRaster(z1_terra, "IDW_interpolated.tif", overwrite=TRUE)
writeRaster(z2_terra, "OK_interpolated.tif", overwrite=TRUE)  

las añadimos al workspace

z1_terra <- rast("IDW_interpolated.tif")
z2_terra <- rast("OK_interpolated.tif")

# Convertir a RasterLayer (si es necesario para compatibilidad con leaflet)
z1_raster <- raster::raster(z1_terra)
z2_raster <- raster::raster(z2_terra)

# Crear el mapa interactivo con leaflet
mapa_interactivo <- leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%  # Agrega un mapa base ligero
  addRasterImage(z1_raster, opacity = 0.7, group = "IDW") %>%
  addRasterImage(z2_raster, opacity = 0.7, group = "OK") %>%
  addLayersControl(
    overlayGroups = c("IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% #añadimos la leyenda
  addLegend("bottomright", pal = paleta, values = range(z1[[1]], na.rm = TRUE),
            title = "Soil organic carbon [%]")
mapa_interactivo

10.2 Validación cruzada

Hemos estimado las superficies de SOC utilizando dos métodos diferentes: IDW y Kriging Ordinario. Aunque es útil examinar y comparar los resultados gráficamente, también necesitamos una forma objetiva de evaluar la precisión de la interpolación. De esta manera, podemos elegir objetivamente el método más preciso entre los métodos de interpolación disponibles.

En la Validación Cruzada Leave-One-Out:

Se elimina un punto de los datos de calibración. Se realiza una predicción para ese punto. Se repite el proceso para todos los puntos. Al final, obtenemos una tabla con un valor observado y un valor predicho para todos los puntos.

Podemos ejecutar la Validación Cruzada Leave-One-Out utilizando la función gstat.cv, que acepta un objeto gstat.

cv2 = gstat.cv(g2)
cv2 = st_as_sf(cv2)
sp::bubble(as(cv2[, "residual"], "Spatial"))

Este es el valor de RMSE para la interpolación OK

sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))
## [1] 22.36934

Crear el modelo IDW

g_idw <- gstat(formula = soc ~ 1, data = samples, nmax = 15, set = list(idp = 2))
cv_idw <- gstat.cv(g_idw)
cv_idw <- st_as_sf(cv_idw)  # Convertir a sf
sp::bubble(as(cv_idw[, "residual"], "Spatial"))

Este es el valor RMSE para la interpolación de IDW

rmse_idw <- sqrt(sum((cv_idw$var1.pred - cv_idw$observed)^2) / nrow(cv_idw))
print(rmse_idw)
## [1] 22.20898

Al comparar los valores de RMSE podemos concluir que el mejor metodo de interpolación para este proyecto es Ordinary Kriging

Referencias

Lizarazo, I. 2025. Spatial interpolation. Available at: https://rpubs.com/ials2un/Spatial_Interpolation