1. Introduccion

Este cuaderno es realizado para presentar un informe de la materia de Geomática Basica donde se ilustran dos técnicas de interpolación espacial: Distancia Inversa Ponderada (IDW) y Kriging Ordinario (OK). IDW es una técnica determinista. OK es una técnica probabilística. Ambas técnicas se utilizan aquí para obtener una superficie continua de carbono orgánico del suelo (SOC) a 15-30 cm a partir de muestras obtenidas de SoilGrids 250 m para el departamento de Cundinamarca.

2. Setup

Limpiamos la memoria y llamamos librerias que posteriormenete fueron instaladas en nuestra consola

rm(list=ls())
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)
library(knitr)
## 
## Adjuntando el paquete: 'knitr'
## The following object is masked from 'package:terra':
## 
##     spin

3. Leer los datos de entrada

list.files(path ="cundinamarca", pattern = "*.gpkg")
## [1] "cundinamarca Aguas_P1.gpkg"          
## [2] "cundinamarca Carreteras_P1.gpkg"     
## [3] "cundinamarca Ciudades_P1.gpkg"       
## [4] "cundinamarca Depto_Cundinamarca.gpkg"
## [5] "cundinamarca mun_Cun.gpkg"           
## [6] "cundinamarca Rios_P1.gpkg"           
## [7] "soc_cundinamarca.gpkg"               
## [8] "Suelos_Cundinamarca.gpkg"

Este paso es muy importante, pues debemos tener guardados los datos de soc del cuaderno anterior en la misma carpeta donde tenemos los demas datos y donde tenemos guardado el Spatial_interpolation.Rmd

samples <- sf::st_read("cundinamarca/soc_cundinamarca.gpkg")
## Reading layer `soc_cundinamarca' from data source 
##   `C:\Users\diazo\OneDrive\Documentos\GEOMATICA\PROYECTO_5\cundinamarca\soc_cundinamarca.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 2000 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.91012 ymin: 3.827538 xmax: -72.91554 ymax: 5.824439
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs

Aca abrimos los datos de SOC (contenido organico en el suelo) con nuestra función st_read.

Como nuestra zona de estudio es el departamento de Cundinamarca, ahora vamos a leer también el geopackage de municipios correspondiente:

munic <- sf::st_read("cundinamarca/cundinamarca mun_Cun.gpkg")
## Reading layer `mgn_adm_mpio_grafico' from data source 
##   `C:\Users\diazo\OneDrive\Documentos\GEOMATICA\PROYECTO_5\cundinamarca\cundinamarca mun_Cun.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 116 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.89063 ymin: 3.730129 xmax: -73.05256 ymax: 5.837258
## Geodetic CRS:  MAGNA-SIRGAS

4. Exploramos datos de entrada

summary(samples)
##       soc                   geom     
##  Min.   : 0.00   POINT        :2000  
##  1st Qu.:17.39   epsg:NA      :   0  
##  Median :22.46   +proj=long...:   0  
##  Mean   :21.13                       
##  3rd Qu.:25.79                       
##  Max.   :43.88

Aca observamos que los datos obtenidos fueron 2000

hist(samples$soc)

Este histograma nos sirve como una representación grafica de la distribución de los datos. En el eje X tenemos valores de SOC en intervalos y en el eje Y tenemos la frecuencia de de los valores. Vamos a dejar los valores de SOC en dos digitos

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

Ahora, vamos a visualizar las muestras para conocer su distribución espacial.

pal <- colorNumeric(
  c("#EF15CA","#EDE574","#FC913A","#FF4E50"),
  domain = samples$soc
)
leaflet() %>%
  addPolygons(
    data = munic,
    color = "black",
    opacity = 0.5,
    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")
## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'

En este mapa observamos la cantidad de carbono organico en el suelo por cada municipio del departamento de Cundinamarca, en este caso observamos valores de 0 en municipios que se encuentran cerca a la capital del pais, por otro lado observamos valores de casi 40 en municipios ubicados al norte del departamento y en municipios ubicados al occidente del departamento se obervan el mayor numero de datos intermedios (20 a 30).

5. Interpolacion

5.1 Creacion del Objeto gstat

La creación del objeto gstat es importante pues permite el análisis de datos y la interpolación mediantes tecnicas de kriking que son un conjunto de métodos geoestadísticos de interpolación que permiten predecir valores en ubicaciones no muestreadas calculando en la estructura espacial de los datos.

Basándose en estos, la función gstat entiende qué tipo de modelo de interpolación queremos utilizar: Sin modelo de variograma → Distancia de ponderación inversa (IDW) Modelo de variograma, sin covariables → Kriging ordinario (OK)

Ahora sabiendo esto vamos a visualizar las muestras para conocer su distribución espacial.

5.2 Interpolacion IDW

Para interpolar el SOC utilizando el método IDW primero se crea el siguiente objeto gstat, especificando sólo la fórmula y los datos:

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

Con nuestro modelo de interpolación g1 ya definido, podemos utilizar la función predecir para interpolar realmente, es decir, para estimar los valores SOc, esto para que el codigo sepa que los datos soc son los estamos usando.

(rrr <- terra::rast(xmin=-74.89063, xmax=-72.91554, ymin=3.730129, ymax=5.837258, nrows=370, ncols = 329, vals = 1, crs = "epsg:4326" ))
## class       : SpatRaster 
## dimensions  : 370, 329, 1  (nrow, ncol, nlyr)
## resolution  : 0.006003313, 0.005694943  (x, y)
## extent      : -74.89063, -72.91554, 3.730129, 5.837258  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : lyr.1 
## min value   :     1 
## max value   :     1

Aca es importante tomar nuestros propios datos de X y Y tanto mínimo como máximo.

Ahora, vamos a convertir este SpatRaster en un objeto de estrellas.

stars.rrr <-stars::st_as_stars(rrr)

La funcion “stars” nos sirve para trabajar con datos espaciales y raster en R convirtiendo datos rrr que pueden ser de elevación y asi poder visualizarlos.

La siguiente expresión interpola los valores SOC según el modelo definido en g1 y la trama definida en stars.rrr

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

Tenemos que ver 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  0.04106756 18.29731 22.44124 21.10835 24.62551 43.85827      0
## var1.var           NA       NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -74.89  0.006003 WGS 84 [x]
## y    1 370  5.837 -0.005695 WGS 84 [y]

Podemos subconjuntar sólo el primer atributo y cambiarle el nombre, y así podemos manejar los valores NA

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

Con este código lo que hacemos es que a los datos NA le vamos a dar valores de o.o

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                  Min.  1st Qu.   Median     Mean  3rd Qu.     Max.   NA's
## var1.pred  0.04106756 18.29731 22.44124 21.10835 24.62551 43.85827      0
## var1.var           NA       NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -74.89  0.006003 WGS 84 [x]
## y    1 370  5.837 -0.005695 WGS 84 [y]
names(z1) = "soc "

Y ya con la superficie interpolada, creamos una paleta de colores:

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

Ahora vemos como quedo la interpolación en un gráfico:

m <- leaflet() %>%
  addTiles() %>%
  leafem::addGeoRaster(
    z1,
    opacity = 0.4,
    colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"),
                                domain = 10:100)
  ) %>%
 addCircleMarkers(
   data = samples,
   radius = 1,
   label = ~soc,
   color = ~paleta(soc),
   fillOpacity = 1,
   stroke = F
 ) %>%
  addLegend("bottomright", pal = paleta, values = z1$soc,
  title = "Interpolacion IDW SOC [SOC]"
  )
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
## Warning in paleta(soc): Some values were outside the color scale and will be
## treated as NA
## Warning in paleta(soc): Some values were outside the color scale and will be
## treated as NA
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m

Este mapa lo que nos muestra son los distintos valores de SOC en Cundinamarca donde observamos valores que van desde 10 hasta los 35, predominando los valores cercanos a 20 y 25, esta función de interpolacion nos permite estar calculando ademas los valores medidos en puntos cercanos donde no se han tomado muestras.

5.3 Interpolacion OK

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 de los datos y asignar pesos en consecuencia a la hora de hacer predicciones.

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

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

Se traza el variograma empirico

plot(v_emp_ok)

Esta gráfica nos explica la estrutura espacial de los datos antes de aplicar metodos Kriging. En el eje X tenemos la distancia o separación entre pares de puntos y en el eje Y tenemos la semivarianza que nos muestra la variabilidad entre puntos separados por una distancia x.

Para ajustar un modelo de variograma a un variograma empírico utilizaremos 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 afina sus parámetros. Puede utilizar show.vgms() para mostrar varios tipos de modelos de variograma.

El modelo ajustado puede representarse con plot:

plot(v_mod_ok)

Esta gráfica representa la distancia entre puntos en el eje X y la semivarianza en el eje Y, en este caso que nuestra curva es creciente interpretamos que los puntos separados son valores mas diferentes a medida a que la distancia crece.

El siguiente condigo es clave para entender la estructura espacial de los datos y verificar que el modelo elegido representa bien la variabilidad del Carbono organico del suelo.

v_mod_ok$var_model
##   model     psill    range kappa
## 1   Nug  2.891776   0.0000   0.0
## 2   Ste 79.765628 157.9428   0.3

Aca el model muestra el tipo de modelo ajustado, psill representa la variabilidad total y range representa la distancia en la que los valores dejan de estar correlacionados espacialmente.Esto es importante pues al tener un rango corto podriamos decir que el SOC varia rapidamente en el espacio.

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

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

Una vez más, vamos a subconjuntar el atributo de valores ya dichos y a cambiarle el nombre

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

Las predicciones OK se muestran en el siguiente grafico.

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 = "Interpolacion OK SOC [%]"
    )
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
## Warning in paleta(soc): Some values were outside the color scale and will be
## treated as NA
## Warning in paleta(soc): Some values were outside the color scale and will be
## treated as NA
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m

Este grafico nos muestra la interpolacion para el carbono organico del suelo en Cundinamarca permitiendonos identificar areas con valores mas altos y mas bajos con mayor presición.

6. Evaluación de resultados

Para nuestra evaluacion de resultados vamos hacer una evaluacion cualitativa y una validacion cruzada.

6.1 Evaluacion cualitativa

colores <- colorOptions(palette = c("orange", "yellow", "cyan", "green"), domain = 10:100, na.color = "transparent")
m <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(z1, colorOptions = colores, opacity = 0.7, group= "IDW")  %>%
  addGeoRaster(z2, colorOptions = colores, opacity = 0.7, group= "OK")  %>%
  # Add layers controls
  addLayersControl(
    overlayGroups = c("IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
    addLegend("bottomright", pal = paleta, values= z1$soc,
    title = "Carbono organico en el suelo [%]"
)
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m

En este mapa podemos observar el raster de los datos de SOC con las dos capas (IDW y OK) donde se evidencian las areas con mayor porcentaje de, en el caso de IDW se observan colores mas oscuros y para OK tenemos colores mas claros o tenues considerando la evaluacion espacial de los datos. En general, al tener suelos con mayor SOC en altitudes mas bajas se habla de mas acumulacion de materia organica y cuando tenemos suelos con menor porcentaje de SOC podemos estar hablando de suelos mas degradados o erosionados

6.2 Validacion cruzada

Hemos estimado las superficies SOC utilizando dos métodos diferentes: IDW y Kriging ordinario. Aunque resulta útil examinar y comparar los resultados gráficamente, también necesitamos una forma objetiva de evaluar la precisión de la interpolación. De este modo, 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:

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 utilizando la función gstat.cv, que acepta un objeto gstat.

## Este codigo lo cargamos en la consola
#cv2 = gstat.cv(g2)

Let’s convert the cv2 object:

# Corremos el codigo en la consola
#cv2 = st_as_sf(cv2)

Let’s plot the residuals:

# run from the console
#sp::bubble(as(cv2[, "residual"], "Spatial"))
#sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))

Con el siguiente codigo revisamos nuetsros datos para asi mirar como esta guardado nuestro objeto residual

#print(cv2)

Gráfico de residuos (grafico de burbujas)

Ya sabiendo eso vamos a visualizar los residuos de un modelo geoestadístico en un mapa, utilizando burbujas de diferentes tamaños para representar la magnitud de los residuos.

include_graphics("C:/Users/diazo/OneDrive/Documentos/GEOMATICA/PROYECTO_5/Residual.jpg")

Segun el gráfico anterior tenemos que -28.273 es nuestro error negativo máximo y 15.189 es nuestro error positivo máximo que nos indica que el modelo tiene problemas en esas áreas donde tenemos las burbujas grandes quizás por datos insuficientes o un mal ajuste del semivariograma en Kriging, ademas como la distribucion de los colores de las burbujas es homogenea podemos hablar de errores aleatorios sin saber exactamente porque.

7. Guardar las salidas de interpolación

Finalmente guardamos nuestro trabajo en la carpeta donde tenemos los demas datos de nuestro proyecto 5

write_stars(
  z1, dsn = "cundinamarca/IDW_soc_cun.tif", layer = 1
)
write_stars(
  z2, dsn = "cundinamarca/OK_soc_cun.tif", layer = 1
)

8. Conclusiones

Los datos de carbono organico en Cundinamarca aunque son variables nos muestran suelos donde predominan buenos porcentajes de material organico y poco erosionados.

La interpolación IDW de SOC es útil para generar mapas de carbono orgánico del suelo a partir de datos de puntos de muestreo. Ayuda en la toma de decisiones en nuestro rol de Ingenieros agrónomos y en la gestión ambiental y estudios de cambio climático , aunque tiene limitaciones en cuanto a su precisión en comparación con métodos geoestadísticos más avanzados como kriging .

Para elegir el mejor metodo de interpolacion debemos tener en cuenta varias cosas pues IDW hace una estimacion mas rapida y sencilla pero OK tiene mas presicion y mayor análisis asi que a nivel general el metodo Ordinary Kriging (OK) es mejor.

9. Bibliografia

Lizarazo, I. 2025. Spatial interpolation. Disponible en: https://rpubs.com/ials2un/Spatial_Interpolation

Hahn, N. Creando mapas con R. Disponible en: https://bookdown.org/nicohahn/making_maps_with_r5/docs/leaflet.html.

Cheng, Karambelkar y Xie. 2015. Disponible en: https://rstudio.github.io/leaflet/

Diaz, O. 2025. Suelos Cuninamarca. Disponible en: https://rpubs.com/osdiaza/1274946

Diaz, A. 2002. Geoestadistica aplicada. Disponible en: https://www.esmg-mx.org/media/courses/geoestadistica/GeoEstadistica.pdf

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.48         leafem_0.2.3       leaflet_2.2.2      RColorBrewer_1.1-3
##  [5] automap_1.1-12     gstat_2.1-3        stars_0.6-8        abind_1.4-8       
##  [9] sp_2.1-4           sf_1.0-19          terra_1.8-15      
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.5            xfun_0.47               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-2         fansi_1.0.6             highr_0.11             
## [19] xts_0.14.0              pkgconfig_2.0.3         KernSmooth_2.23-24     
## [22] lifecycle_1.0.4         compiler_4.4.2          farver_2.1.2           
## [25] FNN_1.1.4.1             munsell_0.5.1           codetools_0.2-20       
## [28] htmltools_0.5.8.1       class_7.3-22            sass_0.4.9             
## [31] yaml_2.3.10             pillar_1.9.0            jquerylib_0.1.4        
## [34] classInt_0.4-10         cachem_1.1.0            tidyselect_1.2.1       
## [37] digest_0.6.37           dplyr_1.1.4             fastmap_1.2.0          
## [40] grid_4.4.2              colorspace_2.1-1        cli_3.6.3              
## [43] magrittr_2.0.3          base64enc_0.1-3         utf8_1.2.4             
## [46] e1071_1.7-16            scales_1.3.0            rmarkdown_2.28         
## [49] jpeg_0.1-10             zoo_1.8-12              png_0.1-8              
## [52] evaluate_1.0.0          rlang_1.1.4             Rcpp_1.0.13            
## [55] glue_1.8.0              DBI_1.2.3               rstudioapi_0.17.1      
## [58] reshape_0.8.9           jsonlite_1.8.9          R6_2.5.1               
## [61] plyr_1.8.9              intervals_0.15.5        units_0.8-5