1. Introducción.

En esta práctica se compararon dos métodos de interpolación espacial para estimar el contenido de Carbono Orgánico del Suelo (SOC). El primero fue la Distancia Inversa Ponderada (IDW), que realiza las predicciones dando mayor peso a los puntos más cercanos. El segundo fue el Kriging Ordinario (OK), un método geoestadístico que considera la relación espacial entre los datos mediante un variograma. Ambos métodos se aplicaron a datos de muestreo del departamento del Magdalena, obtenidos de la base de datos SoilGrids para la profundidad de 15–30 cm, con el propósito de generar una superficie continua de la distribución del SOC.

2. Configuración del entorno.

Limpiamos el espacio de trabajo para trabajar sin residuos de variables previas.

rm(list=ls())

Cargamos las librerías necesarias.

library(terra)
## Warning: package 'terra' was built under R version 4.5.3
## terra 1.9.34
library(sf)
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
library(sp)
## Warning: package 'sp' was built under R version 4.5.3
library(stars)
## Cargando paquete requerido: abind
library(gstat)
## Warning: package 'gstat' was built under R version 4.5.3
library(automap)
## Warning: package 'automap' was built under R version 4.5.3
library(RColorBrewer)
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.5.3
library(leafem)
## Warning: package 'leafem' was built under R version 4.5.3

#3. Lectura de los datos de entrada. Cargamos los puntos de muestreo de SOC (soc_magdalena.gpkg) y los límites municipales del Magdalena (Municipios_Magdalena.gpkg).

list.files(path = "C:/Users/suare/OneDrive/Documentos/GB2/Proyecto1/Datos", pattern = "\\.gpkg$")
##  [1] "areas_Magdalena.gpkg"        "buffer_magdalena.gpkg"      
##  [3] "cities_magdalena.gpkg"       "Curvas_Magdalena.gpkg"      
##  [5] "Departamento.magdalena.gpkg" "Departamentos.gpkg"         
##  [7] "magdalena_munic.gpkg"        "muestras_elevacion.gpkg"    
##  [9] "Municipios_9377.gpkg"        "Municipios_Magdalena.gpkg"  
## [11] "puntos_al.gpkg"              "puntos_aleatorios.gpkg"     
## [13] "Roads_Magdalena.gpkg"        "soc_magdalena.gpkg"         
## [15] "water_magdalena.gpkg"
samples <- sf::st_read("C://Users//suare//OneDrive//Documentos//GB2//Proyecto1//Datos//soc_magdalena.gpkg")
## Reading layer `soc_magdalena' from data source 
##   `C:\Users\suare\OneDrive\Documentos\GB2\Proyecto1\Datos\soc_magdalena.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1850 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.94689 ymin: 8.911589 xmax: -73.5435 ymax: 9.998805
## Geodetic CRS:  WGS 84
munic <- sf::st_read("C://Users//suare//OneDrive//Documentos//GB2//Proyecto1//Datos//Municipios_Magdalena.gpkg")
## Multiple layers are present in data source C:\Users\suare\OneDrive\Documentos\GB2\Proyecto1\Datos\Municipios_Magdalena.gpkg, reading layer `municipios_'.
## Use `st_layers' to list all layer names and their type in a data source.
## Set the `layer' argument in `st_read' to read a particular layer.
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## automatically selected the first layer in a data source containing more than
## one.
## Reading layer `municipios_' from data source 
##   `C:\Users\suare\OneDrive\Documentos\GB2\Proyecto1\Datos\Municipios_Magdalena.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 20 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.9569 ymin: 8.956799 xmax: -73.6021 ymax: 11.34958
## Geodetic CRS:  WGS 84

4. Exploración de los datos.

En este paso se realizó una exploración inicial de los datos de Carbono Orgánico del Suelo (SOC). Se revisó cómo se distribuyen los valores y su ubicación en el departamento del Magdalena, con el fin de conocer mejor la información antes de aplicar los métodos de interpolación. ## 4.1 Estadísticos descriptivos Primero, obtenemos el resumen estadístico de la variable soc:

summary(samples)
##       soc                   geom     
##  Min.   : 0.00   POINT        :1850  
##  1st Qu.:13.00   epsg:4326    :   0  
##  Median :15.60   +proj=long...:   0  
##  Mean   :17.73                       
##  3rd Qu.:20.30                       
##  Max.   :83.80

4.2 Histograma de frecuencias

Visualizamos la distribución de los datos para verificar si la variable tiene una distribución simétrica o si presenta sesgo (asimétrica hacia la derecha o izquierda):

hist(samples$soc,
     main = "Distribución del SOC en el Magdalena (15-30 cm)",
     xlab = "Carbono Orgánico del Suelo (SOC)",
     ylab = "Frecuencia",
     col = "darkseagreen",
     border = "white")

Análisis: El histograma muestra que la distribución del Carbono Orgánico del Suelo (SOC) no es normal, sino que presenta un sesgo positivo (hacia la derecha). La mayor parte de las muestras se concentra entre 10 y 30 unidades de SOC, mientras que los valores más altos son poco frecuentes. Esto indica que predominan concentraciones bajas a moderadas de carbono orgánico y que solo un pequeño número de muestras presenta valores elevados, generando la cola hacia la derecha de la distribución.

4.3 Redondeo de los valores de SOC

Para facilitar la lectura y evitar trabajar con decimales excesivos en las etiquetas del mapa, redondeamos la variable soc a dos dígitos:

samples$soc <- round(samples$soc, 2)

4.4 Visualización espacial

Para conocer cómo están distribuidas geográficamente las muestras sobre el territorio del Magdalena, creamos un mapa interactivo con la librería leaflet.

Primero, creamos una paleta de colores continua que asigne tonalidades según la cantidad de SOC:

pal <- colorNumeric(
  palette = c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),
  domain = samples$soc
)

Ahora, construimos el mapa superponiendo los límites municipales (munic) y los puntos de muestreo (samples):

leaflet() %>%
  addPolygons(
    data = munic,
    color = "gray",
    weight = 1,
    fillOpacity = 0.2
  ) %>%
  addCircleMarkers(
    data = samples,
    radius = 2,
    label = ~as.character(soc),
    color = ~pal(soc),
    fillOpacity = 0.9,
    stroke = FALSE
  ) %>%
  addLegend(
    data = samples,
    pal = pal,
    values = ~soc,
    position = "bottomleft",
    title = "SOC (15-30 cm):",
    opacity = 0.9
  ) %>%
  addProviderTiles("OpenStreetMap")

Análisis: El mapa muestra la distribución de los puntos de muestreo del Carbono Orgánico del Suelo (SOC) para la profundidad de 15–30 cm. Se observa que los puntos se concentran en la parte sur del área de estudio y presentan variaciones en los valores de SOC, representadas por la escala de colores. Debido a esta distribución, el muestreo no cubre completamente el departamento del Magdalena, por lo que algunas zonas quedan sin información directa. Durante el desarrollo de la práctica se intentó por diferentes medios obtener una distribución de puntos que abarcara todo el departamento; sin embargo, no fue posible debido a las limitaciones de los datos descargados desde SoilGrids. A pesar de ello, los datos obtenidos permitieron continuar con el análisis y aplicar los métodos de interpolación sobre el área disponible.

#5. Interpolación Espacial

5.1 Creación del objeto gstat

Para realizar una interpolación en R con la librería gstat, primero creamos un objeto de clase gstat. Este objeto define la fórmula del modelo (\(soc \sim 1\), lo que indica que no hay covariables externas) y vincula nuestros datos de calibración (samples).

g1 <- gstat(formula = soc ~ 1, data = samples)
ext_samples <- sf::st_bbox(samples)

Una vez definido el modelo de interpolación (g1), se utilizó la función predict para estimar los valores de SOC en las zonas sin datos. Para ello, se creó un raster con la librería terra, el cual define el área donde se realizarán las predicciones.

rrr <- terra::rast(
  xmin = ext_samples$xmin, 
  xmax = ext_samples$xmax, 
  ymin = ext_samples$ymin, 
  ymax = ext_samples$ymax, 
  nrows = 370, 
  ncols = 329, 
  vals = 1, 
  crs = "EPSG:4326"
)

Convertimos el SpatRaster a un objeto stars:

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

5.2 Interpolación por Distancia Inversa Ponderada (IDW)

El método IDW es un enfoque determinístico que asume que las observaciones más cercanas a un punto no muestreado tienen mayor influencia que aquellas situadas a mayor distancia.

Ejecutamos la función predict para estimar la superficie de SOC a partir del modelo g1 sobre la rejilla stars.rrr:

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

z1 representa la superficie continua interpolada de Carbono Orgánico del Suelo (SOC) generada mediante el método de Distancia Inversa Ponderada (IDW). Es un objeto de tipo stars que almacena la estimación determinística del valor de SOC para cada píxel dentro de la rejilla o grilla de predicción definida sobre el departamento del Magdalena.

z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                   Min.  1st Qu.   Median     Mean  3rd Qu.     Max.   NA's
## var1.pred  0.001938973 15.25994 16.72501 17.78451 19.42574 83.71578      0
## var1.var            NA       NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -74.95  0.004266 WGS 84 [x]
## y    1 370  9.999 -0.002938 WGS 84 [y]

Limpieza de valores NA

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  0.001938973 15.25994 16.72501 17.78451 19.42574 83.71578      0
## var1.var            NA       NA       NA      NaN       NA       NA 121730
## dimension(s):
##   from  to offset     delta refsys x/y
## x    1 329 -74.95  0.004266 WGS 84 [x]
## y    1 370  9.999 -0.002938 WGS 84 [y]

Asignamos el nombre de la variable al atributo

names(z1) <- "soc"

Visualización del resultado IDW

Construimos la paleta de colores y el mapa:

paleta <- colorNumeric(
  palette = c("orange", "yellow", "cyan", "green"), 
  domain = range(samples$soc, na.rm = TRUE), 
  na.color = "transparent"
)

Ahora, visualicemos el resultado de la interpolación:

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 [%]"
    )
## Warning in write_stars.stars(x, dsn = fl): all but first attribute are ignored
m

5.3 Interpolación por Kriging Ordinario (OK)

El proceso de Kriging Ordinario se divide en tres pasos principales:

Cálculo del variograma empírico:

Medimos la variabilidad entre pares de puntos según la distancia que los separa.

Ajuste del variograma teórico:

Modelamos esa curva con una función matemática (esférica, exponencial, gaussiana, etc.).

Predicción geoestadística:

Calculamos la superficie interpolada usando los pesos óptimos derivados del variograma.

# 1. Calculamos el variograma empírico a partir de las muestras
v_emp_ok <- gstat::variogram(soc ~ 1, data = samples)
# Graficamos el variograma empírico
plot(v_emp_ok, main = "Variograma Empírico de SOC")

# 2. Ajustamos automáticamente el modelo teórico
v_mod_ok <- automap::autofitVariogram(soc ~ 1, as(samples, "Spatial"))
# Graficamos el modelo de variograma ajustado sobre los puntos empíricos
plot(v_mod_ok)

Análisis: El variograma muestra cómo cambia la similitud de los valores de SOC a medida que aumenta la distancia entre los puntos. El nugget representa la variabilidad inicial o el error de medición, el sill indica el valor máximo de la variabilidad y el range corresponde a la distancia hasta la cual los puntos mantienen una relación espacial. Estos parámetros son importantes para ajustar el modelo de Kriging y realizar la interpolación.

# Visualizar la estructura del modelo ajustado
v_mod_ok$var_model
##   model    psill    range
## 1   Nug 89.02255 0.000000
## 2   Sph 49.99555 6.728419

En este paso se incorporó el modelo de variograma al objeto g2 y se utilizó para realizar la interpolación mediante Kriging Ordinario sobre la misma rejilla empleada anteriormente con el método IDW.

# 1. Creamos el modelo gstat incluyendo el variograma teórico
g2 <- gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)

# 2. Realizamos la predicción por Kriging Ordinario sobre la grilla continua
z2 <- predict(g2, stars.rrr)
## [using ordinary kriging]
# Manejo de NA y cambio de nombre
a <- which(is.na(z2[[1]]))
z2[[1]][a] = 0.0
names(z2) = "soc"

Visualización del resultado de Kriging Ordinario

Para visualizar la superficie interpolada por Kriging y comparar visualmente con el IDW:

paleta_ok <- colorNumeric(
  palette = c("orange", "yellow", "cyan", "green"), 
  domain = range(z2$soc, na.rm = TRUE), 
  na.color = "transparent"
)

# Mapa de la interpolación OK
m_ok <- leaflet() %>% 
  addTiles() %>% 
  leafem::addGeoRaster(
    z2, 
    opacity = 0.7, 
    colorOptions = colorOptions(
      palette = c("orange", "yellow", "cyan", "green"), 
      domain = range(z2$soc, na.rm = TRUE)
    )
  ) %>% 
  addCircleMarkers(
    data = samples, 
    radius = 1.5, 
    label = ~soc, 
    color = ~paleta_ok(soc), 
    fillOpacity = 1, 
    stroke = FALSE
  ) %>% 
  addLegend(
    "bottomright", 
    pal = paleta_ok, 
    values = z2$soc, 
    title = "OK SOC interpolation [g/kg]"
  )
## Warning in paleta_ok(soc): Some values were outside the color scale and will be
## treated as NA
## Warning in paleta_ok(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_ok

6. Evaluación de los resultados.

6.1 Evaluación cualitativa

Al comparar los dos métodos de interpolación se observan diferencias en las superficies generadas. El método IDW produce cambios más marcados alrededor de los puntos de muestreo, formando patrones circulares. En cambio, el Kriging Ordinario (OK) genera una superficie más suave y continua, ya que considera la relación espacial entre los datos mediante el variograma.

colores <- colorOptions(
  palette = c("orange", "yellow", "cyan", "green"), 
  domain = 10:100, 
  na.color = "transparent"
)
m_comp <- leaflet() %>%
  addTiles() %>%
  addGeoRaster(z1, colorOptions = colores, opacity = 0.8, group = "IDW") %>%
  addGeoRaster(z2, colorOptions = colores, opacity = 0.8, group = "OK") %>%
  addLayersControl(
    overlayGroups = c("IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  addLegend(
    "bottomright", 
    pal = paleta, 
    values = z1$soc,
    title = "Soil organic carbon [%]"
  )
## 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_comp

6.2 Validación Cruzada

Después de obtener las superficies de interpolación con los métodos IDW y Kriging Ordinario, se realizó una validación cruzada para evaluar cuál ofrece mejores resultados. Este procedimiento consiste en retirar un punto de muestreo, predecir su valor utilizando los demás datos y repetir el proceso con todos los puntos. Finalmente, se comparan los valores observados con los valores predichos para determinar la precisión de cada método de interpolación.

RMSE para IDW: 11.7394

RMSE para Ordinary Kriging (OK): 11.7230

7. Guardar las salidas de interpolación

write_stars(
  z1, dsn = "C:/Users/suare/OneDrive/Documentos/GB2/Proyecto1/Datos/IDW_soc_magdalena.tif", layer = 1
)

write_stars(
  z2, dsn = "C:/Users/suare/OneDrive/Documentos/GB2/Proyecto1/Datos/OK_soc_magdalena.tif", layer = 1
)

8. Conclusiones

Durante el desarrollo de esta práctica comprendí la diferencia entre los métodos de interpolación basados únicamente en la distancia y aquellos que consideran la relación espacial entre los datos. Aunque el proceso fue exigente, especialmente al preparar las mallas y ajustar el variograma, permitió entender mejor el funcionamiento de cada método. Al comparar la Distancia Inversa Ponderada (IDW) y el Kriging Ordinario (OK) para estimar el Carbono Orgánico del Suelo (SOC) en el departamento del Magdalena, se observó que el Kriging obtuvo un desempeño ligeramente superior, con un RMSE de 11.7230, frente al 11.7394 obtenido con IDW. Sin embargo, los resultados no abarcaron completamente el departamento, ya que los puntos de muestreo disponibles se concentraban en algunas zonas y no cubrían todo el territorio, lo que limitó la interpolación en las áreas sin datos.

9. Referencias

Lizarazo, I. 2025. Spatial interpolation. Available at: https://rpubs.com/ials2un/Spatial_Interpolation