“Muestras de carbono organico de los suelos de El Carmen, Norte de Santander”

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/galle/OneDrive - Universidad Nacional de Colombia/Universidad Nacional/UN 2024-2/Geomatica/R/El_carmen_r/el_carmen.shp"))

## Reading layer `el_carmen' from data source 
##   `C:\Users\galle\OneDrive - Universidad Nacional de Colombia\Universidad Nacional\UN 2024-2\Geomatica\R\El_carmen_r\el_carmen.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 91 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 4936439 ymin: 2486485 xmax: 4999704 ymax: 2584646
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional

## Simple feature collection with 1 feature and 91 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 4936439 ymin: 2486485 xmax: 4999704 ymax: 2584646
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
##   DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR MPIO_CDPMP VERSION       AREA LATITUD
## 1         54        245  EL CARMEN      54245    2018 1727742958 8.86799
##    LONGITUD STCTNENCUE STP3_1_SI STP3_2_NO STP3A_RI STP3B_TCN STP4_1_SI
## 1 -73.34586       5461       424      5037      424         0       424
##   STP4_2_NO STP9_1_USO STP9_2_USO STP9_3_USO STP9_4_USO STP9_2_1_M STP9_2_2_M
## 1      5037       4390        137        931          3          1        126
##   STP9_2_3_M STP9_2_4_M STP9_2_9_M STP9_3_1_N STP9_3_2_N STP9_3_3_N STP9_3_4_N
## 1         10          0          0          2         90         78         20
##   STP9_3_5_N STP9_3_6_N STP9_3_7_N STP9_3_8_N STP9_3_9_N STP9_3_10 STP9_3_99
## 1        139        526         26          0          4        46         0
##   STVIVIENDA STP14_1_TI STP14_2_TI STP14_3_TI STP14_4_TI STP14_5_TI STP14_6_TI
## 1       4527       4152         78        179        104          9          5
##   STP15_1_OC STP15_2_OC STP15_3_OC STP15_4_OC TSP16_HOG STP19_EC_1 STP19_ES_2
## 1       3591          5        385        546      3630       2852        739
##   STP19_EE_1 STP19_EE_2 STP19_EE_3 STP19_EE_4 STP19_EE_5 STP19_EE_6 STP19_EE_9
## 1       2068        764          6          4          0          2          8
##   STP19_ACU1 STP19_ACU2 STP19_ALC1 STP19_ALC2 STP19_GAS1 STP19_GAS2 STP19_GAS9
## 1       1530       2061       1387       2204         21       3545         25
##   STP19_REC1 STP19_REC2 STP19_INT1 STP19_INT2 STP19_INT9 STP27_PERS STPERSON_L
## 1       1505       2086        142       3418         31      12001         42
##   STPERSON_S STP32_1_SE STP32_2_SE STP34_1_ED STP34_2_ED STP34_3_ED STP34_4_ED
## 1      11959       6261       5740       2306       2493       1834       1566
##   STP34_5_ED STP34_6_ED STP34_7_ED STP34_8_ED STP34_9_ED STP51_PRIM STP51_SECU
## 1       1367        999        729        448        259       5552       2615
##   STP51_SUPE STP51_POST STP51_13_E STP51_99_E Shape_Leng Shape_Area Codigo_Mun
## 1        450         82       2030        190   3.214756  0.1419942      54245
##                         geometry
## 1 POLYGON ((4999571 2583470, ...

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 
## -8222911.2   935342.9 -8166003.6  1034252.3

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 
## -8222911.2  1034252.3 -8166003.6   935342.9

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/galle/OneDrive - Universidad Nacional de Colombia/Universidad Nacional/UN 2024-2/Geomatica/rsoc_igh_15_30.tif"

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

gdal_translate(paste0(sg_url,datos), file ,
               tr=c(250,250),
               projwin=bbx, 
               projwin_srs =repro)

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 
## dimensions  : 396, 228, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8223000, -8166000, 935500, 1034500  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source(s)   : memory
## varname     : rsoc_igh_15_30 
## name        : rsoc_igh_15_30 
## min value   :           10.5 
## max value   :           76.2

##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] 27.63008

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)

##  rsoc_igh_15_30 
##  Min.   :10.50  
##  1st Qu.:23.10  
##  Median :26.90  
##  Mean   :27.63  
##  3rd Qu.:31.10  
##  Max.   :76.20  
##  NA's   :27

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

## Warning: sf layer is not long-lat data

## Warning: sf layer has inconsistent datum (+proj=tmerc +lat_0=4 +lon_0=-73 +k=0.9992 +x_0=5000000 +y_0=2000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs).
## Need '+proj=longlat +datum=WGS84'

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 
## dimensions  : 407, 267, 1  (nrow, ncol, nlyr)
## resolution  : 0.002184878, 0.002184878  (x, y)
## extent      : -73.58484, -73.00148, 8.403826, 9.293072  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        :      soc 
## min value   : 10.61527 
## max value   : 73.71278

(samples = terra::spatSample(geog.soc, 2000, "random", as.points=TRUE))

##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 2000, 1  (geometries, attributes)
##  extent      : -73.58375, -73.00257, 8.404919, 9.291979  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       :   soc
##  type        : <num>
##  values      : 43.15
##                   NA
##                46.16

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: -73.58375 ymin: 8.404919 xmax: -73.00257 ymax: 9.291979
## 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  43.14563 POINT (-73.09871 9.110634)
## 2        NA POINT (-73.02005 8.968617)
## 3  46.15673 POINT (-73.17736 8.442061)
## 4  26.51857 POINT (-73.30408 8.813491)
## 5        NA POINT (-73.58156 9.097525)
## 6        NA POINT (-73.57282 9.248282)
## 7  27.21408 POINT (-73.42644 8.699877)
## 8        NA POINT (-73.57719 8.763239)
## 9  36.05132 POINT (-73.05501 8.782902)
## 10 31.55754 POINT (-73.39585 8.636516)

muestras %>% as_tibble()

## # A tibble: 2,000 × 2
##      soc             geometry
##    <dbl>          <POINT [°]>
##  1  43.1 (-73.09871 9.110634)
##  2  NA   (-73.02005 8.968617)
##  3  46.2 (-73.17736 8.442061)
##  4  26.5 (-73.30408 8.813491)
##  5  NA   (-73.58156 9.097525)
##  6  NA   (-73.57282 9.248282)
##  7  27.2 (-73.42644 8.699877)
##  8  NA   (-73.57719 8.763239)
##  9  36.1 (-73.05501 8.782902)
## 10  31.6 (-73.39585 8.636516)
## # ℹ 1,990 more rows

Omitimos los valores de NA

(nmuestras <- na.omit(muestras))

## Simple feature collection with 1748 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -73.57719 ymin: 8.404919 xmax: -73.00694 ymax: 9.291979
## 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  43.14563 POINT (-73.09871 9.110634)
## 3  46.15673 POINT (-73.17736 8.442061)
## 4  26.51857 POINT (-73.30408 8.813491)
## 7  27.21408 POINT (-73.42644 8.699877)
## 9  36.05132 POINT (-73.05501 8.782902)
## 10 31.55754 POINT (-73.39585 8.636516)
## 11 26.04888 POINT (-73.15333 8.798197)
## 12 30.82602 POINT (-73.24509 8.485759)
## 13 20.44965  POINT (-73.5619 8.404919)
## 14 23.28090 POINT (-73.28879 8.494499)

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. 
##   12.53   23.76   27.27   28.03   31.52   63.85

length(soc)

## [1] 1748

¿Cuántas muestras no contienen datos validos?

id <- seq(1,1748,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/galle/OneDrive - Universidad Nacional de Colombia/Universidad Nacional/UN 2024-2/Geomatica/R/El_carmen_r/soc_stder.gpkg", driver = "GPKG", append = F)

7. Leer los datos de entrada

list.files(path="C:/Users/galle/OneDrive - Universidad Nacional de Colombia/Universidad Nacional/UN 2024-2/Geomatica/R/El_carmen_r", pattern = "soc_stder.gpkg")

## [1] "soc_stder.gpkg"

Entonces, leemos las muestras usando sf:

samples <- sf::st_read("C:/Users/galle/OneDrive - Universidad Nacional de Colombia/Universidad Nacional/UN 2024-2/Geomatica/R/El_carmen_r/soc_stder.gpkg")

## Reading layer `soc_stder' from data source 
##   `C:\Users\galle\OneDrive - Universidad Nacional de Colombia\Universidad Nacional\UN 2024-2\Geomatica\R\El_carmen_r\soc_stder.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1748 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -73.57719 ymin: 8.404919 xmax: -73.00694 ymax: 9.291979
## 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/galle/OneDrive - Universidad Nacional de Colombia/Universidad Nacional/UN 2024-2/Geomatica/R/El_carmen_r/el_carmen.gpkg")

## Reading layer `el_carmen' from data source 
##   `C:\Users\galle\OneDrive - Universidad Nacional de Colombia\Universidad Nacional\UN 2024-2\Geomatica\R\El_carmen_r\el_carmen.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1 feature and 91 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 4936439 ymin: 2486485 xmax: 4999704 ymax: 2584646
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional

8. Explorar los datos de entrada

Obtenemos un resumen de las muestras

summary(samples)

##       soc                   geom     
##  Min.   :12.53   POINT        :1748  
##  1st Qu.:23.76   epsg:NA      :   0  
##  Median :27.27   +proj=long...:   0  
##  Mean   :28.03                       
##  3rd Qu.:31.52                       
##  Max.   :63.85

Un histograma puede ayudar de nuevo

hist(samples$soc)

Redondeamos los valores de SOC a dos digitos

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

Ahora podemos visualizar la distribucion espacial de las muestras

pal <- colorNumeric(
  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)

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

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=-73.577194, xmax=-73.006941, ymin=8.404919, ymax=9.291979, nrows=370, ncols = 329, vals = 1, crs = "epsg:4326"))

## class       : SpatRaster 
## dimensions  : 370, 329, 1  (nrow, ncol, nlyr)
## resolution  : 0.001733292, 0.002397459  (x, y)
## extent      : -73.57719, -73.00694, 8.404919, 9.291979  (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  12.75162 25.79533 28.04482 27.73713 29.53867 62.43863      0
## var1.var         NA       NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -73.58  0.001733 WGS 84 [x]
## y    1 370  9.292 -0.002397 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  12.75162 25.79533 28.04482 27.73713 29.53867 62.43863      0
## var1.var         NA       NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -73.58  0.001733 WGS 84 [x]
## y    1 370  9.292 -0.002397 WGS 84 [y]

names(z1) = "soc"

Para imprimir la superficie interpolada, necesitaremos una paleta de colores

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

Ahora visualicemos el resultado de la interpolación:

z1 <- z1["soc"]

m <- leaflet() %>%
  addTiles() %>%
  leafem:::addGeoRaster(
    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 = ~paleta(soc),
    fillOpacity = 1,
    stroke = FALSE
  ) %>%
  addLegend(
    "bottomright",
    pal = paleta,
    values = samples$soc,
    title = "IDW SOC interpolation [%]"
  )
m #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  0.8028901 0.000000
## 2   Exp 39.6389427 2.559963

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

m <- leaflet() %>%
  addTiles() %>%
  leafem:::addGeoRaster(
    z2, 
    opacity = 0.7,
    colorNumeric = colorNumeric(
      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 write_stars.stars(x, dsn = fl): all but first attribute are ignored

m

10. Evaluación de resultados

10.1 Evaluación cualitativa

Comparamos las graficas

colores <- colorNumeric(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

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] 3.903031

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] 4.023636

####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