Interpolación espacial del carbono organico del suelo
Gabriel Andres Gutierrez Gomez
20/11/2023
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.2library(terra)## terra 1.7.46library(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 TRUElibrary(stars)## Loading required package: abindlibrary(gstat)## Warning: package 'gstat' was built under R version 4.3.2library(automap)## Warning: package 'automap' was built under R version 4.3.2library(leaflet)## Warning: package 'leaflet' was built under R version 4.3.2library(leafem)
library(ggplot2)
library(dplyr)## Warning: package 'dplyr' was built under R version 4.3.2##
## 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, unionlibrary(curl)## Using libcurl 7.84.0 with Schannelh <- 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 : 1069, 1173, 1 (nrow, ncol, nlyr)
## resolution : 250, 250 (x, y)
## extent : -8326000, -8032750, 518250, 785500 (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine
## source : soc_igh_15_30.tif
## name : soc_igh_15_30Ahora, se convierten los datos del SOC en porcentaje.
soc.perc <- soc/10¿Cuál es el CRS de los datos del mundo real? los datos de SoilGrids suelen utilizar el sistema de referencia espacial (CRS, por sus siglas en inglés) WGS 84 (World Geodetic System 1984). WGS 84 es un sistema de coordenadas geográficas ampliamente utilizado en sistemas de información geográfica (SIG) y sistemas de posicionamiento global (GPS).
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 : 1091, 1256, 1 (nrow, ncol, nlyr)
## resolution : 0.002199529, 0.002199529 (x, y)
## extent : -74.71029, -71.94768, 4.65658, 7.056267 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## source(s) : memory
## name : soc_igh_15_30
## min value : 8.874964
## max value : 293.181030Convirtamos 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 map4. 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 : -74.70919, -71.94878, 4.65768, 7.050768 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names : soc_igh_15_30
## type : <num>
## values : 56.71
## 39.4
## 47.88Describa las características principales del objeto de muestra.
Es un objeto de tipo SpatVector que contiene datos espaciales. Sus características detalladas se presentan a continuación: Clase (class): SpatVector. Esto indica que el objeto es un vector espacial, que probablemente contenga información geoespacial.
Geometría (geometry): Los datos están representados como puntos. En el contexto geoespacial, esto significa que los datos están asociados con ubicaciones específicas en la superficie de la Tierra.
Dimensiones (dimensions): El vector tiene dimensiones de 500 geometrías y 1 atributo. Esto significa que hay 500 puntos espaciales, y cada punto tiene asociado al menos un atributo.
Extensión (extent): La extensión geográfica de los datos está definida por los valores: -74.70919 (xmin), -71.94878 (xmax), 4.65768 (ymin), 7.050768 (ymax). Estos valores representan las coordenadas geográficas mínimas y máximas en longitud y latitud.
- Coordenadas de referencia (coord. ref.): Las coordenadas de referencia indican el sistema de referencia espacial utilizado. En este caso, se utiliza el sistema de coordenadas geográficas WGS84 (+proj=longlat +datum=WGS84 +no_defs).
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_30id <- 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 mapAhora 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, a 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)Que es rrr?
rrr## class : SpatRaster
## dimensions : 273, 314, 1 (nrow, ncol, nlyr)
## resolution : 0.008798118, 0.008798118 (x, y)
## extent : -74.71029, -71.94768, 4.65438, 7.056267 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## source(s) : memory
## name : soc_igh_15_30
## min value : 9.67114
## max value : 147.83141Definir nuevos valores:
values(rrr) <-1Definir nuevos nombres:
names(rrr) <- "valor"Que es rrr ahora?
rrr## class : SpatRaster
## dimensions : 273, 314, 1 (nrow, ncol, nlyr)
## resolution : 0.008798118, 0.008798118 (x, y)
## extent : -74.71029, -71.94768, 4.65438, 7.056267 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## source(s) : memory
## name : valor
## min value : 1
## max value : 1stars.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]Que 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 11.92467 33.07971 38.62312 39.06538 45.96558 89.63706 0
## var1.var NA NA NA NaN NA NA 85722
## dimension(s):
## from to offset delta refsys x/y
## x 1 314 -74.71 0.008798 +proj=longlat +datum=WGS8... [x]
## y 1 273 7.056 -0.008798 +proj=longlat +datum=WGS8... [y]Tome nota de los nombres de los dos atributos incluidos dentro del objeto z1.
Columna x: Esta columna parece representar la variable x asociada con cada geometría. En este caso, los valores son 1, 314, y -74.71 para las tres geometrías.
Columna y: Similar a la columna x, esta columna parece representar la variable y asociada con cada geometría. Los valores son 1, 273, y 7.056 para las tres geometrías.
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 map5.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_modelEn tu cuaderno, explica el significado de cada elemento en el modelo anterior.
-Columna model: Tipo de modelo (
-Columna psill: Nivel de sill (
Columna range: Rango (
): El rango del variograma indica la distancia a la cual la variabilidad espacial alcanza aproximadamente el 95% del sill. Para el modelo “Nug”, el valor es 0.00000 (indicando que la variabilidad espacial se alcanza instantáneamente), y para el modelo “Ste”, el valor es 64.23925. Columna kappa: Parámetro de suavidad (
): Este parámetro controla la forma del variograma esférico y afecta la rapidez con la que la correlación disminuye a medida que aumenta la distancia. Para el modelo “Nug”, el valor es 0.0 (indicando que no hay suavidad adicional), y para el modelo “Ste”, el valor es 0.3.
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 mapa6. 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 mapa6.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.
- 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 usando la función gstat.cv, que acepta un objeto gstat.
Al escribir el siguiente fragmento, oculte el mensaje y los resultados.
### correr el siguiente codigo 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 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 : 478
## extent : -74.685, -72.01257, 4.65768, 7.050768 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 6
## names : var1.pred, var1.var, observed, residual, zscore, fold
## min values : 15.7896278955756, NA, 11.7452383041382, -23.0864999703684, NA, 1
## max values : 62.7843156942805, NA, 90.4286422729492, 46.3445951971944, NA, 478Ahora 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):
#Este es el valor RMSE para la interpolación IDW
sqrt(sum((cv1$var1.pred - cv1$observed)^2) / nrow(cv1))## [1] 10.65261Ahora, repetimos el proceso con los resultados correctos:
Tiempo de conversión:
cv2 = gstat.cv(g2)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.618802Escribe un párrafo que resuma tus resultados. Según ellos, qué técnica de interpolación parece más precisa. Explicar.
Los resultados de las técnicas de interpolación indican que la técnica de interpolación de Kriging Ordinario (OK) muestra una mayor precisión en comparación con la técnica de Inverso de la Distancia Ponderado (IDW), según el valor del Error Cuadrático Medio (RMSE). El RMSE para la interpolación OK es de 9.618802, mientras que para la interpolación IDW es de 10.65261. El RMSE es una medida que cuantifica la diferencia entre los valores observados y los valores predichos, y en este caso, valores más bajos indican una mayor precisión. Por lo tanto, el RMSE más bajo asociado con la interpolación OK sugiere que esta técnica proporciona predicciones más precisas en comparación con IDW. La elección entre técnicas de interpolación depende mucho de la naturaleza específica de los datos y los patrones espaciales, para este caso, los resultados por Kriging Ordinario muestran un método más preciso para la interpolación.
7. Referencias
Lizarazo, I., 2023. Spatial interpolation of soil organic carbon. obtenido de: https://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 11 x64 (build 22621)
##
## 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.0.2 dplyr_1.1.3 ggplot2_3.4.3 leafem_0.2.3 leaflet_2.2.0
## [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-46 sp_2.1-1
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.3 xfun_0.40 bslib_0.5.1 raster_3.6-26
## [5] htmlwidgets_1.6.2 lattice_0.21-8 vctrs_0.6.3 tools_4.3.1
## [9] crosstalk_1.2.0 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.4
## [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.6 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.2 scales_1.2.1 rmarkdown_2.24 zoo_1.8-12
## [49] png_0.1-8 evaluate_0.21 knitr_1.43 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