1. Introduccion

Este cuaderno ilustra dos técnicas de interpolación espacial:Inverse Distance Weighted (IDW) y Ordinary Kriging, (OK)IDW es una técnica determinista. OK es probabilístico. Ambas técnicas se utilizan aquí para obtener una superficie continua de SOC a 15-30 cm de muestras obtenidas de SoilGrids 250 m del departamento de Antioquia.

2. Configuracion

Primero limpiamos la memoria

rm(list=ls())

Luego cargamos las librerias

library(sp)
## Warning: package 'sp' was built under R version 4.2.3
library(terra)
## Warning: package 'terra' was built under R version 4.2.3
library(sf)
## Warning: package 'sf' was built under R version 4.2.3
library(stars)
## Warning: package 'stars' was built under R version 4.2.3
library(gstat)
## Warning: package 'gstat' was built under R version 4.2.3
library(automap)
## Warning: package 'automap' was built under R version 4.2.3
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.2.3
library(leafem)
## Warning: package 'leafem' was built under R version 4.2.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.2
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.2.2
library(curl)
## Warning: package 'curl' was built under R version 4.2.3
h <- new_handle()
handle_setopt(h, http_version = 2)

3. Lectura de los datos de entrada

leamos la capa SOC que descargamos de ISRIC) usando la biblioteca terra:

archivo <- ("G:/GB_2/interp_SOC/soc_ant.tif")
(soc <- rast(archivo))
## class       : SpatRaster 
## dimensions  : 2511, 1770, 1  (nrow, ncol, nlyr)
## resolution  : 0.002259887, 0.002389486  (x, y)
## extent      : -77.1483, -73.1483, 4.5209, 10.5209  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : soc_ant.tif 
## name        : soc_ant

Ahora convertimos los datos SOC en porcentaje y revisamos el factor de escala de SOC en SoidGrids:

soc.perc <-  soc/10

Luego transformamos de CRS en el conocido WGS84 CRS:

geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(soc.perc, geog))
## class       : SpatRaster 
## dimensions  : 2556, 1704, 1  (nrow, ncol, nlyr)
## resolution  : 0.00234726, 0.00234726  (x, y)
## extent      : -77.1483, -73.14857, 4.521303, 10.5209  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        :  soc_ant 
## min value   :   0.0000 
## max value   : 273.5996

Convertimos la capa SpatRaster en un stars object:

ngeog.soc = aggregate(geog.soc, 4)
stars.soc = st_as_stars(ngeog.soc)

Observamos el mapa:

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

4. Muestreo del mundo

Cogemos una muestra de aprox. 500 sitios de los datos del mundo real utilizando una muestra ubicada aleatoriamente:

set.seed(123456)

Muestreo aleatorio de 500 points

(samples <- spatSample(geog.soc, 500, "random", as.points=TRUE))
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 500, 1  (geometries, attributes)
##  extent      : -77.14713, -73.14974, 4.527172, 10.51268  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : soc_ant
##  type        :   <num>
##  values      :       0
##                  81.98
##                      0

convertimos el spatVector object en simple feature object para poder describir las características principales del objeto de muestras:

(muestras <- sf::st_as_sf(samples))
## Simple feature collection with 500 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -77.14713 ymin: 4.527172 xmax: -73.14974 ymax: 10.51268
## 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_ant                   geometry
## 1    0.00000 POINT (-73.84218 9.972815)
## 2   81.97841 POINT (-75.41485 6.050543)
## 3    0.00000 POINT (-75.03929 8.853172)
## 4   68.39283 POINT (-73.89617 5.287684)
## 5  227.27423 POINT (-76.62134 5.707843)
## 6    0.00000 POINT (-75.89134 9.620726)
## 7   65.94017 POINT (-76.14015 7.296938)
## 8    0.00000  POINT (-73.56521 8.96584)
## 9   37.74040 POINT (-74.51115 6.336909)
## 10  52.52910 POINT (-74.10977 7.728834)

Eliminamos las filas que tienen valores faltantes

nmuestras <- na.omit(muestras)

Visualizamos las muestras

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

creamos un vector llamado “id” que contiene una secuencia de números enteros desde 1 hasta 500

id <- seq(1,500,1)

Calculamos los sitios

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

Eliminamos los valores de NA:

sitios <- na.omit(sitios)

Observamos el resultado

Resumen de sitios

head(sitios)
##   id    longit    latit       soc
## 1  1 -73.84218 9.972815   0.00000
## 2  2 -75.41485 6.050543  81.97841
## 3  3 -75.03929 8.853172   0.00000
## 4  4 -73.89617 5.287684  68.39283
## 5  5 -76.62134 5.707843 227.27423
## 6  6 -75.89134 9.620726   0.00000

Resumen de sitios

tail(sitios)
##      id    longit    latit      soc
## 495 495 -74.05344 8.357900 61.91631
## 496 496 -73.68727 5.888582 84.75270
## 497 497 -75.88899 6.811056 88.13267
## 498 498 -74.40083 9.709922  0.00000
## 499 499 -75.15665 6.038807 59.35070
## 500 500 -75.25054 6.172601 98.34066

Ahora visualizamos las muestras

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

Ahora, estamos listos para realizar las tareas de interpolación:

5. Interpolación

5.1 Creación del gstat object

Para interpolar, primero necesitamos crear un objeto de clase gstat, usando una función del mismo nombre: gstat.

fórmula: a.k.a covariates data: a.k.a the train data

5.2 IDW interpolacion

creamos el siguiente objeto gstat, especificando solo la fórmula y los datos:

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

Ahora que nuestro modelo de interpolación g1 está definido, podemos usar la función actually interpolate para estimar los valores de precipitación:

creamos un objeto ráster con valores de celda iguales a 1:

rrr = aggregate(ngeog.soc, 4)

Visualizamos:

rrr
## class       : SpatRaster 
## dimensions  : 160, 107, 1  (nrow, ncol, nlyr)
## resolution  : 0.03755616, 0.03755616  (x, y)
## extent      : -77.1483, -73.12979, 4.511914, 10.5209  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_ant 
## min value   :   0.000 
## max value   : 238.116

Definimos nuevos valores:

values(rrr) <- 1

Definimos nuevos nombres:

names(rrr) <- "valor"

Visualizamos el resultado:

rrr
## class       : SpatRaster 
## dimensions  : 160, 107, 1  (nrow, ncol, nlyr)
## resolution  : 0.03755616, 0.03755616  (x, y)
## extent      : -77.1483, -73.12979, 4.511914, 10.5209  (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)

la siguiente expresión interpola los valores SOC según el modelo definido en g1 y la plantilla ráster definida en stars.rrr:

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

Observamos 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.01220758 19.34815 53.31169 56.52249 81.47947 218.9322     0
## var1.var           NA       NA       NA      NaN       NA       NA 17120
## dimension(s):
##   from  to   offset      delta                       refsys x/y
## x    1 107 -77.1483  0.0375562 +proj=longlat +datum=WGS8... [x]
## y    1 160  10.5209 -0.0375562 +proj=longlat +datum=WGS8... [y]

creamos un subconjunto solo del primer atributo y cambiarle el nombre a “soc”:

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

Configuramos paleta de colores:

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

El ráster SOC interpolado, usando IDW, se muestra en la siguiente figura:

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = 11:55)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$soc, clusterOptions = markerClusterOptions()) %>%
    addLegend("bottomright", pal=paleta, values= z1$soc,
    title = "IDW SOC interpolacion de antioquia [%]"
    )
## 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 OK interpolacion

Los métodos 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.

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

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

obtenemos el variograma:

plot(v_emp_ok)

Ajustamos el modelo de variograma usando la función autofitVariogram del paquete automap:

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

El modelo ajustado se puede trazar con plot:

plot(v_mod_ok)

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

v_mod_ok$var_model
##   model       psill    range kappa
## 1   Nug    159.0497     0.00   0.0
## 2   Ste 140414.9588 84191.34   0.4
  • El modelo NUG (Nugget) y STE (Sill-To-Effective) se utiliza en geostadística para analizar la variabilidad espacial de un fenómeno. El NUG representa la variabilidad mínima o los errores aleatorios en los datos, mientras que el STE describe la variabilidad estructural a medida que aumenta la distancia entre los puntos. Juntos, proporcionan una descripción completa de la variabilidad espacial de los datos, permitiendo ajustar el variograma y realizar interpolaciones para predecir los valores en ubicaciones no muestreadas.

Ahora, el modelo de variograma se puede pasar a la función gstat, y podemos continuar con the Ordinary Kriging interpolation:

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

Nuevamente, subdividiremos el atributo de valores predichos y lo renombraremos:

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

Las Ordinary Kriging predictions se muestran en la siguiente figura:

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

6 Evaluacion de resultados

6.1 Evaluación cualitativa

Otra vista de las tres salidas de interpolación:

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

Observamos el resultado:

m <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(stars.soc, 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")  %>%

  addLayersControl(
    overlayGroups = c("RealWorld", "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
m

6.2 Cross-validation

podemos elegir objetivamente el método más preciso entre los métodos de interpolación disponibles, a continuacion escribiremos el siguiente fragmento, ocultando el mensaje y los resultados:

cv1 = gstat.cv(g1)

El resultado de gstat.cv tiene los siguientes atributos:

var1.pred—Predicted value var1.var—Variance (only for Kriging) observed—Observed value residual—Observed-Predicted zscore—Z-score (only for Kriging) fold—Cross-validation ID

cv1 = na.omit(cv1)

Observamos el resultado

cv1
## class       : SpatialPointsDataFrame 
## features    : 500 
## extent      : -77.14713, -73.14974, 4.527172, 10.51268  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 6
## names       :         var1.pred, var1.var,         observed,          residual, zscore, fold 
## min values  : 0.275784606749383,       NA,                0, -75.7345820607976,     NA,    1 
## max values  :  193.753858321137,       NA, 227.274230957031,  103.004855858435,     NA,  500

Convirtamos el objeto cv1:

cv1 = st_as_sf(cv1)

Ahora, grafiquemos los residuos:

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

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

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

calculamos el cv2

cv2 = gstat.cv(g2)

Tiempo de conversión:

cv2 = st_as_sf(cv2)

Calcule RSME para obtener resultados correctos:

sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))
## [1] 19.32786
  • La técnica de interpolación IDW es especialmente adecuada para observar los datos de carbono orgánico en el suelo en el departamento de Antioquia. Con su topografía variada y terreno montañoso, Antioquia presenta una distribución espacial irregular de las mediciones del carbono orgánico. El IDW asigna un mayor peso a las mediciones cercanas al punto de interpolación, lo que refleja la influencia de la proximidad geográfica en los valores estimados. Además, teniendo en cuenta los datos obtenidos en SoilGrids, brinda una mayor cantidad de información para la interpolación. La técnica IDW es fácil de implementar y computacionalmente eficiente, lo que la hace práctica a gran escala.El IDW es una opción idónea para obtener una representación precisa y detallada del carbono orgánico en el suelo de Antioquia, lo que resulta fundamental para comprender su distribución y su impacto en la salud y calidad del suelo.

Referencias

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

sessionInfo()
## R version 4.2.1 (2022-06-23 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252   
## [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C                  
## [5] LC_TIME=Spanish_Spain.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] curl_5.0.1    dplyr_1.1.0   ggplot2_3.4.1 leafem_0.2.0  leaflet_2.1.2
##  [6] automap_1.1-9 gstat_2.1-1   stars_0.6-1   abind_1.4-5   sf_1.0-12    
## [11] terra_1.7-23  sp_1.6-0     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.10        lattice_0.20-45    FNN_1.1.3.2        png_0.1-8         
##  [5] class_7.3-20       zoo_1.8-11         digest_0.6.31      utf8_1.2.3        
##  [9] R6_2.5.1           plyr_1.8.8         evaluate_0.20      e1071_1.7-13      
## [13] highr_0.10         pillar_1.8.1       rlang_1.1.0        rstudioapi_0.14   
## [17] raster_3.6-20      jquerylib_0.1.4    rmarkdown_2.20     rgdal_1.6-5       
## [21] htmlwidgets_1.6.2  munsell_0.5.0      proxy_0.4-27       compiler_4.2.1    
## [25] xfun_0.37          pkgconfig_2.0.3    base64enc_0.1-3    htmltools_0.5.4   
## [29] tidyselect_1.2.0   tibble_3.2.1       intervals_0.15.3   codetools_0.2-18  
## [33] reshape_0.8.9      fansi_1.0.4        spacetime_1.3-0    withr_2.5.0       
## [37] grid_4.2.1         jsonlite_1.8.4     lwgeom_0.2-11      gtable_0.3.1      
## [41] lifecycle_1.0.3    DBI_1.1.3          magrittr_2.0.3     units_0.8-1       
## [45] scales_1.2.1       KernSmooth_2.23-20 cli_3.6.0          cachem_1.0.7      
## [49] farver_2.1.1       bslib_0.4.2        ellipsis_0.3.2     xts_0.13.1        
## [53] generics_0.1.3     vctrs_0.6.1        tools_4.2.1        glue_1.6.2        
## [57] crosstalk_1.2.0    parallel_4.2.1     fastmap_1.1.1      yaml_2.3.7        
## [61] colorspace_2.1-0   classInt_0.4-9     knitr_1.42         sass_0.4.5