1 Introducción

Este notebook presenta un análisis de interpolación espacial del carbono orgánico del suelo (SOC) en el municipio del Meta, Colombia. El análisis se realiza siguiendo la guía académica del profesor Iván Lizarazo, en la Facultad de Ciencias Agrarias de la Universidad Nacional de Colombia. Se aplican los métodos de interpolación IDW (Distancia Inversa Ponderada) y Kriging Ordinario.

2 Carga de librerías

library(sf)
library(gstat)
library(leaflet)
library(leafem)
library(terra)
library(stars)
library(automap)

3 Lectura y preparación de datos

3.1 Lectura de archivos espaciales

municipio <- st_read("C:/Users/vidam/Downloads/Municipio Meta.gpkg")
## Reading layer `mgn_adm_mpio_grafico' from data source 
##   `C:\Users\vidam\Downloads\Municipio Meta.gpkg' using driver `GPKG'
## Simple feature collection with 29 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.93335 ymin: 1.617732 xmax: -71.07753 ymax: 4.899101
## Geodetic CRS:  WGS 84
samples <- st_read("C:/Users/vidam/Downloads/socmeta.gpkg")
## Reading layer `puntos_vectoriales' from data source 
##   `C:\Users\vidam\Downloads\socmeta.gpkg' using driver `GPKG'
## Simple feature collection with 2339260 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.9318 ymin: 1.621266 xmax: -71.0795 ymax: 4.897252
## Geodetic CRS:  WGS 84

3.2 Transformación y filtrado

samples <- st_transform(samples, crs = st_crs(municipio))
samples_meta <- samples[municipio, ]
set.seed(123)
samples <- samples[sample(1:nrow(samples), 2000), ]

4 Exploración de datos

4.1 Estadísticas básicas

nrow(samples)
## [1] 2000
summary(samples$SOC.PUNTOS)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     0.0   107.0   104.9   155.0   871.0

4.2 Histograma

hist(samples$SOC.PUNTOS, main = "Distribución de SOC (%)",
     xlab = "SOC", col = "lightblue")
Distribución del SOC

Distribución del SOC

Interpretación: El histograma muestra la frecuencia de valores de SOC en los puntos muestreados. Se observa una distribución moderadamente sesgada, con valores más frecuentes entre 20% y 40% de SOC.

4.3 Mapa de puntos

samples$SOC.PUNTOS <- round(samples$SOC.PUNTOS, 2)
pal <- colorNumeric(
  palette = c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),
  domain = samples$SOC.PUNTOS
)

leaflet() %>%
  addProviderTiles("OpenStreetMap") %>%
  addPolygons(data = municipio, color = "gray", weight = 1, fillOpacity = 0.2) %>%
  addCircleMarkers(data = samples, radius = 1.5, label = ~SOC.PUNTOS,
                   color = ~pal(SOC.PUNTOS), fillOpacity = 1, stroke = FALSE) %>%
  addLegend("bottomleft", pal = pal, values = samples$SOC.PUNTOS,
            title = "SOC (%)", opacity = 0.9)

Interpretación: Este mapa muestra la distribución espacial de los puntos de muestreo y su concentración de SOC. Se observan zonas con mayor acumulación de carbono en áreas específicas del municipio.

5 Interpolación con IDW

5.1 Preparación del raster

g1 <- gstat(formula = SOC.PUNTOS ~ 1, data = samples)
ext <- st_bbox(municipio)
rrr <- terra::rast(
  xmin = ext["xmin"], xmax = ext["xmax"],
  ymin = ext["ymin"], ymax = ext["ymax"],
  nrows = 300, ncols = 300, vals = 1,
  crs = "EPSG:4326"
)
stars.rrr <- stars::st_as_stars(rrr)

5.2 Interpolación y visualización

z1 <- predict(g1, stars.rrr)
## [inverse distance weighted interpolation]
z1[[1]][is.na(z1[[1]])] <- 0
names(z1) <- "soc"

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

leaflet() %>%
  addTiles() %>%
  addGeoRaster(z1, opacity = 0.7,
               colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"),
                                           domain = 10:100)) %>%
  addCircleMarkers(data = samples, radius = 1.5, label = ~SOC.PUNTOS,
                   color = ~paleta(SOC.PUNTOS), fillOpacity = 1, stroke = FALSE) %>%
  addLegend("bottomright", pal = paleta, values = z1$soc,
            title = "IDW SOC interpolation [%]")
## 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.PUNTOS): Some values were outside the color scale and
## will be treated as NA
## Warning in paleta(SOC.PUNTOS): 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

Interpretación: La interpolación IDW refleja la influencia de los puntos cercanos sobre los valores estimados. Las zonas cercanas a concentraciones altas de SOC muestran estimaciones elevadas, pero sin continuidad espacial sofisticada.

6 Interpolación con Kriging Ordinario

6.1 Cálculo del variograma

v_emp_ok <- variogram(SOC.PUNTOS ~ 1, data = samples)
plot(v_emp_ok, main = "Variograma empírico (SOC)")

Interpretación: El variograma empírico describe la dependencia espacial entre los puntos. Su forma indica que existe correlación espacial, justificada para aplicar Kriging.

6.2 Ajuste del modelo y Kriging

v_mod_ok <- automap::autofitVariogram(SOC.PUNTOS ~ 1, as(samples, "Spatial"))
plot(v_mod_ok)

v_mod_ok$var_model
##   model      psill    range kappa
## 1   Nug   569.2035  0.00000   0.0
## 2   Ste 20622.0525 78.07375   0.4
g2 <- gstat(formula = SOC.PUNTOS ~ 1,
            model = v_mod_ok$var_model,
            data = samples)

z2 <- predict(g2, stars.rrr)
## [using ordinary kriging]
z2[[1]][is.na(z2[[1]])] <- 0
names(z2) <- "soc"

leaflet() %>%
  addTiles() %>%
  addGeoRaster(z2, opacity = 0.7,
               colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"),
                                           domain = 10:100)) %>%
  addCircleMarkers(data = samples, radius = 1.5, label = ~SOC.PUNTOS,
                   color = ~paleta(SOC.PUNTOS), fillOpacity = 1, stroke = FALSE) %>%
  addLegend("bottomright", pal = paleta, values = z2$soc,
            title = "OK SOC interpolation [%]")
## 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.PUNTOS): Some values were outside the color scale and
## will be treated as NA
## Warning in paleta(SOC.PUNTOS): 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

Interpretación: El Kriging considera la estructura espacial detectada por el variograma, resultando en una superficie más continua y realista, especialmente útil en análisis de suelos.

7 Comparación de métodos

7.1 Visualización comparativa

leaflet() %>%
  addTiles() %>%
  addGeoRaster(z1, colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"),
                                               domain = 10:100), opacity = 0.8, group = "IDW") %>%
  addGeoRaster(z2, colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"),
                                               domain = 10:100), opacity = 0.8, group = "OK") %>%
  addCircleMarkers(data = samples, radius = 1.5, label = ~SOC.PUNTOS,
                   color = ~paleta(SOC.PUNTOS), fillOpacity = 1, stroke = FALSE) %>%
  addLayersControl(
    overlayGroups = c("IDW", "OK"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  addLegend("bottomright", pal = paleta, values = z1$soc,
            title = "Soil Organic Carbon [%]")
## 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.PUNTOS): Some values were outside the color scale and
## will be treated as NA
## Warning in paleta(SOC.PUNTOS): 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

Interpretación: Comparando ambos métodos, se aprecia que el Kriging suaviza mejor las transiciones y ofrece interpolaciones más continuas, mientras que IDW tiende a formar “zonas” influenciadas por puntos cercanos.

8 Exportación de resultados

8.1 Guardar como archivos raster

write_stars(z1, dsn = "C:/Users/vidam/Downloads/IDW_SOC_Meta.tif", layer = 1)
write_stars(z2, dsn = "C:/Users/vidam/Downloads/OK_SOC_Meta.tif", layer = 1)

9 Conclusiones

  • Se implementaron y compararon dos métodos de interpolación espacial sobre los datos de SOC en el municipio del Meta: IDW y Kriging Ordinario.
  • IDW resultó útil para visualizaciones rápidas, mientras que Kriging ofreció una superficie más realista y continua.
  • Las diferencias observadas reflejan la importancia de seleccionar el método apropiado según el objetivo del análisis.
  • Los resultados exportados permiten su uso posterior en SIG y análisis agroambientales.

10 Referencias

  • Lizarazo, I. (2023). Notas de clase de Geomática. Facultad de Ciencias Agrarias, Universidad Nacional de Colombia.