1. Introducción

En este cuaderno ilustraremos el contenido de carbono en suelo del departamento de Huila, para ello utilizaremos dos técnicas de interpolación espacial: Inverse Distance Weighted (IDW) que es una técnica determinística y Ordinary Kriging (OK). que es una técnica probabilística.Para ambas técnicas se utilizaron datos para obtener una superficie continua del contenido de carbono orgánico del suelo (SOC,) a una profundidad de 15-30 cm, a partir de muestras obtenidas de SoilGrids 250 m.

2. Requerimientos

Para iniciar con las ilustraciones empezaremos por limpiar (variables y datos cargados) el espacio de trabajo:

# Uncomment at first-time run
rm(list=ls())

Empezaremos con la instalación de paquetes requeridos para trabajar las funciones necesarias para nuestras ilustraciones:

# DO NOT INSTALL  PACKAGES THAT YOU ALREADY INSTALLED
#paquetes = c('terra', 'sf', 'sp', 'stars', 'gstat', 'automap', 'leaflet', 'leafem')
#install.packages(paquetes)
#

Cargamos nuestras librerías:

library(terra)
## terra 1.8.21
library(sf)
## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(stars)
## Cargando paquete requerido: abind
library(sp)
library(gstat)
library(automap)
library(RColorBrewer)
library(leaflet)
library(leafem)

3. Lectura de datos

Ya con las respectivas librerías cargadas empezamos a cagar nuestros datos, para nuestro caso necesitamos un archivo .gpkg(ya creado anteriormente) y para ello lo primero es saber el nombre específico dentro de la carpeta en la que este depositado:

list.files(path="datos", pattern = "*.gpkg")
## [1] "Huila.gpkg"     "soc_HUila.gpkg"

Ya obtenido el contenido de archivos .gpkg disponibles tomamos el deseado, en nuestro caso “soc_HUila.gpkg” y lo llamaremos con el siguiente código:

samples <- sf::st_read("datos/soc_HUila.gpkg")
## Reading layer `soc_HUila' from data source 
##   `C:\Users\usuario\Documents\GB2\INFORME_2\datos\soc_HUila.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 2000 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -76.49652 ymin: 1.530938 xmax: -74.50194 ymax: 3.527839
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

Habiendo leído los datos para a librería sf y nombrándolos samples, haremos lo mismos para leerlo como varíale munic:

munic <- sf::st_read("datos/soc_HUila.gpkg")
## Reading layer `soc_HUila' from data source 
##   `C:\Users\usuario\Documents\GB2\INFORME_2\datos\soc_HUila.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 2000 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -76.49652 ymin: 1.530938 xmax: -74.50194 ymax: 3.527839
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

Para ambos casos tenemos que estamos leyendo 2000 datos con geometría de punto.

4. Reconocimiento de datos

Lo primero que aremos con los datos es identificar los valores entre los que están los datos:

summary(samples)
##       soc                    geom     
##  Min.   :  0.00   POINT        :2000  
##  1st Qu.: 23.65   epsg:NA      :   0  
##  Median : 39.98   +proj=long...:   0  
##  Mean   : 40.55                       
##  3rd Qu.: 55.37                       
##  Max.   :100.62

Ya sabiendo que los datos rondan entre los 0.00 y los 100.62 datos de contenido de carbono orantico en suelo generaremos un histograma para ver que tan frecuentes son los datos que estamos viendo:

hist(samples$soc)

Con lo anterior vemos que la frecuencia de datos mas recurrente esta entre los 20-30 y resto de los datos entre 10-60 se presentan casi con la misma frecuencia.

Posterior a esto ilustraremos los puntos del SOC para nuestro departamento del Huila:

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

Agregaremos una paleta de colores:

pal <- colorNumeric(
  c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),
  # colors depend on the count variable
  domain = samples$soc,
  )

Traeremos la ilustración con los parámetros que indicamos a continuación:

leaflet() %>%
  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")

Como ya nos lo indicaba el histograma, la coloración de puntos mas presente es la de los valores entre 20-30 y los siguientes colores mas presentes están los de 10-60:

5. Interpolación de datos

Como la ilustración anterior nos arrojo solo la ubicación y valores de los puntos de muestra contenidos y el resto del espacio solo están en blanco pues son espacios no muestreados, utilizaremos el método de interpolación de datos para lograr estimar los valores desconocidos de los puntos no muestreados a partir de los que ya conocemos y tener una superficie continua.

5.1 Interpolación IDW

gstat: Es una función del paquete gstat que se utiliza para crear modelos de interpolación espacial.

formula = soc ~ 1: Especifica la variable a interpolar (soc, que es el contenido de carbono orgánico del suelo) y el modelo a usar. En este caso, ~ 1 indica que no hay variables explicativas adicionales (solo se usa la ubicación de las muestras).

data = samples: Es el conjunto de datos que contiene las muestras de SOC y sus coordenadas.

g1 = gstat(formula = soc ~ 1, data = samples)
 terra::rast: Crea un objeto raster (rejilla) utilizando el paquete terra.

xmin, xmax, ymin, ymax: Define los límites geográficos de la rejilla (en este caso, el departamento del Huila).

nrows y ncols: Especifican el número de filas y columnas de la rejilla. Esto determina la resolución espacial del raster.

vals = 1: Asigna un valor inicial a todas las celdas del raster (en este caso, 1).

crs = “epsg:4326”: Define el sistema de referencia de coordenadas (CRS) como WGS84, que es común para datos geográficos.

rrr es un objeto raster que cubre el área de interés (Huila) y servirá como base para la interpolación.

(rrr <- terra::rast(xmin=-76.49652, xmax=-74.50194, ymin=1.530938, ymax=3.527839, nrows=370, ncols = 329, vals = 1, crs = "epsg:4326"))
## class       : SpatRaster 
## dimensions  : 370, 329, 1  (nrow, ncol, nlyr)
## resolution  : 0.006062553, 0.00539703  (x, y)
## extent      : -76.49652, -74.50194, 1.530938, 3.527839  (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::st_as_stars: Convierte el objeto raster (rrr) en un objeto stars, que es un formato compatible con el paquete gstat para realizar la interpolación.

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

predict: Utiliza el modelo g1 (IDW) para predecir valores de SOC en cada celda de la rejilla stars.rrr.

g1: Es el modelo de interpolación IDW creado en el Paso 1.

stars.rrr: Es la rejilla de interpolación creada en el Paso 3.

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

z1 es un objeto stars que contiene los valores interpolados de SOC en cada celda de la rejilla.

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                 Min.  1st Qu.   Median     Mean  3rd Qu.     Max.   NA's
## var1.pred  0.4084332 30.79741 40.68857 40.72309 51.20584 96.58159      0
## var1.var          NA       NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329  -76.5  0.006063 WGS 84 [x]
## y    1 370  3.528 -0.005397 WGS 84 [y]

is.na(z1[[1]]): Identifica las celdas de la rejilla que tienen valores NA (faltantes).

which: Encuentra los índices de las celdas con valores NA.

z1[[1]][a] = 0.0: Asigna un valor de 0.0 a las celdas que tienen NA.

a <- which(is.na(z1[[1]]))
z1[[1]][a] = 0.0

Se eliminan los valores NA de la rejilla interpolada, asignándoles un valor de 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.4084332 30.79741 40.68857 40.72309 51.20584 96.58159      0
## var1.var          NA       NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329  -76.5  0.006063 WGS 84 [x]
## y    1 370  3.528 -0.005397 WGS 84 [y]

names(z1): Asigna un nombre a la capa de datos en el objeto z1.

“soc”: El nombre asignado es soc, que representa el contenido de carbono orgánico del suelo.

names(z1) = "soc"

La capa interpolada en z1 ahora se llama soc. Agregaremos una paleta de colores:

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  # Print the map

5.2 Interpolación OK

variogram: Es una función del paquete gstat que calcula el variograma experimental, que describe cómo varía la similitud entre los valores de SOC en función de la distancia.

soc ~ 1: Especifica que la variable a analizar es soc (contenido de carbono orgánico del suelo) y que no hay variables explicativas adicionales.

data = samples: Es el conjunto de datos que contiene las muestras de SOC y sus coordenadas.

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

v_emp_ok es un objeto que contiene el variograma experimental, que muestra la variabilidad espacial de los datos.

Visualizaremos los datos.

plot(v_emp_ok)

automap::autofitVariogram: Es una función del paquete automap que ajusta automáticamente un modelo de variograma a los datos.

soc ~ 1: Especifica la variable a analizar (soc).

as(samples, “Spatial”): Convierte el objeto samples (que es un sf) a un objeto Spatial (requerido por autofitVariogram).

#make sure you understand the parameters needed to run this fitting function
v_mod_ok = automap::autofitVariogram(soc ~ 1, as(samples, "Spatial"))

v_mod_ok es un objeto que contiene el modelo de variograma ajustado, que describe la estructura espacial de los datos.

plot(v_mod_ok)

v_mod_ok$var_model: Extrae el modelo de variograma ajustado del objeto v_mod_ok.

v_mod_ok$var_model
##   model     psill    range
## 1   Nug  22.10619  0.00000
## 2   Sph 300.64894 71.83277

gstat: Es una función del paquete gstat que crea un modelo de interpolación.

formula = soc ~ 1: Especifica la variable a interpolar (soc) y que no hay variables explicativas adicionales.

model = v_mod_ok$var_model: Usa el modelo de variograma ajustado para la interpolación.

data = samples: Es el conjunto de datos que contiene las muestras de SOC y sus coordenadas.

g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2= predict(g1, stars.rrr)
## [inverse distance weighted interpolation]

g2 es un objeto que contiene el modelo de Kriging Ordinario (OK), que se utilizará para predecir valores de SOC en puntos no muestreados.

is.na(z2[[1]]): Identifica las celdas de la rejilla que tienen valores NA (faltantes).

which: Encuentra los índices de las celdas con valores NA.

z2[[1]][a] = 0.0: Asigna un valor de 0.0 a las celdas que tienen NA.

names(z2) = “soc”: Asigna un nombre a la capa de datos en el objeto z2.

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

6. Evaluación

6.1 Evaluación grafica

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

6.2 Evaluación numérica

# Validación cruzada para IDW
#cv_idw <- gstat.cv(g1)
#print("Resultados de validación cruzada para IDW:")
#print(summary(cv_idw))

# Validación cruzada para Kriging
#cv_ok <- gstat.cv(g2)
#print("Resultados de validación cruzada para Kriging:")
#print(summary(cv_ok))

7

write_stars(
  z1, dsn = "IDW_SOC_Huila.tif", layer = 1
)
write_stars(
  z2, dsn = "OK_SOC_Huila.tif", layer = 1
)