Interpolación espacial del carbono orgánico del suelo

1. INTRODUCCIÓN

Este cuaderno ilustra dos técnicas de interpolación espacial: distancia ponderada inversa (IDW) y Kriging ordinario (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 a partir de muestras obtenidas de SoilGrids 250 m.

2. CONFIGURACIÓN

-Se limpia la memoria

rm(list=ls())

-Ahora se instalan las librerias:

library(sp)

## Warning: package 'sp' was built under R version 4.3.2

library(terra)

## Warning: package 'terra' was built under R version 4.3.2

## terra 1.7.55

library(sf)

## Warning: package 'sf' was built under R version 4.3.2

## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE

library(stars)

## Warning: package 'stars' was built under R version 4.3.2

## Loading required package: abind

library(gstat)

## Warning: package 'gstat' was built under R version 4.3.2

library(automap)

## Warning: package 'automap' was built under R version 4.3.2

library(leaflet)

## Warning: package 'leaflet' was built under R version 4.3.2

library(leafem)

## Warning: package 'leafem' was built under R version 4.3.2

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.3.2

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:terra':
## 
##     intersect, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(curl)

## Warning: package 'curl' was built under R version 4.3.2

## Using libcurl 8.3.0 with Schannel

h <- new_handle()
handle_setopt(h, http_version = 2)

3. LEER LOS DATOS DE ENTRADA

Necesitamos leer un conjunto de datos para imitar datos del mundo real. Por lo tanto, leamos la capa SOC que descargamos de (ISRIC) usando la biblioteca terra:

archivo <- ("soc_igh_15_30.tif")
(soc <- rast(archivo))

## class       : SpatRaster 
## dimensions  : 1090, 712, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8476000, -8298000, 319750, 592250  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine 
## source      : soc_igh_15_30.tif 
## name        : soc_igh_15_30

Ahora, se convierten los datos del SOC en porcentaje.

soc.perc <- soc/10

-¿Cuál es el CRS de los datos del mundo real?

Parece que necesitamos una transformación de dicho CRS al conocido CRS WGS84:

geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(soc.perc, geog))

## class       : SpatRaster 
## dimensions  : 1104, 758, 1  (nrow, ncol, nlyr)
## resolution  : 0.002216352, 0.002216352  (x, y)
## extent      : -76.11119, -74.4312, 2.87342, 5.320272  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_igh_15_30 
## min value   :      11.02765 
## max value   :     120.17535

Convirtamos la capa SpatRaster en un objeto stars:

stars.soc = st_as_stars(geog.soc)

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

4. MUESTREO DEL MUNDO

-Consigamos una muestra de aprox. 500 sitios a partir de datos del mundo real utilizando una muestra ubicada aleatoriamente:

set.seed(123456)

# Random sampling of 500 points
(samples <- spatSample(geog.soc, 500, "random", as.points=TRUE))

##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 500, 1  (geometries, attributes)
##  extent      : -76.11008, -74.4323, 2.876745, 5.314731  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : soc_igh_15_30
##  type        :         <num>
##  values      :         68.33
##                        52.37
##                        28.41

-Describa las características principales del objeto de muestra.

• Tipo: Esto indica que los valores son de tipo numérico.

• Valores: Estos son los valores numéricos específicos que resultaron del muestreo aleatorio. En este caso, los valores son: 68.33, 52.37 y 28.41.

Ahora, necesitamos convertir el objeto spatVector en un objeto de característica simple:

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

nmuestras <- na.omit(muestras)

-Ahora visualicemos las muestras:

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

id <- seq(1,500,1) 

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

-Ahora removemos los elementos con valores de NA, lo que significa que no existen.

sitios <- na.omit(sitios)

head(sitios)

Ahora visualizamos las muestras por medio de un mapa:

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

-Ahora estamos listos para realizar las tareas de interpolación.

5. INTERPOLACIÓN

5.1 Creación del objeto gstat

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

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

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

Según sus argumentos, la función gstat “entiende” qué tipo de modelo de interpolación queremos usar:

• Sin modelo de variograma → IDW
• Modelo de variograma, sin covariables → Kriging ordinario

Vamos a utilizar tres parámetros de la función gstat:

• Fórmula: la “fórmula” de predicción que especifica las variables dependientes e independientes (también conocidas como covariables)
• Datos: los datos de calibración (también conocidos como los datos del tren)
• Modelo: el modelo de variograma

5.2 interpolación IDW

Para interpolar SOC usando el método IDW creamos el siguiente objeto gstat, especificando solo la fórmula y los datos:

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

Ahora que nuestro modelo de interpolación g1 está definido, podemos usar la función de predicción para interpolar, es decir, estimar los valores de precipitación.

La función de predicción acepta:

• Un objeto ráster: estrellas, como dem
• Un modelo: objeto gstat, como g1

El ráster tiene dos propósitos:

• Especificar las ubicaciones donde queremos hacer predicciones (en todos los métodos)
• Especificación de valores de covariables (solo en Universal Kriging)

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

# Una copia sencilla
rrr = aggregate(geog.soc, 4)

¿Qué es rrr?

rrr

## class       : SpatRaster 
## dimensions  : 276, 190, 1  (nrow, ncol, nlyr)
## resolution  : 0.008865406, 0.008865406  (x, y)
## extent      : -76.11119, -74.42676, 2.87342, 5.320272  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : soc_igh_15_30 
## min value   :      11.97661 
## max value   :     103.90206

Definir nuevos valores:

values(rrr) <-1

Definir nuevos nombres:

names(rrr) <- "valor"

¿Que es rrr ahora?

rrr

## class       : SpatRaster 
## dimensions  : 276, 190, 1  (nrow, ncol, nlyr)
## resolution  : 0.008865406, 0.008865406  (x, y)
## extent      : -76.11119, -74.42676, 2.87342, 5.320272  (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 SOC según el modelo definido en g1 y la plantilla ráster definida en stars.rrr:

## [interpolación ponderada por distancia inversa]
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  14.00324 32.48897 41.49431 42.04787 50.96143 95.0423     0
## var1.var         NA       NA       NA      NaN       NA      NA 52440
## dimension(s):
##   from  to offset     delta                       refsys x/y
## x    1 190 -76.11  0.008865 +proj=longlat +datum=WGS8... [x]
## y    1 276   5.32 -0.008865 +proj=longlat +datum=WGS8... [y]

Tome nota de los nombres de los dos atributos incluidos dentro del objeto z1.

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

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

Necesitamos una paleta de color:

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

El ráster SOC interpolado, utilizando 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 interpolation [%]"
    )
m  # Print the map

5.3 OK interpolación

Los 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 ponderaciones en consecuencia al hacer predicciones.

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

La función requiere dos argumentos:

• Fórmula: especifica la variable dependiente y las covariables, como en gstat
• Datos: la capa de puntos con la variable dependiente y covariables como atributos de puntos Por ejemplo, la siguiente expresión calcula el variograma empírico de muestras, sin covariables:

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

Ahora se plotea el variograma:

plot(v_emp_ok)

Hay varias formas de ajustar un modelo de variograma a un variograma empírico. Usaremos el más simple: ajuste automático usando la función autofitVariogram del paquete automap:

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

La función elige el tipo de modelo que mejor se adapta y también ajusta los parámetros. Puede utilizar show.vgms() para mostrar los tipos de modelos de variograma.

Tenga 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 ajustado se puede trazar con plot:

plot(v_mod_ok)

El objeto resultante es en realidad una lista con varios componentes, incluido 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

En tu cuaderno, explica el significado de cada elemento en el modelo anterior.

• NUG: El nugget representa variabilidad a una escala más fina que la resolución de tus datos o a variabilidad debida a errores de medición. En tu caso, el valor del nugget es 31.53533.

• PSILL: El sill es el valor límite que se alcanza a medida que la distancia entre los puntos aumenta indefinidamente. Representa la máxima variabilidad espacial. En tu caso, el valor del sill es 361.35783.

• RANGE: El rango es la distancia a la cual se alcanza el sill o un porcentaje significativo del sill. Representa la “longitud” de la correlación espacial. En tu caso, el rango es de 0.00000 a 94.09982 unidades de distancia.

Estos componentes son fundamentales para modelar el variograma y, por ende, para realizar interpolaciones espaciales o predicciones en lugares no muestreados. Dependiendo del comportamiento de estos componentes, puedes elegir diferentes modelos de variograma que se ajusten mejor a la estructura espacial de tus datos.

-Ahora, el modelo de variograma se puede pasar a la función gstat y podemos continuar con la interpolación Kriging ordinaria:

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

## [using ordinary kriging]

Nuevamente, crearemos un subconjunto del atributo de valores predichos y le cambiaremos el nombre:

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

Las predicciones de Kriging ordinario 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 interpolation [%]"
    )
m  # Imprimir el mapa

6. EVALUACIÓN 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")

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")  %>%
  # Add layers controls
  addLayersControl(
    overlayGroups = c("RealWorld", "IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
    addLegend("bottomright", pal = paleta, values= z1$soc,
    title = "Soil organic carbon [%]"
)
m  # Imprimir el mapa

6.2 Validación cruzada

Hemos estimado las superficies climáticas 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 esa 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 nosotros:

• Sacamos un punto de los datos de calibración
• Se hizo una predicción para ese punto
• Repita para todos los puntos. Al final, lo que obtenemos es una tabla con un valor observado y un valor predicho para todos los puntos
• Se puede ejecutar la validación cruzada Leave-One-Out usando la función gstat.cv, que acepta un objeto gstat

### run the following code from the console -- line-by-line
cv1 = gstat.cv(g1)
#cv2 = gstat.cv(g2)

El resultado de gstat.cv tiene los siguientes atributos:

• var1.pred: valor previsto
• var1.var: varianza (solo para Kriging)
• observado: valor observado
• residual: observado-predicho
• zscore: puntuación Z (solo para Kriging)
• plegar: ID de validación cruzada

cv1 = na.omit(cv1)

cv1

## class       : SpatialPointsDataFrame 
## features    : 477 
## extent      : -76.0857, -74.4722, 2.876745, 5.314731  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 6
## names       :        var1.pred, var1.var,         observed,          residual, zscore, fold 
## min values  : 17.2291317140456,       NA,  13.805643081665, -23.1561670705783,     NA,    1 
## max values  : 74.7313777581187,       NA, 95.9597930908203,  47.1583280514459,     NA,  477

Ahora convertimos el objeto cv1:

cv1 = st_as_sf(cv1)

Ahora, grafiquemos los residuos:

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

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

# This is the RMSE value for the IDW interpolation
sqrt(sum((cv1$var1.pred - cv1$observed)^2) / nrow(cv1))

## [1] 10.86099

Ahora, repetimos el proceso con los resultados correctos:

cv2 = gstat.cv(g2)

Tiempo de conversión:

cv2 = st_as_sf(cv2)

Calcule RSME para obtener resultados correctos:

# Este es el valor RMSE para la interpolación OK
sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))

## [1] 9.34399

Escribe un párrafo que resuma tus resultados. Según ellos, qué técnica de interpolación parece más precisa. Explicar.

Los resultados obtenidos al utilizar gstat.cv para evaluar las superficies climáticas generadas por los métodos de interpolación IDW y Kriging Ordinario revelan diferencias en la precisión de las predicciones. El cálculo del error cuadrático medio (RMSE) proporciona una medida objetiva de la discrepancia entre los valores predichos y observados. En este caso, se obtuvo un RMSE inicial de 10.86099, y después de una recalibración, el RMSE se redujo a 9.34399. Un RMSE menor indica una mayor precisión en las predicciones, ya que refleja una menor magnitud de errores entre las estimaciones y las observaciones reales. En consecuencia, basándonos en estos resultados, el método de interpolación que produce el RMSE más bajo, específicamente 9.34399 después de la recalibración, sugiere que el Kriging Ordinario es potencialmente más preciso en la predicción de las superficies climáticas en comparación con el método IDW. La elección del método de interpolación más preciso es crucial para garantizar resultados confiables en aplicaciones relacionadas con el clima y puede respaldar decisiones informadas en diversas disciplinas.

7. REFERENCIAS

Lizarazo, I., 2023. Spatial interpolation of soil organic carbon. obtenido de: rehttps://rpubs.com/ials2un/soc_interp.

sessionInfo()

## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## 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] curl_5.1.0    dplyr_1.1.3   ggplot2_3.4.4 leafem_0.2.3  leaflet_2.2.1
##  [6] automap_1.1-9 gstat_2.1-1   stars_0.6-4   abind_1.4-5   sf_1.0-14    
## [11] terra_1.7-55  sp_2.1-1     
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.4       xfun_0.40          bslib_0.5.1        raster_3.6-26     
##  [5] htmlwidgets_1.6.3  lattice_0.21-8     vctrs_0.6.3        tools_4.3.1       
##  [9] crosstalk_1.2.1    generics_0.1.3     parallel_4.3.1     tibble_3.2.1      
## [13] proxy_0.4-27       spacetime_1.3-0    fansi_1.0.4        highr_0.10        
## [17] xts_0.13.1         pkgconfig_2.0.3    KernSmooth_2.23-21 lifecycle_1.0.3   
## [21] farver_2.1.1       compiler_4.3.1     FNN_1.1.3.2        munsell_0.5.0     
## [25] codetools_0.2-19   htmltools_0.5.7    class_7.3-22       sass_0.4.7        
## [29] yaml_2.3.7         pillar_1.9.0       jquerylib_0.1.4    ellipsis_0.3.2    
## [33] classInt_0.4-10    cachem_1.0.8       tidyselect_1.2.0   digest_0.6.33     
## [37] fastmap_1.1.1      grid_4.3.1         colorspace_2.1-0   cli_3.6.1         
## [41] magrittr_2.0.3     base64enc_0.1-3    utf8_1.2.3         e1071_1.7-13      
## [45] withr_2.5.1        scales_1.2.1       rmarkdown_2.25     zoo_1.8-12        
## [49] png_0.1-8          evaluate_0.23      knitr_1.45         rlang_1.1.1       
## [53] Rcpp_1.0.11        glue_1.6.2         DBI_1.1.3          rstudioapi_0.15.0 
## [57] reshape_0.8.9      jsonlite_1.8.7     R6_2.5.1           plyr_1.8.9        
## [61] intervals_0.15.4   units_0.8-4