1. Introducción:

En este cuaderno aprenderemos a interpolar datos obtenidos de SoilGrids, tales como variables y propiedades con interés agronómico. Para el ejemplo del presente cuaderno, tomaremos la propiedad del pH del agua en el departamento del Huila; emplearemos dos métodos de interpolación, el IDW (Inverso De la Distancia), el cual es un método determinista, y el Kriging, el cual es un método estocástico, este último se considera más preciso que el anterior debido a que el Kriging, utiliza modelos estadísticos para analizar los datos y estimar los valores en ubicaciones donde no hay muestreo. En relación con los anteriores cuadernos, utilizaremos nuevas librerias, las cuales son “stars” y “automap”.

2. Parte Práctica

En primer lugar, y como es costumbre antes de realizar un cuaderno, limpiaremos la memoria en la consola, esto es importante para que no hallan errores a causa de conflictos con objetos o variables asignadas con el mismo nombre.

## Escribelo o copialo directamente en la consola
## rm(list=ls())
  • (Como dato curioso, otra alternativa es hacer click, en la escoba que se encuentra en la parte superior derecha donde se almacenan los objetos y variables asignadas).

Posteriormente de limpiar la consola, cargamos las librerias que vamos a utilizar, para el presente cuaderno, ya deben estar descargadas y instaladas en su Rstudio con anterioridad.

library(sp)
library(terra)
library(sf)
library(stars)
library(gstat)
library(automap)
library(leaflet)
library(leafem)
library(ggplot2)
library(dplyr)
library(curl)
h <- new_handle()
handle_setopt(h, http_version = 2)

3. Lectura de los Datos

Leemos el archivo raster que obuvimos por medio del cuaderno pasado o ya sea descargado directamente de SoilGrids.

archivo <- ("phh2o_igh_15_30.tif")
(phh2o <- rast(archivo))
## class       : SpatRaster 
## size        : 1020, 967, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8531500, -8289750, 173000, 428000  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source      : phh2o_igh_15_30.tif 
## name        : phh2o_igh_15_30

Convertimos los datos del pH del agua en un porcentaje, y anotamos el factor de escala que aparece en SoilGrids de la propiedad.

phh2o.perc <- phh2o/10

Para no tener inconvenientes a futuro, que no hallan errones ni confusiones del programa, transformamos los datos a un sistema WGS84.

geog = "+proj=longlat +datum=WGS84"
(geog.phh2o = project(phh2o.perc, geog))
## class       : SpatRaster 
## size        : 1030, 998, 1  (nrow, ncol, nlyr)
## resolution  : 0.002224706, 0.002224706  (x, y)
## extent      : -76.63117, -74.41092, 1.553343, 3.844789  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : phh2o_igh_15_30 
## min value   :        4.412293 
## max value   :        7.418550

Usando la libreria “stars”, (la cual es una libreria para trabajar con datos espacio-temporales tanto raster como vectoriales), transformamos la capa “SpatRaster” en un objeto “stars”

stars.phh2o =st_as_stars(geog.phh2o)

Hacemos un ploteo para observar nuestro objeto “stars” recién creado, este objeto tiene los datos que necesitamos para hacer la interpolación,por último, se tiene que revisar los valores minímos y máximos que comprende el objeto para colocar ese rango en el dominio del ploteo.

m <- leaflet() %>%
  addTiles() %>%
  leafem:::addGeoRaster(
    stars.phh2o,
    opacity = 0.8,
    colorOptions = colorOptions(palette = c("orange", "yellow","cyan", "green"),
                                domain = 4.3:7.5)
  )
m

4. Muestreo al azar

Realizaremos un muestreo de 800 sitios al azar reales utilizando una muestra aleatoria

set.seed(123456)

(samples <- spatSample(geog.phh2o, 800, "random", as.points=TRUE))
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 800, 1  (geometries, attributes)
##  extent      : -76.63006, -74.41425, 1.558904, 3.839228  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : phh2o_igh_15_30
##  type        :           <num>
##  values      :           5.066
##                            5.5
##                          5.284

Podemos observar que el muestreo arroja un objeto de clase “SpatVector”, con geometria de puntos, un sistema de referencia de coordenadas WGS84, con 800 sitios al azar seleccionados de 1 atributo y coordenadas en decimales.

Ahora convertiremos el objeto “SpatVector” en un “Simple Feature”.

(muestras <- sf::st_as_sf(samples))

Y omitimos las muestras que no tienen datos.

nmuestras <- na.omit(muestras)

Visualizamos las muestras.

longit <- st_coordinates(muestras)[,1]
latit <- st_coordinates(muestras)[,2]
phh2o <- muestras$phh2o_igh_15_30
id <- seq(1,800,1) 
(sitios <- data.frame(id, longit, latit, phh2o))

Y nuevamente omitimos las muestras que no tienen valores.

sitios <- na.omit(sitios)

Con la función “head”, podemos visualizar unicamente la 6 primeras muestras.

head(sitios)

Y ploteamos para observar las muestras con un gráfico más elaborado e interactivo en referencia con el anterior.

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      stars.phh2o,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 8:130)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$phh2o, clusterOptions = markerClusterOptions())
m

Luego de completar esta parte del cuaderno, ya estamos lo suficientemente preparados para empezar a hacer las interpolaciones por IDW y Kriging Ordinario.

5. Interpolación

5.1 Creación del objeto gstat

Para poder realizar la interpolación necesitamos emplear la función gstat, con esta función crearemos un objeto de clase gstat y podremos proceder.

Los objetos gstat contienen la información necesaria para realizar la interpolación espacial.

-la definición del modelo. -la calibración de los datos.

Con base en esto, la interpolación espacial predice que tipo de interpolación queremos utilizar.

-El IDW, el cual es un modelo sin variogramas. -El Kriging Ordinario, el cual es modelo con variogramas pero sin covariables

Vamos a utilizar 3 parametros de la función gstat.

  • fórmula: la “fórmula” de predicción que especifica las variables dependientes e independientes (que también se conoce como covariables). -datos: los datos de calibración (a.k.a datos del tren) -modelo: el modelo de variograma

5.2 Interpolación por el método IDW

Ahora crearemos un objeto gstat especificando la formula y los datos para interpolar el pH del agua de Huila por medio del método IDW.

g1 = gstat(formula = phh2o_igh_15_30 ~ 1, data = nmuestras)

Ahora que definimos la variable g1 podemos usar la función de predicción para interpolar. La función de predicción comprende:

-Un objeto raster-stars como un dem. -Un modelo-gstat como un g1.

El raster sirve para dos propositos:

-Especificar los lugares donde queremos hacer las predicciones en todos los métodos. -Especificar de valores de covariables (unicamente en Universal Kriging).

Ahora crearemos un objeto raster con valores de celda iguales a 1:

rrr = aggregate(geog.phh2o, 4)

llamamos al objeto rrr para ver sus características.

rrr
## class       : SpatRaster 
## size        : 258, 250, 1  (nrow, ncol, nlyr)
## resolution  : 0.008898822, 0.008898822  (x, y)
## extent      : -76.63117, -74.40647, 1.548893, 3.844789  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : phh2o_igh_15_30 
## min value   :        4.572673 
## max value   :        6.823683

Definiremos nuevos valores

values(rrr) <-1

Definiremos nuevos nombres

names(rrr) <- "valor"

Y observamos como cambio el objeto rrr

rrr
## class       : SpatRaster 
## size        : 258, 250, 1  (nrow, ncol, nlyr)
## resolution  : 0.008898822, 0.008898822  (x, y)
## extent      : -76.63117, -74.40647, 1.548893, 3.844789  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : valor 
## min value   :     1 
## max value   :     1
stars.rrr = st_as_stars(rrr)

Por ejemplo, la siguiente expresión interpola los valores del pH del agua en el departamento del huila según el modelo definido en g1 y la plantilla raster definida en stars.rrr

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

Llamamos al objeto z1 para observar sus características

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median     Mean  3rd Qu.     Max.  NA's
## var1.pred  4.619415 5.243824 5.372093 5.345633 5.453519 6.099254     0
## var1.var         NA       NA       NA      NaN       NA       NA 64500
## dimension(s):
##   from  to offset     delta                       refsys x/y
## x    1 250 -76.63  0.008899 +proj=longlat +datum=WGS8... [x]
## y    1 258  3.845 -0.008899 +proj=longlat +datum=WGS8... [y]

Para esta parte es importante tomar nota de los nombres de cada atributo del objeto z1. Podemos renombrar el primer atributo y cambiar su nombre a “phh2o”.

z1 = z1["var1.pred",,]
names(z1) = "phh2o"

Escogemos una paleta de color que más se acomode a lo que queremos visualizar.

paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), domain = 10:100, na.color = "transparent")

Ploteamos para visualizar la capa raster que interpolamos por medio del método IDW.

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 4.6:6.1)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$phh2o, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal=paleta, values= z1$phh2o,
    title = "pH Water"
    )
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
m

5.3 Interpolación por el método de Kriging Ordinario

Para hacer la interpolación con el método de kriging ordinario, es necesario entender que se requiere. El modelo de variograma, es una forma objetiva de cuantificar el patrón de autocorrelación de los datos y asignar pesos en consecuencia a la hora de hacer predicciones.

En primer lugar vamos a calcular y examinar el variograma empiríco utilizando la función de variograma.

La función tienes dos requerimientos -La formula, que precisa la variable dependiente y las covariables al igual que gstat. -Los datos, como una capa de puntos con la variable dependiente y las covariables como atributos de los puntos.

Por ejemplo, la siguiente expresión calcula el variograma empírico de las muestras, sin las covariables

v_emp_ok = variogram(phh2o_igh_15_30 ~ 1, data=nmuestras)

Ploteamos el variograma

plot(v_emp_ok)

Utilizaremos la función autofitVariogram de la libreria automap, para realizar un ajuste automático, de este modo podemos ajustar un modelo de variograma a un modelo de variograma empírico.

v_mod_ok = autofitVariogram(phh2o_igh_15_30 ~ 1, as(nmuestras, "Spatial"))

Esta función es muy útil ya que esta elige el modelo que más se ajusta ya ajusta sus parámetros, podemos utilizar show.vgms() para mostrar los tipos de modelo de variograma.

Tenemos que tener en cuenta que la función autofitVariogram no funciona en objetos sf, por lo que convertimos el objeto en un SpatialPointsDataFrame (paquete sp).

El modelo que más se ajuste puede ser gráficado con plot

plot(v_mod_ok)

El objeto que resulta es realmente una lista con varios componentes incluyendo el variograma empírico y el modelo del variograma ajustado, el componente $var_model contiene el modelo real.

v_mod_ok$var_model

La tabla muestra los parametros de los semivariogramas, ““Nug” hace referencia a Nugget (o pepita, en español) el cual es la variación a distancias muy pequeñas o error de medición

Ste puede hacer referencia al modelo Stein que modela la continuidad espacial a modelos mayores.

psill es la varianza parcial atribuida a cada modelo correspondiente.

range es la distancia apartir de la cual ya no hay distancia espacial significativa, para “Nug” es 0 porque se reprensenta la variación a distancia 0

kappa controla la suavidad del modelo

Una vez más, vamos a subconjuntar el atributo de valores predichos y cambiarle el nombre:

Ahora, podemos pasar el modelo del variograma a la función gstat, y podemos continuar con la interpolación Kriging Ordinario:

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

Una vez más, vamos a subconjuntar el atributo de valores predichos y cambiarle el nombre:

z2 = z2["var1.pred",,]
names(z2) = "phh2o"

y ploteamos para visualizar las predicciones del Kriging Ordinario

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z2,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 4.6:6.1)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$phh2o, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal = paleta, values= z2$phh2o,
    title = "OK pH del agua del Huila, interpolacion [%]"
    )
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
m

6. Evaluación de los resultados

6.1 Evaluación cualitativa

Podemos realizar un ploteo para obtener un mapa interactivo de nuestros datos de interés, y las interpolaciones anteriormente realizadas.

colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), domain = 4.6:6.1, na.color = "transparent")
m <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(stars.phh2o, opacity = 0.8, colorOptions = colores, group="RealWorld") %>%
  addGeoRaster(z1, colorOptions = colores, opacity = 0.8, group= "IDW")  %>%
  addGeoRaster(z2, colorOptions = colores, opacity = 0.8, group= "OK")  %>%
  # Add layers controls
  addLayersControl(
    overlayGroups = c("RealWorld", "IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
    addLegend("bottomright", pal = paleta, values= z1$phh2o,
    title = "pH Water"
)
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
m

6.2 Validación cruzada

Hemos estimado el pH del agua del departamento del Huila en zonas que no tienen datos por medio de dos métodos diferentes de interpolación, IDW y Kriging Ordinario, aunque resulta útil examinar y comparar los resultados de manera gráfica necesitamos una forma objetiva de evaluar la precisión de la interpolación y asi podremos elegir el método más preciso entre los diferentes métodos de interpolación que existen.

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

-Sacamos un punto de los datos de calibración -Hacemos una predicción para ese punto -Repetimos para todos los puntos Al final, lo que obtenemos es 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.

Comenta el siguiente código cuando vallas a correr el cuaderno, como se indica en el chunk, corre las lineas de código en la consola.

### corre los siguientes codigos en la consola linea por linea
cv1 = gstat.cv(g1)
cv2 = gstat.cv(g2)

El resultado de gstat.cv tiene los siguientes atributos:

-var1.pred-Valor predicho. -var1.var-Varianza (sólo para Kriging). -observed-Valor observado. -residual-Observed-Predicted. -zscore-Z-score (sólo para Kriging). -fold-Cross-validation ID.

 cv1 = na.omit(cv1)
 cv1
## class       : SpatialPointsDataFrame 
## features    : 780 
## extent      : -76.61894, -74.4187, 1.558904, 3.839228  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 6
## names       :        var1.pred, var1.var,         observed,           residual, zscore, fold 
## min values  :  4.8326780158247,       NA, 4.59999990463257, -0.483285165645715,     NA,    1 
## max values  : 6.06737366020752,       NA, 6.14971733093262,  0.593011435127221,     NA,  780

Convertimos el cv1 en objeto

cv1 = st_as_sf(cv1)

Y ploteamos los residuos

sp::bubble(as(cv1[, "residual"], "Spatial"))

A continuación, calculamos los índices de precisión de la predicción, como el error cuadrático medio (RMSE):

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

Ahora, repitamos el proceso con los resultados de Kriging Ordinario:

cv2 = na.omit(cv2)
cv2
## class       : SpatialPointsDataFrame 
## features    : 780 
## extent      : -76.61894, -74.4187, 1.558904, 3.839228  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 6
## names       :        var1.pred,           var1.var,         observed,           residual,            zscore, fold 
## min values  : 4.80403044256832, 0.0127572720524457, 4.59999990463257, -0.454440290141085, -3.27775585926781,    1 
## max values  : 5.94601854561173, 0.0265664863303328, 6.14971733093262,  0.533988015161604,  3.61934128053091,  780
cv2 = st_as_sf(cv2)
sp::bubble(as(cv2[, "residual"], "Spatial"))

Ahora calculamos los índices de precisión de la predicción nuevamente, el error cuadrático medio (RMSE):

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

En base a la evidencia y los resultados, resulta más confiable el Kriging Ordinario, analizando la validación cruzada, exactamente el bubble plot de residuos espaciales, el tamaño de las “burbujas” es más pequeño, lo cual indica errores más chicos, por otro lado en la evaluación cualitiva, se ve más definido y acorde al dato real el Kriging que el IDW.

Referencias

Lizarazo, I., 2023. Spatial interpolation of soil organic carbon. Available at https://rpubs.com/ials2un/soc_interp.

sessionInfo()
## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## 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] curl_6.2.3     dplyr_1.1.4    ggplot2_3.5.2  leafem_0.2.4   leaflet_2.2.2 
##  [6] automap_1.1-16 gstat_2.1-3    stars_0.6-8    abind_1.4-8    sf_1.0-21     
## [11] terra_1.8-54   sp_2.2-0      
## 
## loaded via a namespace (and not attached):
##  [1] generics_0.1.4     sass_0.4.10        class_7.3-23       KernSmooth_2.23-26
##  [5] lattice_0.22-6     digest_0.6.37      magrittr_2.0.3     evaluate_1.0.3    
##  [9] grid_4.5.0         RColorBrewer_1.1-3 fastmap_1.2.0      plyr_1.8.9        
## [13] jsonlite_2.0.0     e1071_1.7-16       reshape_0.8.10     DBI_1.2.3         
## [17] crosstalk_1.2.1    scales_1.4.0       codetools_0.2-20   jquerylib_0.1.4   
## [21] cli_3.6.5          rlang_1.1.6        units_0.8-7        intervals_0.15.5  
## [25] withr_3.0.2        base64enc_0.1-3    cachem_1.1.0       yaml_2.3.10       
## [29] FNN_1.1.4.1        raster_3.6-32      tools_4.5.0        parallel_4.5.0    
## [33] spacetime_1.3-3    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    classInt_0.4-11   
## [41] htmlwidgets_1.6.4  pkgconfig_2.0.3    pillar_1.10.2      bslib_0.9.0       
## [45] gtable_0.3.6       glue_1.8.0         Rcpp_1.0.14        tidyselect_1.2.1  
## [49] tibble_3.2.1       xfun_0.52          rstudioapi_0.17.1  knitr_1.50        
## [53] farver_2.1.2       htmltools_0.5.8.1  rmarkdown_2.29     xts_0.14.1        
## [57] compiler_4.5.0