Interpolación espacial del Meta

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 utilizan métodos como IDW y Kriging Ordinario.

2 Carga de librerías

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

# Leer shapefile del municipio
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

# Leer puntos de SOC recortados al municipio
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

# Asegurar sistema de referencia espacial común
samples <- st_transform(samples, crs = st_crs(municipio))

# Filtrar puntos que están dentro del municipio
samples_meta <- samples[municipio, ]
# Reducir a 2000 puntos aleatoriamente
set.seed(123)
samples <- samples[sample(1:nrow(samples), 2000), ]

# Revisar estructura y resumen
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

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

# Redondear valores
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)

# Crear objeto gstat para IDW
g1 <- gstat(formula = SOC.PUNTOS ~ 1, data = samples)

# Crear raster vacío en el bbox del municipio
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"
)

# Convertir raster a stars
stars.rrr <- stars::st_as_stars(rrr)

# Interpolación IDW
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() %>%
  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.PUNTOS,
                   color = ~paleta(SOC.PUNTOS), fillOpacity = 1, stroke = FALSE) %>%
  addLegend("bottomright", pal = paleta, values = z1$soc,
            title = "IDW SOC interpolation [%]")

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

# Ajustar modelo automático
v_mod_ok <- automap::autofitVariogram(SOC.PUNTOS ~ 1, as(samples, "Spatial"))
plot(v_mod_ok)

# Modelo ajustado
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

# Kriging
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() %>%
  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.PUNTOS,
                   color = ~paleta(SOC.PUNTOS), fillOpacity = 1, stroke = FALSE) %>%
  addLegend("bottomright", pal = paleta, values = z2$soc,
            title = "OK SOC interpolation [%]")

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 [%]")

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)