R Notebook

1. Introducción

Este cuaderno ilustra dos técnicas de interpolación espacial aplicadas a datos de carbono orgánico del suelo (SOC) a una profundidad de 15-30 cm en el departamento del Putumayo. Se emplean la técnica de ponderación de distancia inversa (IDW), que es determinista, y el kriging ordinario (OK), que es probabilístico. Ambas metodologías permiten generar superficies continuas a partir de datos puntuales obtenidos de SoilGrids a una resolución de 250 m.

El análisis y procesamiento de estos datos siguen la estructura de un cuaderno guía proporcionado por el profesor de Geomática Básica. Sin embargo, este documento amplía la documentación del código y proporciona una interpretación de los resultados obtenidos. El cuaderno guia ejecutado en R,se encuentra aquí

2. Preparación del espacio de trabajo

Antes de comenzar con la interpolación espacial en el departamento del Putumayo, es fundamental organizar nuestro entorno de trabajo en R. Esto incluye limpiar el espacio de trabajo y cargar las bibliotecas necesarias.

2.1 ¿Por qué es importante limpiar el espacio de trabajo?

Al limpiar el entorno, eliminamos variables y objetos creados en sesiones anteriores que podrían interferir con el análisis actual. Esto ayuda a:

Evitar conflictos entre objetos con nombres similares. Liberar memoria y mejorar el rendimiento de R. Garantizar reproducibilidad, asegurando que el código funcione desde cero sin dependencias ocultas. Podemos limpiar el espacio de trabajo con:

rm(list=ls())

2.2 ¿Por qué es importante cargar las bibliotecas?

Las bibliotecas contienen funciones especializadas que permiten manejar datos geoespaciales y realizar análisis avanzados. Puntualmente, para el desarrollo de este cuaderno, se utilizaron:

• Lectura y escritura de datos geoespaciales: terra (datos ráster) y sf (datos vectoriales).
• Interpolación espacial: gstat y automap.
• Visualización de resultados: ggplot2, leaflet y tmap.

Antes de continuar, es importante asegurarse de tenerlas instaladas desde la consola de R (no en este cuaderno). Luego, cárgarlas con:

library(terra)

terra 1.8.21

library(sf)

Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2()
is TRUE

library(sp)
library(stars)

Loading required package: abind

library(gstat)
library(automap)
library(RColorBrewer)
library(leaflet)
library(leafem)

3. Lectura de datos de entrada

Como se mencionó en la introducción, las muestras utilizadas corresponden a valores de SOC a una profundidad de 15-30 cm, obtenidos de SoilGrids 250 m a través de un cuaderno previo, . Estos datos fueron almacenados en formato GeoPackage dentro del directorio de datos, siguiendo la nomenclatura soc_{putu}.gpkg.

Antes de continuar con el análisis, es importante verificar que el archivo se encuentra disponible y listo para su uso.

list.files(path="DATOS", pattern = "*.gpkg")

[1] "Putumayo.gpkg" "soc_putu.gpkg"

3.1 Lectura de datos de entrada

summary(samples)

hist(samples$soc)

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

pal <- colorNumeric(
  c("#E1F5C4", "#EDE574", "#F9D423", "#FC913A", "#FF4E50"),
  # colors depend on the count variable
  domain = samples$soc,
  )

munic <- munic[sf::st_geometry_type(munic) %in% c("POLYGON", "MULTIPOLYGON"), ]

# This a simple visualization
# Change the code to suit your preferences
leaflet() %>%
  addPolygons(
    data = munic,
    color = "gray",
    # set the opacity of the outline
    opacity = 1,
    # set the stroke width in pixels
    weight = 1,
    # set the fill opacity
    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")

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

#Extensión  -77.1833333329999931,-0.5583333329999991 : -73.8416666669999984,1.4666666669999999
#Anchura    401
#Altura 243

#(rrr <- terra::rast(xmin=-74.55524, xmax=-72.38985, ymin=5.708308, ymax=8.145474, nrows=370, ncols = 329, vals = 1, crs = "epsg:4326"))
(rrr <- terra::rast(xmin=-77.183, xmax=-73.842, ymin=-0.558, ymax=1.47, nrows=243, ncols = 401, vals = 1, crs = "epsg:4326"))

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

## [inverse distance weighted interpolation]
## running this code takes several minutes
z1 = predict(g1, stars.rrr)

# dealing with NA values
a <- which(is.na(z1[[1]]))
z1[[1]][a] = 0.0

z1

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

## [inverse distance weighted interpolation]
## running this code takes several minutes
z1 = predict(g1, stars.rrr)

z1

# dealing with NA values
a <- which(is.na(z1[[1]]))
z1[[1]][a] = 0.0

z1

names(z1) = "soc"

# change this chunk to meet your preferences
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  # Print the map

range(z1$soc, na.rm = TRUE)

# Obtener el rango de valores reales
soc_range <- range(z1$soc, na.rm = TRUE)

# Ajustar la paleta al rango correcto
paleta <- colorNumeric(palette = c("orange", "yellow", "cyan", "green"), 
                       domain = soc_range, na.color = "transparent")

leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      z1,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("orange", "yellow", "cyan", "green"), 
                                  domain = soc_range)  # Aquí aplicamos el rango real
    ) %>%
   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 [%]"
    )

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

plot(v_emp_ok)

#make sure you understand the parameters needed to run this fitting function
v_mod_ok = automap::autofitVariogram(soc ~ 1, as(samples, "Spatial"))

plot(v_mod_ok)

v_mod_ok$var_model

## [using ordinary kriging]
## This takes several minutes 
## (or even ages depending on your raster cell size and your computer resources)
g2 = gstat(formula = soc ~ 1, model = v_mod_ok$var_model, data = samples)
z2= predict(g2, stars.rrr)

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

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  # Print the map

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  # Print the map