1. Introducción

En este cuaderno se presentan dos técnicas de interpolación: Inverse Distance Weighted (IDW) y Ordinary Kriging (OK).

Inverse Distance Weighting (IDW)

Esta técnica determinística, también conocida como inverse distance-based weighted interpolation. Es la estimación del valor z en la locación x por una media ponderada de las observaciones ercanas.

Como se observa en la imagen anterior, los puntos rojos tienen valores de elevación conocidos. Así que los otros puntos serán interpolados:

El IDW es un método muy flexible de interpolación espacial. Se puede configurar de diferentes formas. Si se especifíca el radio de búsqueda, la interpolación solo usará el número de puntos conocidos dentro de ese radio de búsqueda.

Ordinary Kriging (OK)

El método probabilístico Ordinay Kriging es el kriging más utilizado. Se usa para estimar un valor en un punto de una región para la cual se conoce un variograma, utilizando datos de los vecinos de la ubicación de estimación.

La relación geoespacial (disimilitud entre puntos en función de la distancia/dirección) dentro del área de interés puede describirse mediante un variograma empírico. Se calculan las distancias entre los puntos de muestreo. Luego, se calculan los valores que conforman el variograma, ordenando las muestras cuyas distancias entre puntos caen dentro de un intervalo de distancia (h1, h2, etc.); calculando el valor como el promedio de la raíz cuadrada de la diferencia en las muestras γ(h).

Ambas técnicas se usarán para obtener una superficie plana del carbono orgánico del suelo de 15 a 30 cm de profundidad. Gracias a las muestras obtenidas de SoilGrids 250 m.

2. Configuración

Limpiamos el espacio de trabajo:

rm(list = ls())

Para la lectura y escritura de los datos geoespaciales, usaremos terra (para los ráster) y sf (para los vectores). Para la interpolación espacial, usaremos las librerías de gstat y automap. Finalmente, para la visualización, usaremos ggplot, leaflet y tmap.

Primero, necesitamos instalar las librerías requeridas desde la consola.

Una vez instaladas, vamos a cargar las librerías:

library(terra)
## terra 1.8.15
library(sf)
## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
library(sp)
library(stars)
## Cargando paquete requerido: abind
library(gstat)
library(automap)
library(RColorBrewer)
library(leaflet)
library(leafem)

3. Leer los datos de entrada

Nuestras muestras son valores de Carbono Orgánico del suelo entre 15-30 cm de profundidad. Fueron obtenidas desde SoilGrids 250 m usando este cuaderno.

Estas muestras se guardaron en formato de GeoPackage (GPKG), dentro del directorio de datos, usando soc_{depto}.gpkg como nombre del archivo.

Para empezar, validemos que los datos estén disponibles para usarse:

list.files(path="D:/GB2_UNAL/Proyecto_5/INFORME_2", pattern = "*.gpkg")
## [1] "Mun_meta.gpkg"  "soc_stder.gpkg"

Ahora, leamos las muestras usando sf.

samples = sf::st_read("D:/GB2_UNAL/Proyecto_5/INFORME_2/soc_stder.gpkg")
## Reading layer `soc_stder' from data source 
##   `D:\GB2_UNAL\Proyecto_5\INFORME_2\soc_stder.gpkg' using driver `GPKG'
## Simple feature collection with 1937 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.93558 ymin: 1.619848 xmax: -71.02165 ymax: 4.894763
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

Como nuestra área de estudio es El Meta,es necesario leer también el GeoPackage correspondiente de los municipios:

munic = sf::st_read("D:/GB2_UNAL/Proyecto_5/INFORME_2/mun_meta.gpkg")
## Reading layer `municipios_colombia' from data source 
##   `D:\GB2_UNAL\Proyecto_5\INFORME_2\mun_meta.gpkg' using driver `GPKG'
## Simple feature collection with 29 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.93335 ymin: 1.617732 xmax: -71.07753 ymax: 4.899101
## Geodetic CRS:  MAGNA-SIRGAS

Como mi gpkg de municipios tiene el sistema de coordenadas MAGNA-SIRGAS, lo transformarémos a WGS 84:

munic_wgs84 = st_transform(munic, crs=4326)

4. Explorar los datos de entrada

Vamos generar un resumen de las muestras:

summary(samples)
##       soc                    geom     
##  Min.   : 8.449   POINT        :1937  
##  1st Qu.:12.170   epsg:NA      :   0  
##  Median :14.803   +proj=long...:   0  
##  Mean   :20.642                       
##  3rd Qu.:20.841                       
##  Max.   :86.841

En total tenemos 1937 muestras de Carbono Orgánico del suelo.

Representación gráfica con un histograma:

hist(samples$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)

El hitograma anterior presenta la frecuncia de g/kG de carbono orgánico en el suelo, observamos que se presenta mayor frecuencia en 10 - 15 g/kg.

Vamos a redondear los valores de Carbono Orgánico del suelo a dos dígitos:

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

Visualicemos las muestras para saber su distribución espacial, usando leaflet.

pal <- colorNumeric(
  c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),

  domain = samples$soc,
  )
leaflet() %>%
  addPolygons(
    data = munic_wgs84,
    color = "gray",
    # set the opacity of the outline
    opacity = 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 = ~soc,
    position = "bottomleft",
    title = "SOC:",
    opacity = 0.7) %>%
  addProviderTiles("OpenStreetMap")

Con el color de los marcadores podemos obtener el valor del contenido de Carbono Orgánico del suelo según dónde ubiquemos el puntero.

Según esta visualización, la mayor cantidad está en la zona de la Cordillera Oriental en los Parques Nacionales Naturales Sumapaz y Cordillera de los Picachos.

5. Interpolación

5.1 Creación del objeto gstat

Para interpolar, necesitams primero crear un objeto de clase gstat, usando la función con el mismo nombre gstat.

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

  • La definición del modelo
  • Los datos de calibración

Basándonos en sus argumentos, la función gstats “entiende” qué modelo de interpolación queremos usar:

  • Modelo sin variograma → Inverse Distance Weighting (IDW)

  • Modelo de variograma, sin covariables → Ordinary Kriging (OK)

Usaremos 3 parámetros de la función gstat:

- Fómula: La predicción “fórmula” especifíca las variables dependientes e intependientes (también conocidas como covariables)

- Datos: Los datos de calibración (también conocidos como datos de entrenamiento)

- Modelo: El modelo de variograma

5.2 Interpolación IDW

Para interpolar el contenido de Carbono Orgánico usando el método IDW, creamos el siguiente objeto gstat especificando solo la fórmula Y los datos:

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

Ya que nuestro modelo de interpolación g1 está definido, podemos usar la función predecir para interpolar realmente, es decir, para estimar los valores de carbono orgánico del suelo.

La función predecir acepta:

  • Un objeto raster-stats, como el dem

  • Un objeto model-gstat, como el g1

El ráster sirve para dos propósitos:

  • Especificar las ubicaciones donde queremos hacer las predicciones (en todos los métodos)

  • Especificar los valores de las covariables (Solo en Kriging Universal)

Utilicemos la biblioteca 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.

(rrr = terra::rast(xmin=-74.93335, xmax=-71.07753, ymin=1.617732 , ymax= 4.899101, nrows=370, ncols=329, vals= 1, crs= "epsg:4326"))
## class       : SpatRaster 
## dimensions  : 370, 329, 1  (nrow, ncol, nlyr)
## resolution  : 0.01171982, 0.008868565  (x, y)
## extent      : -74.93335, -71.07753, 1.617732, 4.899101  (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 valores de carbono orgánico del suelo de acuerdo al modelo definido en g1 y en el ráster definido en stars.rrr:

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

z1

Son los valores interpolados de carbono orgánico del suelo mediante el modelo IDW y reflejados en el ráster del departamento del Meta.

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min. 1st Qu.   Median     Mean  3rd Qu.     Max.   NA's
## var1.pred  8.474784   13.79 16.22458 20.96207 22.05424 83.40085      0
## var1.var         NA      NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -74.93   0.01172 WGS 84 [x]
## y    1 370  4.899 -0.008869 WGS 84 [y]

Los atributos del objeto stars presentan las variables o medidas asociadas con el objeto, en este caso, hay dos atributos:

1. var1.pred:

Este atributo contiene valores numéricos que representan las predicciones o estimaciones.

2. var1.var:

Este atributo representa la varianza asociada con las predicciones.

Podemos subconjuntar el primer atributo y cambiarle el nombre:

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

Este chunk reemplaza todos los valores NA en el primer atributo de Z1 y los reemplaza por 0.0. El resto de los valores en el primer atributo permanecen sin cambios.

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min. 1st Qu.   Median     Mean  3rd Qu.     Max.   NA's
## var1.pred  8.474784   13.79 16.22458 20.96207 22.05424 83.40085      0
## var1.var         NA      NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -74.93   0.01172 WGS 84 [x]
## y    1 370  4.899 -0.008869 WGS 84 [y]

Les damos nombre a los valores de la interpolación:

names(z1) = "soc"

Para trazar la superficie interpolada, necesitamos una paleta de colores:

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

Ahora, visualicemos el resultado de la interpolación.

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 [%]"
    )
m

5.3 Interpolación OK

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

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

Esta función requiere dos argumentos:

  • Fórmula: Especifíca la variable dependiente y las cobariables, justo como en gstat

  • Datos: 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 las muestras, sin covariables:

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

Tracemos el variograma empírico:

plot(v_emp_ok)

Hay muchas formas de ajustar un modelo de variograma a un variograma empírico. Aquí utilizaremios la más sencilla: el ajuste automático mediante la función autofitVariogram del paquete automap:

v_mod_ok = automap::autofitVariogram(soc ~ 1, as(samples, "Spatial"))

La función elige el tipo de modelo que mejor se ajusta y también ajusta finamente sus parámetros.

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

El modelo ajustado puede ser trazado con 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 propiamente dicho.

v_mod_ok$var_model
##   model      psill    range
## 1   Nug   9.667854   0.0000
## 2   Exp 268.740402 147.1198

Ahora, el modelo de variograma se puede pasar 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, subconjuntaremos el atributo de valores predichos y le cambiaremos el nombre.

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

Las predicciones de Kriging Ordinario (OK) se muestran en el siguiente gráfico.

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 [%]"
    )
m

6. Evaluación de los resultados

6.1 Evaluación cualitativa

Finalmente, vamos a comparar los dos métodos entre sí:

colores <- colorOptions(palette = c("orange", "yellow","cyan", "green"), domain = 10:100, na.color = "transparent")
m <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(z1, colorOptions = colores, opacity = 0.75, group= "IDW")  %>%
  addGeoRaster(z2, colorOptions = colores, opacity = 0.75, 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

6.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 seleccionar el método más preciso entre las opciones disponibles.

En Leave-One-Out Cross Validation (LOOCV):

  • Se extrae 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 valores observados y valores predichos para todos los puntos.

Podemos ejecutar LOOCV usando la función gstat.cv, que acepta un objeto de gstat.

El resultado de *g,stat.cv tiene los siguientes atributos:

  • var1.pred: Valor predicho
  • var1.var. Varianza (solo para el Kriging)
  • Observados: Valores observados
  • Resuduales: Obvervados-Predecidos
  • zscore: Z-score (solo para el Kriging)
  • Fold: ID de validación cruzada”

Valores residuales del OK

Valor RMSE del OK: 5.262043

Valores residuales del IDW

Valor REMSE de IDW: 6.384484

6.3 ¿Qué es el RMSE?

El error cuadrático medio (RMSE) mide la cantidad de error que hay entre dos conjuntos de datos. En otras palabras, compara un valor predicho y un valor observado o conocido.

Los valores más bajos de RMSE indican un mejor ajuste. RMSE es una buena medida de la precisión con que el modelo predice la respuesta, y es el criterio más importante para ajustar si el propósito principal del modelo es la predicción.

7. Interpretación del RMSE

En nuestro caso:

El RMSE, al ser una medida del error cuadrático medio, indica qué tan bien cada modelo predice los valores en la región de estudio. Un valor más bajo sugiere un mejor ajuste de los datos, por lo que el modelo OK presentó un desempeño superior en comparación con el IDW. Esto se debe a que el Kriging Ordinario no solo pondera las distancias entre los puntos de muestreo, sino que también incorpora la estructura espacial y la correlación entre ellos a través del semivariograma. En contraste, el IDW asigna pesos a los puntos cercanos sin modelar explícitamente la variabilidad espacial, lo que puede llevar a predicciones menos precisas.

Dado que el estudio se basó en 1937 muestras, es probable que exista una estructura espacial significativa en la distribución del carbono orgánico del suelo, lo que permitió que el Kriging captara mejor la variabilidad y redujera el error de predicción. Por esta razón, el modelo OK se considera más adecuado para la interpolación en la región del departamento del Meta.

8. Guardamos los resultados de interpolación

write_stars(
  z1, dsn = "D:/GB2_UNAL/Proyecto_5/INFORME_2/Modelos de interpolación/IDW_soc_stder.tif", layer = 1
)
write_stars(
  z1, dsn = "D:/GB2_UNAL/Proyecto_5/INFORME_2/Modelos de interpolación/OK_soc_stder.tif", layer = 1
)

Referencias

sessionInfo()
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## 
## 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] leafem_0.2.3       leaflet_2.2.2      RColorBrewer_1.1-3 automap_1.1-16    
##  [5] gstat_2.1-3        stars_0.6-8        abind_1.4-8        sp_2.1-4          
##  [9] sf_1.0-19          terra_1.8-15      
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6            xfun_0.49               bslib_0.8.0            
##  [4] ggplot2_3.5.1           raster_3.6-31           htmlwidgets_1.6.4      
##  [7] lattice_0.22-6          leaflet.providers_2.0.0 vctrs_0.6.5            
## [10] tools_4.4.2             crosstalk_1.2.1         generics_0.1.3         
## [13] parallel_4.4.2          tibble_3.2.1            proxy_0.4-27           
## [16] spacetime_1.3-3         fansi_1.0.6             xts_0.14.1             
## [19] pkgconfig_2.0.3         KernSmooth_2.23-24      lifecycle_1.0.4        
## [22] compiler_4.4.2          farver_2.1.2            FNN_1.1.4.1            
## [25] munsell_0.5.1           codetools_0.2-20        htmltools_0.5.8.1      
## [28] class_7.3-22            sass_0.4.9              yaml_2.3.10            
## [31] pillar_1.9.0            jquerylib_0.1.4         classInt_0.4-10        
## [34] cachem_1.1.0            tidyselect_1.2.1        digest_0.6.37          
## [37] dplyr_1.1.4             fastmap_1.2.0           grid_4.4.2             
## [40] colorspace_2.1-1        cli_3.6.3               magrittr_2.0.3         
## [43] base64enc_0.1-3         utf8_1.2.4              e1071_1.7-16           
## [46] scales_1.3.0            rmarkdown_2.29          zoo_1.8-13             
## [49] png_0.1-8               evaluate_1.0.1          knitr_1.49             
## [52] rlang_1.1.4             Rcpp_1.0.13             glue_1.8.0             
## [55] DBI_1.2.3               rstudioapi_0.17.1       reshape_0.8.9          
## [58] jsonlite_1.8.9          R6_2.5.1                plyr_1.8.9             
## [61] intervals_0.15.5        units_0.8-5

Hartmann, K., Krois, J., Rudolph, A. (2023): Statistics and Geodata Analysis using R (SOGA-R). Department of Earth Sciences, Freie Universitaet Berlin.

Wackernagel, H. (1995). Ordinary Kriging. In: Multivariate Geostatistics. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-662-03098-1_11