Este cuadernos podemos observar y explorar dos tecnicas de interpolacion espacial ampliamente utilizadas: Inverse Distance Weighted (IDW) y Ordinary Kriging (OK). Mientras que IDW asigna valores a partir de una relación de proximidad ponderada, OK utiliza modelos estadísticos para estimar la incertidumbre en las predicciones. Ambas metodologías se aplican aquí para generar una superficie continua del carbono orgánico en el suelo (SOC) a una profundidad de 15-30 cm, empleando datos de SoilGrids 250 m.
Procedemos a instalar las bibliotecas que utilizaremos en este codigo.
# paquetes = c('terra', 'sf', 'sp', 'stars', 'gstat', 'automap', 'leaflet', 'leafem')
#install.packages(paquetes)
Despues de haber instalado procedemos a cargar las bibliotecas.
library(terra)
library(sf)
library(sp)
library(stars)
library(gstat)
library(automap)
library(RColorBrewer)
library(leaflet)
library(leafem)
library(knitr)
El siguiente código permite listar y cargar archivos en formato GeoPackage (.gpkg) dentro del entorno de R, un formato de almacenamiento de datos geoespaciales que es eficiente y ampliamente utilizado en SIG (Sistemas de Información Geográfica).
Se obtuvieron muestras de SOC a 15-30 cm de profundidad. dicha muestra la obtuvimos de SoilGrids 250m, lo cual procedemos a guardar las muestras en formato geopackage en el directorio de datos usando de nombre soc_Sucre.gpkg
list.files(path = "data", pattern = "*.gpkg")
## [1] "COL_water_Areas_sucre.gpkg"
## [2] "depto_sucre_corregido_roads.gpkg"
## [3] "depto_sucre_corregido_water_lines.gpkg"
## [4] "Mun_Sucre.gpkg"
## [5] "Municipio_Sucre.gpkg"
## [6] "Soc_15_30cm_Sucre.gpkg"
## [7] "Soc_15_30cm_Sucre.gpkg.aux.xml"
## [8] "soc_stder.gpkg"
## [9] "soc_Sucre.gpkg"
## [10] "Sucre_Mun.gpkg"
## [11] "Sucre_Mun.gpkg.tif"
## [12] "World_Cities_sucre.gpkg"
Este código se utiliza para cargar un archivo GeoPackage (.gpkg) en R, específicamente “soc_Sucre.gpkg”, que contiene datos espaciales relacionados con el carbono orgánico en el suelo (SOC) en la región de Sucre.
samples <- sf::st_read("data/soc_Sucre.gpkg")
## Reading layer `soc_Sucre' from data source
## `C:\Users\carlo\Desktop\INFORME 2\data\soc_Sucre.gpkg' using driver `GPKG'
## Simple feature collection with 1527 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -75.79167 ymin: 8.278027 xmax: -74.42446 ymax: 10.14538
## Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs
Este codigo carga un archivo GeoPackage (.gpkg) llamado “Mun_Sucre.gpkg”, que contiene información espacial sobre los municipios del departamento de Sucre.
munic <- sf::st_read("data/Mun_Sucre.gpkg")
## Reading layer `mapa' from data source
## `C:\Users\carlo\Desktop\INFORME 2\data\Mun_Sucre.gpkg' using driver `GPKG'
## Simple feature collection with 26 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -75.70603 ymin: 8.277882 xmax: -74.53271 ymax: 10.14548
## Geodetic CRS: MAGNA-SIRGAS
El siguiente código genera un resumen estadístico de los datos contenidos en la variable samples, que previamente cargó la capa espacial de muestras de carbono orgánico en el suelo (SOC).
Interpretacion de los resultados
SOC (carbono orgánico en el suelo): - Mínimo de 10.39 y máximo de 123.19, indicando el rango de valores de carbono en el suelo. La media (24.90) y la mediana (20.40) ayudan a entender la distribución de los datos.
Geometry: - Muestra que los datos están en formato espacial con geometría de tipo POINT, es decir, cada muestra es un punto georreferenciado.
summary(samples)
## soc geom
## Min. : 10.39 POINT :1527
## 1st Qu.: 15.82 epsg:NA : 0
## Median : 20.40 +proj=long...: 0
## Mean : 24.90
## 3rd Qu.: 30.12
## Max. :123.19
El histograma revela una distribución sesgada a la derecha, donde la mayoría de las muestras tienen niveles bajos a moderados de SOC, mientras que existen algunos valores altos menos frecuentes. Este análisis es clave para tomar decisiones informadas sobre el procesamiento de los datos y la aplicación de métodos de interpolación en estudios geoespaciales.
samples$soc <- as.numeric(as.character(samples$soc))
hist(samples$soc)
El siguiente código redondea los valores de carbono orgánico en el suelo (SOC) a dos decimales esto ayuda a mejorar la presentacion de datos evitando numeros con demaciados decimales facilitando la interpretacion en tablas y graficos
samples$soc = round(samples$soc,2)
Este código permite visualizar la distribución espacial del SOC en una región específica, facilitando el análisis de patrones geoespaciales.
Es útil para: ✅ Identificar áreas con altos y bajos niveles de SOC. ✅ Comparar la distribución del SOC con los límites municipales. ✅ Facilitar la interpretación mediante un mapa interactivo.
pal <- colorNumeric(
c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),
domain = samples$soc,
)
leaflet() %>%
addPolygons(
data = munic,
color = "gray",
opacity = 1,
weight = 1,
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.9) %>%
addProviderTiles("OpenStreetMap")
Para realizar una interpolación espacial, primero necesitamos crear un objeto de clase gstat utilizando la función gstat(). Este objeto contiene toda la información necesaria para llevar a cabo el proceso de interpolación, incluyendo:
Tipos de modelos de interpolacion esta función gstat() determina el método de interpolación en función de los argumentos proporcionados:
Sin modelo de variograma → Inverse Distance Weighting (IDW) Enfoque determinista que asigna más peso a los puntos más cercanos.
Con modelo de variograma y sin covariables → Ordinary Kriging (OK) Enfoque probabilístico basado en la estructura espacial de los datos.
Al configurar la función gstat(), utilizamos tres parámetros principales:
En este fragmento de código, se crea un objeto gstat para realizar Interpolación Ponderada por la Inversa de la Distancia (IDW) sobre los valores de carbono orgánico en el suelo (SOC).
este código prepara la base para interpolar el SOC en un área de estudio, generando un modelo espacial continuo basado en los puntos de muestreo disponibles. 🌍📊
g1 = gstat(formula = soc ~ 1, data = samples)
En este paso, se crea un objeto de tipo raster utilizando la biblioteca terra en R. Este raster define el área geográfica donde se realizarán las predicciones de contenido de carbono orgánico en el suelo (SOC). Se establece una malla de celdas con valores iniciales de 1, abarcando la extensión espacial de interés.
¿Para qué sirve?
Define la ubicación de las predicciones: Especifica las coordenadas espaciales en las que se calcularán los valores interpolados de SOC.
Estandariza la estructura espacial: Permite que la interpolación se realice sobre una grilla homogénea con un sistema de referencia geográfica definido (crs = “epsg:4326”).
Facilita el uso de modelos geoespaciales: El raster generado puede ser utilizado en combinación con el modelo de interpolación (gstat) para predecir valores en ubicaciones sin datos de muestreo.
(rrr <- terra::rast(xmin=-75.70603, xmax=-74.53271, ymin=8.277882, ymax=10.14548, nrows=370, ncols = 329, vals = 1, crs = "epsg:4326"))
## class : SpatRaster
## dimensions : 370, 329, 1 (nrow, ncol, nlyr)
## resolution : 0.003566322, 0.005047562 (x, y)
## extent : -75.70603, -74.53271, 8.277882, 10.14548 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : lyr.1
## min value : 1
## max value : 1
En este paso, se convierte un objeto de tipo SpatRaster (rrr) en un objeto stars utilizando la función st_as_stars() de la librería stars. Esta transformación es necesaria porque algunos métodos de interpolación espacial y visualización en R trabajan mejor con objetos stars en lugar de SpatRaster.
¿Para qué sirve?
Optimiza el procesamiento espacial: Los objetos stars son más eficientes para ciertas operaciones de análisis espacial.
Facilita la compatibilidad con otros paquetes: Algunos paquetes en R requieren datos en formato stars para la interpolación y la manipulación geoespacial.
Reduce el uso de memoria: Para sistemas con recursos limitados, la conversión a stars puede ser útil para mejorar el rendimiento en cálculos espaciales.
stars.rrr <- stars::st_as_stars(rrr)
En este proceso, se emplea el modelo de interpolación definido en g1 y el raster stars.rrr para estimar los valores de SOC (Carbono Orgánico del Suelo) en ubicaciones no muestreadas. La interpolación se realiza mediante el método Inverse Distance Weighted (IDW), que asigna valores a los puntos en función de la distancia a los datos observados.
¿Para qué sirve?
Estimación de valores en áreas no muestreadas: Permite obtener datos interpolados en ubicaciones donde no se realizaron mediciones directas.
Análisis espacial de suelos: Útil en estudios ambientales, agrícolas y de manejo de recursos naturales.
Visualización y generación de mapas: Se puede usar el resultado para crear mapas de distribución de SOC y analizar patrones espaciales.
z1 = predict(g1, stars.rrr)
## [inverse distance weighted interpolation]
z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## var1.pred 10.78552 19.24331 23.27773 25.68444 31.15314 121.4536 0
## var1.var NA NA NA NaN NA NA 121730
## dimension(s):
## from to offset delta refsys x/y
## x 1 329 -75.71 0.003566 WGS 84 [x]
## y 1 370 10.15 -0.005048 WGS 84 [y]
El objeto z1 generado en la interpolación contiene dos atributos principales:
Podemos extraer únicamente el primer atributo (var1.pred), que representa los valores interpolados de SOC, y asignarle un nuevo nombre para facilitar su uso en análisis posteriores. Esto permite una manipulación más clara y específica de los datos interpolados, optimizando su visualización y aplicación.
a <- which(is.na(z1[[1]]))
z1[[1]][a] = 0.0
z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## var1.pred 10.78552 19.24331 23.27773 25.68444 31.15314 121.4536 0
## var1.var NA NA NA NaN NA NA 121730
## dimension(s):
## from to offset delta refsys x/y
## x 1 329 -75.71 0.003566 WGS 84 [x]
## y 1 370 10.15 -0.005048 WGS 84 [y]
names(z1) = "soc"
El código proporciona una visualización de la interpolación espacial de SOC (Carbono Orgánico del Suelo) utilizando la biblioteca leaflet en R. Para ello, se define una paleta de colores que representa los valores interpolados y se genera un mapa interactivo con los datos procesados.
¿Para qué sirve?
Este script sirve para representar gráficamente los resultados de la interpolación espacial mediante un mapa interactivo. Utiliza la función addGeoRaster para mostrar la superficie interpolada y addCircleMarkers para superponer los puntos de muestreo. Además, se agrega una leyenda (addLegend) que permite interpretar los valores en el mapa, facilitando el análisis visual de la distribución espacial del SOC.
paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), domain = 10:120, na.color = "transparent")
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
El variograma es una herramienta fundamental en los métodos de interpolación Kriging, ya que permite cuantificar la autocorrelación espacial de los datos. En este caso, el código utiliza la función variogram de la librería gstat en R para calcular el variograma empírico de la variable soc contenida en el conjunto de datos samples. Posteriormente, se genera su visualización con la función plot.
Este procedimiento es útil para entender la estructura espacial de los datos y determinar cómo varía la similitud entre observaciones a diferentes distancias. El variograma es un paso clave en la construcción de modelos Kriging, ya que permite definir un modelo de variograma teórico que se utilizará en la interpolación espacial de datos.
v_emp_ok = variogram(soc ~ 1, data=samples)
plot(v_emp_ok)
El código utiliza la función autofitVariogram del paquete automap en R para ajustar automáticamente un modelo de variograma a un variograma empírico. Esto simplifica el proceso de selección de un modelo adecuado, evitando la necesidad de ajustar manualmente los parámetros del variograma.
Este proceso es fundamental en la interpolación Kriging, ya que el modelo de variograma describe la dependencia espacial de los datos. Ajustar un modelo adecuado permite realizar predicciones espaciales más precisas y mejorar la calidad de los resultados en análisis geoespaciales.
v_mod_ok = automap::autofitVariogram(soc ~ 1, as(samples, "Spatial"))
La función autofitVariogram selecciona automáticamente el tipo de modelo de variograma que mejor se ajusta a los datos y ajusta sus parámetros óptimos. Para visualizar los distintos tipos de modelos de variograma disponibles, se puede usar la función show.vgms().
Dado que autofitVariogram no es compatible con objetos sf, es necesario convertir el objeto samples en un SpatialPointsDataFrame utilizando el paquete sp.
El modelo ajustado se puede visualizar mediante la función plot().
plot(v_mod_ok)
La tabla presentada muestra los parámetros del modelo de variograma ajustado, extraídos de la componente $var_model del objeto v_mod_ok. Este modelo incluye distintos parámetros clave que describen la estructura espacial de la variable analizada (soc), tales como:
El modelo de variograma ajustado es esencial para la interpolación Kriging, ya que permite caracterizar la variabilidad espacial de los datos y predecir valores en ubicaciones no muestreadas. Con esta información, se pueden generar mapas de interpolación más precisos, optimizando la estimación de variables espaciales.
v_mod_ok$var_model
## model psill range kappa
## 1 Nug 5.15323 0.000 0.0
## 2 Ste 875.94471 3884.466 0.2
En este codigo muestra la implementación de la interpolación espacial mediante Kriging Ordinario utilizando la función gstat en R. En primer lugar, se construye un modelo geoestadístico (g2) basado en el variograma previamente ajustado (v_mod_ok$var_model) y en los datos de muestreo (samples). Posteriormente, se usa la función predict para generar una superficie interpolada, aplicando el modelo de Kriging a un objeto ráster (stars.rrr).
Este procedimiento permite predecir valores de una variable en ubicaciones no muestreadas a partir de datos espaciales existentes, aprovechando la autocorrelación espacial. Es una técnica ampliamente utilizada en análisis geoespacial.
g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2= predict(g2, stars.rrr)
## [using ordinary kriging]
Nuevamente, crearemos un subconjunto del 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 correctas 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
El código emplea la librería leaflet en R para visualizar de manera interactiva los resultados de dos técnicas de interpolación espacial: IDW (Inverse Distance Weighting) y OK (Ordinary Kriging). Utiliza addGeoRaster para agregar capas ráster con diferentes métodos de interpolación y addLayersControl para permitir la alternancia entre ellas. Además, se implementa una paleta de colores personalizada y una leyenda explicativa para facilitar la interpretación de los datos.
Esta representación visual interactiva es útil para comparar la calidad y diferencias entre los métodos de interpolación utilizados. Ayuda a identificar cuál de las técnicas produce mejores estimaciones de la variable analizada (en este caso, carbono orgánico en el suelo)
colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), domain = 10:100, na.color = "transparent")
m <- leaflet() %>%
addTiles() %>%
addGeoRaster(z1, colorOptions = colores, opacity = 0.8, group= "IDW") %>%
addGeoRaster(z2, colorOptions = colores, opacity = 0.8, 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
Este apartado presenta la validación cruzada Leave-One-Out (LOOCV) como un método para evaluar la precisión de dos técnicas de interpolación: IDW (Inverse Distance Weighting) y Kriging Ordinario. Este procedimiento consiste en eliminar temporalmente un punto de la base de datos, predecir su valor y comparar con el dato real, repitiendo el proceso para todos los puntos.
El análisis se realiza mediante la función gstat.cv en R, que devuelve una tabla con valores observados, predichos y otros indicadores estadísticos como la varianza y el Z-score en el caso del Kriging.
Este método permite comparar de manera objetiva la exactitud de distintos modelos de interpolación y seleccionar el más adecuado para la estimación de variables espaciales, como el carbono orgánico del suelo (SOC).
# este código debe ejecutarse desde la consola, no desde aquí
# cv2 = gstat.cv(g2)
En este apartado se presentan los atributos generados por la función gstat.cv en R, utilizada para evaluar la precisión de interpolaciones espaciales. Los resultados incluyen valores predichos, observados, residuales, varianza (para Kriging) y otros indicadores de error.
Para facilitar su análisis y representación geoespacial, se convierte el objeto cv2 al formato espacial sf mediante la función st_as_sf(cv2).
# Ejecutar desde la consola
# cv2 = st_as_sf(cv2)
Trazamos los residuos
# Ejecute desde la consola
# sp::bubble(as(cv2[, "residual"], "Spatial"))
include_graphics("C:/Users/carlo/Desktop/INFORME 2/data/Residual.png.jpg")
# Este es el valor RMSE para la interpolación OK
# ejecutar desde la consola
# sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))
En este apartado, se almacenan los resultados de las interpolaciones realizadas mediante los métodos IDW (Inverse Distance Weighting) y Kriging Ordinario (OK). Se utilizan archivos en formato TIFF para guardar las capas interpoladas (IDW_soc_Sucre.tif y OK_soc_Sucre.tif)
Guardar las interpolaciones en formato raster permite su posterior análisis en software SIG como QGIS o ArcGIS. Estos archivos pueden ser utilizados para la evaluación de la calidad del modelo, visualización en mapas y comparación de distintos métodos de interpolación.
write_stars(
z1, dsn = "data/IDW_soc_Sucre.tif", layer = 1
)
write_stars(
z2, dsn = "data/OK_soc_Sucre.tif", layer = 1
)
Al escribir este cuaderno aprendimos y aplicamos dos técnicas fundamentales de interpolación espacial: Inverse Distance Weighted (IDW) y Ordinary Kriging (OK). Aprendimos que IDW es un método determinista que asigna mayor peso a los puntos más cercanos, mientras que OK es un enfoque geoestadístico que modela la estructura espacial de los datos a través de un variograma, permitiendo estimaciones más precisas y con incertidumbre cuantificada. Mediante la manipulación de datos geoespaciales, el análisis estadístico y la visualización en mapas interactivos, comprendimos la importancia de seleccionar adecuadamente la metodología de interpolación según la naturaleza y distribución de los datos. Este proceso no solo fortaleció nuestras habilidades en R sino que también nos permitió obtener una representación espacial detallada. Además, comprendí la importancia del preprocesamiento de datos geoespaciales, la visualización interactiva con leaflet y el uso de herramientas como gstat y automap para la modelización geoestadística.
Principalmente se uso Lizarazo, I., 2025.Spatial interpolation. Available at: https://rpubs.com/ials2un/Spatial_Interpolation
sessionInfo()
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
##
## 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] knitr_1.49 leafem_0.2.3 leaflet_2.2.2 RColorBrewer_1.1-3
## [5] automap_1.1-16 gstat_2.1-3 stars_0.6-8 abind_1.4-8
## [9] sp_2.1-4 sf_1.0-19 terra_1.8-10
##
## 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-30 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 jpeg_0.1-10
## [49] zoo_1.8-12 png_0.1-8 evaluate_1.0.1
## [52] rlang_1.1.4 Rcpp_1.0.13-1 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