A través de técnicas de análisis espacial, procesamiento de datos raster y visualización cartográfica, este trabajo busca caracterizar la distribución espacial de estas propiedades edáficas, identificar patrones regionales y establecer relaciones entre las variables analizadas. La relación C/N emerge como un indicador particularmente relevante, ya que refleja el equilibrio entre los procesos de descomposición de la materia orgánica y la disponibilidad de nutrientes, aspectos cruciales para la productividad agrícola y la salud del suelo.
Como paso inicial para este cuaderno, se limpia la memoria en nuestro trabajo para evitar la confusión en los datos.
rm(list=ls())
Para el trabajo se van a utilizar las siguientes librerias:
library(XML)
## Warning: package 'XML' was built under R version 4.5.2
library(sf)
## Warning: package 'sf' was built under R version 4.6.0
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(gdalUtilities)
##
## Adjuntando el paquete: 'gdalUtilities'
## The following object is masked from 'package:sf':
##
## gdal_rasterize
library(terra)
## Warning: package 'terra' was built under R version 4.6.0
## terra 1.8.70
library(dplyr)
##
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:terra':
##
## intersect, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.6.0
library(RColorBrewer)
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.5.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.1 ✔ readr 2.1.5
## ✔ ggplot2 4.0.0 ✔ stringr 1.5.2
## ✔ lubridate 1.9.4 ✔ tibble 3.3.0
## ✔ purrr 1.1.0 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks terra::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(tidyterra)
## Warning: package 'tidyterra' was built under R version 4.5.2
##
## Adjuntando el paquete: 'tidyterra'
##
## The following object is masked from 'package:stats':
##
## filter
**XML:*:Esta librería proporciona herramientas para analizar y generar XML (Lenguaje de Marcado Extensible) en R. Se utiliza para leer y manipular documentos XML y HTML, lo cual es útil para extraer datos de páginas web o trabajar con formatos de datos basados en XML.
**sf*:Es una librería para trabajar con datos espaciales vectoriales. Proporciona una forma estandarizada de representar objetos espaciales y permite realizar operaciones geométricas, manipulaciones y transformaciones de sistemas de referencia espacial (CRS).
**ggalutilities*: Incluye componentes adicionales para ggplot2, como geometrías y estadísticas extendidas. Incluye funciones para crear gráficos alternativos, por ejemplo, gráficos de barras con intervalos, splines, etc.
**terra*:Es una librería para el análisis de datos raster y vectoriales. Reemplaza y mejora la funcionalidad de raster, proporcionando un manejo más eficiente de grandes conjuntos de datos espaciales, como operaciones algebraicas, manipulación de capas y análisis espacial.
**dplyr*: ofrece un conjunto de verbos coherentes para resolver problemas comunes de manipulación de datos. Incluye funciones que permiten transformar y resumir datos de manera intuitiva.
**leaflet*: permite crear mapas interactivos. Permite visualizar datos espaciales en mapas base con capas interactivas, popups, leyendas y controles.
**RColorBrewer*: Ofrece paletas de colores predefinidas que son adecuadas para gráficos temáticos y mapas. Incluye paletas secuenciales, divergentes y cualitativas, que son útiles para representar datos de manera efectiva.
**tidyverse*: Es un meta-paquete que carga un conjunto de librerías diseñadas para ciencia de datos. Incluye dplyr, ggplot2, tidyr, readr, purrr, tibble, stringr, forcats, entre otras. Proporciona un ecosistema coherente para importar, transformar, visualizar y modelar datos.
SoilGrids genera mapas de suelos mediante modelos globales; si se está buscando información de suelos en niveles nacionales se debe comparar las predicciones de SoilGrids con mapas de suelos derivados de agencias nacionales de suelos. En Colombia.
Los mapas nacionales de suelos suelen elaborarse a partir de información más detallada sobre las características del suelo, lo que les permite ofrecer una mayor precisión que SoilGrids en las zonas donde tienen cobertura. Tanto la media como la mediana D.S. quantity pueden emplearse para estimar una propiedad del suelo en una celda específica. La media refleja el valor esperado, por lo que constituye una estimación no sesgada de dicha propiedad.
La relación carbono/nitrógeno (C/N) es un indicador fundamental de la calidad y del estado de descomposición de la materia orgánica del suelo. Este valor expresa la proporción entre el carbono —principal fuente de energía para los microorganismos— y el nitrógeno, que es esencial para la síntesis de proteínas y enzimas durante el proceso de descomposición. Cuando la relación C/N es alta, significa que hay mucho carbono en comparación con el nitrógeno, lo que provoca una descomposición más lenta, ya que los microorganismos deben inmovilizar nitrógeno del suelo para poder degradar la materia orgánica.
En cambio, una relación C/N baja indica una mayor disponibilidad relativa de nitrógeno, lo que favorece una descomposición más rápida y una liberación acelerada de nutrientes al suelo. Este equilibrio influye directamente en la fertilidad del suelo, la disponibilidad de nitrógeno para las plantas y la dinámica del carbono en los ecosistemas. Por ello, el control y monitoreo de la relación C/N es clave para comprender los procesos de ciclo de nutrientes y mantenimiento de la materia orgánica en distintos sistemas agrícolas y naturales.
Primero, definimos la URL correspondiente al sitio webDAD -servicio para acceder a datos de suelos de forma remota- de ISRIC -International Soil Reference and Information Centre-.
url = "https://files.isric.org/soilgrids/latest/data/"
Una vez obtenida la información se crean los siguientes objetos:
voi_soc <- "soc"
voi_nit <- "nitrogen"
depth <- "15-30cm"
quantile <- "mean"
variable_soc <- paste(url, voi_soc, sep = "")
variable_nit <- paste(url, voi_nit, sep = "")
A objetos -voi_soc <- “soc”- y -voi_nit <- “nitrogen”- les queremos atribuir la inforcaión del ISREC, por lo que primero con (depth <- “15-30cm”) le indicamos la profundidad del suelo a la cual se quieren obtener los datos. En este caso, es la capa de 15 a 30 centímetros de profundidad. Y con (quantile <- “mean”) Especifica el tipo de valor estadístico que quieres usar. “mean” significa que vas a descargar el valor promedio estimado de la propiedad del suelo. y por ultimo con la función “paste()” galvanizamos el url de las variables de interes (soc o nitrogen).
Ahora definimos las propiedades del suelo que queremos descargar:
layer_soc <- paste(variable_soc, depth, quantile, sep = "_")
layer_nit <- paste(variable_nit, depth, quantile, sep = "_")
Nuevamente, con la función “paste()”unimos, en este caso la profundidad de 15-30cm, a las variables antes realizadas.
vrt_layer_soc <- paste(layer_soc, ".vrt", sep = "")
vrt_layer_nit <- paste(layer_nit, ".vrt", sep = "")
Y por ultimo, con la misma función, ahora unimos todo en un “.vrt”. Virtual Raster Tile, un formato que sirve como un enlace remoto al archivo raster (no descarga todavía los datos, solo los referencia).
Ahora, necesitamos leer un GeoPackage de un departamento descargado previamente de DANE. Para ello, utilizaremos la biblioteca sf:
Para ello listamos los objetos en la capeta de trabajo:
list.files()
## [1] "cund_CN_ratio.tif" "data"
## [3] "MagnaCund_munic_dane.gpkg" "nitrogen_igh_15_30.tif"
## [5] "soc_igh_15_30.tif" "spatial_interpolation.Rmd"
## [7] "WGS84_Cund_munic.qgz" "WGS84Cund_munic.gpkg"
Y seleccionamos el archivo GeoPackage de nuestro departamento:
cund <- st_read("WGS84Cund_munic.gpkg") %>%
st_as_sf() %>%
mutate(MPIO_CNMBR_MAY <- str_to_upper(MPIO_CNMBR))
## Reading layer `WGS84Cund_munic' from data source
## `E:\GEOMATICA\GB\proyect_5\WGS84Cund_munic.gpkg' using driver `GPKG'
## Simple feature collection with 116 features and 53 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.89063 ymin: 3.730129 xmax: -73.05256 ymax: 5.837258
## Geodetic CRS: WGS 84
cund %>% select(DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR) %>%
head()
## Simple feature collection with 6 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.8191 ymin: 4.175736 xmax: -74.23773 ymax: 4.952534
## Geodetic CRS: WGS 84
## DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR geom
## 1 25 019 ALBÁN MULTIPOLYGON (((-74.47852 4...
## 2 25 035 ANAPOIMA MULTIPOLYGON (((-74.55201 4...
## 3 25 040 ANOLAIMA MULTIPOLYGON (((-74.43043 4...
## 4 25 053 ARBELÁEZ MULTIPOLYGON (((-74.42848 4...
## 5 25 086 BELTRÁN MULTIPOLYGON (((-74.76139 4...
## 6 25 095 BITUIMA MULTIPOLYGON (((-74.51062 4...
Queremos atribuir los atributos de del geopackage a una nueva variable “cund”. Y la encadenamos a “st_as_sf()” que asegura que el objeto resultante conserve la estructura espacial compatible con sf.con las función “mutate()” de dplyr (parte del tidyverse), queremos crear o modificar columnas en un data frame, en este caso el nombre de los municipios; y con “str_to_upper()” hacemos que aparezcan en mayúsculas. Posteriomente, con la función “select()” sirve para seleccionar columnas específicas del conjunto de datos.
Como se puede apreciar en el archivo utiliza un CRS: “WGS84 (EPSG:4326)” es el sistema de referencia geodésico global que usan la mayoría de los datos satelitales, incluyendo los de SoilGrids (ISRIC), MODIS, Sentinel, Landsat, etc.
Con el siguiente codigo:
igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
cund_igh <- st_transform(cund, igh)
Se busca definir las coordenadas para nuestro “Homolosine Projectión”. vamos a atribuirle a un objeto “igh” con (+proj=igh) para indicarle que se usará la proyección “Interrupted Goode Homolosine”, una proyección cartográfica pseudocilíndrica que minimiza la distorsión en las áreas continentales mediante interrupciones en los océanos.Los parámetros (+lat_0=0 y +lon_0=0) establecen el centro de la proyección en el ecuador (latitud 0) y el meridiano de Greenwich (longitud 0).(+datum=WGS84) especifica el datum geodésico WGS84, el sistema de referencia global más utilizado, que con (+units=m) establece las unidades en metros .
“st_transform(cund, igh)” toma el objeto (cund), que contiene datos geográficos como polígonos, puntos o líneas del departamento; y lo reproyecta desde su sistema de coordenadas original al sistema Interrupted Goode Homolosine.
Junto a este codigo:
(bbox <- st_bbox(cund_igh))
## xmin ymin xmax ymax
## -8344613.5 415236.1 -8142417.0 649800.5
necesitamos una caja delimitadora (bounding box) que contenga nuestro objeto “cund_igh”. Un rectángulo mínimo que contiene completamente todos los elementos geométricos de tu dataset espacial. con el argumento “(st_bbox())” lo que se busca es extraer los limites espaciales del objeto y se lo atribuimos a la bounding box que necesitamos.
ulx = bbox$xmin
uly = bbox$ymax
lrx= bbox$xmax
lry = bbox$ymin
(bb <- c(ulx, uly, lrx, lry))
## xmin ymax xmax ymin
## -8344613.5 649800.5 -8142417.0 415236.1
adicionalmente, lo que se busca con el codigo anterior es extrae y reorganizar las coordenadas de la bounding box en un formato específico, particularmente útil para ciertas aplicaciones geoespaciales.”(ulx = bbox\(xmin y lrx = bbox\)xmax )” toman las coordenadas minimas y maximas de x, respectivamente y las asignan a “ulx” - upper left- y “lrx” - lower right- ; y “(lry = bbox\(ymin y lry = bbox\)ymin)” toman las coordenadas minimas y maximasde y, respectivamente y las asignan a “lry” - lower right- y “lry” - upper left-.
Ahora, lo que se busca es descargar un raster de la variablo, o digase propiedad, de interes -nitrogeno y carbono orgánico del suelo-:
sg_url="/vsicurl/https://files.isric.org/soilgrids/latest/data/"
soc_data <- 'soc/soc_15-30cm_mean.vrt'
soc_file <- "soc_igh_15_30.tif"
nitrogen_data <- 'nitrogen/nitrogen_15-30cm_mean.vrt'
nitrogen_file <- "nitrogen_igh_15_30.tif"
Se quiere con “sg_url=”/vsicurl/https://files.isric.org/soilgrids/latest/data/“” generar una variable unida a un URL base donde están almacenados los datos de SoilGrids, una base mundial de propiedades del suelo (del ISRIC). El propocito del protocolo de comunicaciÓn (https:) es guiar al servidor donde están los archivos .vrt (raster virtuales) de propiedades del suelo; que junato a “/vsicurl/” le dice a R que acceda directamente a un archivo remoto a través de HTTP sin descargarlo completamente.
“soc_data <- ‘soc/soc_15-30cm_mean.vrt’” y “nitrogen_data<-‘nitrogen/nitrogen_15-30cm_mean.vrt’” buscan aguardar la ruta dentro del servidor del Contenido de carbono orgánico del suelo (SOC) Y Nitrógeno del suelo, de unos 15 a 30 cm de profundidad. y los archivos raster seluntantes los asociamos a “soc_file” y “nitrogen_file” respectivamente.
En conjunto:
#gdal_translate(paste0(sg_url, soc_data), soc_file,
# tr = c(250, 250),
# projwin = bb,
#projwin_srs = igh)
#gdal_translate(paste0(sg_url, nitrogen_data), nitrogen_file,
#tr = c(250, 250),
#projwin = bb,
#projwin_srs = igh)
Junto a la libreria gdalUtilities, queremos extraer (recortar) y convertir una parte de un raster remoto (.vrt en este caso) en un nuevo archivo .tif local.”paste0(sg_url, datos)” une el URL base con el dataset remoto tanto del SOC, como del nitrogeno. y con el sufijo file, se la da el nombre del archivo de salida local (el .tif que vas a guardar); “tr = c(250, 250)” Define la resolución espacial (en metros o unidades del CRS destino). En este caso, crea un raster con píxeles de 250x250 m. “projwin = bb” representa el bounding box (ventana espacial) que define el recorte y bb suele ser un vector con los límites en coordenadas del CRS destino; que con “projwin_srs = igh” Define el sistema de referencia de coordenadas (CRS) en el que está expresada esa ventana (bb).
(cund_soc <- terra::rast(soc_file)/10)
## class : SpatRaster
## size : 940, 810, 1 (nrow, ncol, nlyr)
## resolution : 250, 250 (x, y)
## extent : -8344750, -8142250, 415000, 650000 (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine
## source(s) : memory
## varname : soc_igh_15_30
## name : soc_igh_15_30
## min value : 8.9
## max value : 322.9
(cund_nitrogen <- terra::rast(nitrogen_file)/100)
## class : SpatRaster
## size : 940, 810, 1 (nrow, ncol, nlyr)
## resolution : 250, 250 (x, y)
## extent : -8344750, -8142250, 415000, 650000 (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine
## source(s) : memory
## varname : nitrogen_igh_15_30
## name : nitrogen_igh_15_30
## min value : 0.87
## max value : 13.54
(terra::rast(soc_file) y terra::rast(nitrogen_file)) cargan los archivos TIFF descargados como objetos SpatRaster del paquete terra; convirtiendo los datos del disco en objetos manipulables en R.
para el SOC se dividen todos los valores del raster por 10Convierte de dg/kg (decigramos/kilogramo) a g/kg (gramos/kilogramo) –> “136.3 dg/kg ÷ 10 = 13.63 g/kg”. para el nitrogeno se Divide todos los valores del raster por 100, Convierte de cg/kg (centigramos/kilogramo) a g/kg (gramos/kilogramo) –> 50 cg/kg ÷ 100 = 0.5 g/kg
ahora generamos un histograma para cada variable de interes, (terra::hist()) crea un histograma de los valores de un objeto ráster en este caso del Carbono Orgánico del SUelo y el nitrogeno, para “xlab” en cada histograma, se encuentra la cantidad de soc y nitrogeno en g/kg y e la “ylab” se encontraria su frecuencia en los pexeles del raster. los atributos “col”y “border” dan color tanto al histograma como a las barras.
terra::hist(cund_soc,
main = "Valor Modal Carbono Orgánico del Suelo (SOC)\nProfundidad: 15-30 cm - Departamento de Cundinamarca",
xlab = "Carbono Orgánico (g/kg)",
ylab = "Frecuencia de Pixeles",
col = "steelblue",
border = "white")
terra::hist(cund_nitrogen,
main = "Valor Modal Nitrógeno Total del Suelo\nProfundidad: 15-30 cm - Departamento de Cundinamarca",
xlab = "Nitrógeno Total (g/kg)",
ylab = "Frecuencia de Pixeles",
col = "darkgreen",
border = "white")
El histograma correspondiente al carbono orgánico del suelo muestra que la mayor parte de los valores se concentran entre 20 y 60 g/kg, con un pico alrededor de los 30–40 g/kg. Esta distribución evidencia la existencia de un número reducido de zonas con contenidos altos.Este comportamiento sugiere que la mayor parte de los suelos del departamento de Cundinamarca presentan niveles medios de materia orgánica, mientras que algunas áreas particulares —probablemente asociadas a condiciones de mayor humedad, acumulación de residuos vegetales o prácticas agrícolas conservacionistas— muestran valores más elevados.
El promedio estimado de 39,6 g/kg indica que, en general, los suelos poseen una reserva orgánica significativa, la cual contribuye a mejorar la estructura del suelo, aumentar su capacidad de retención de agua y potenciar la actividad biológica. Estos valores son característicos de suelos con un buen grado de fertilidad natural o con una cobertura vegetal permanente que favorece la incorporación continua de materia orgánica.
En el histograma del nitrógeno total se observa una mayor concentración de valores en el rango de 1,5 a 4 g/kg, con una disminución progresiva hacia valores más altos. El máximo se sitúa cerca de los 2–3 g/kg, lo que sugiere que la mayoría de los suelos presentan contenidos bajos a moderados de nitrógeno disponible. La distribución revela la presencia de pocas áreas con valores más altos (por encima de 8 g/kg), que probablemente coinciden con zonas de mayor contenido de materia orgánica o donde la mineralización del nitrógeno es más activa.
El promedio general, cercano a 2,9 g/kg, refleja que la disponibilidad de nitrógeno podría ser un factor limitante para la fertilidad en varios sectores del departamento, especialmente en aquellos con alta relación C/N o donde la descomposición de la materia orgánica es más lenta.
También pedimos un resumen de las dos variables;
summary(cund_soc)
## soc_igh_15_30
## Min. : 9.10
## 1st Qu.: 25.90
## Median : 38.30
## Mean : 39.65
## 3rd Qu.: 51.40
## Max. :209.80
## NA's :1810
summary(cund_nitrogen)
## nitrogen_igh_15_30
## Min. : 0.910
## 1st Qu.: 1.660
## Median : 3.050
## Mean : 2.935
## 3rd Qu.: 3.860
## Max. :11.620
## NA's :1810
Min. : 9.10 # Valor mínimo encontrado 1st Qu. : 25.90 # Primer cuartil (25% de los datos son ≤ este valor) Median : 38.30 # Mediana (50% de los datos son ≤ este valor) Mean : 39.65 # Promedio de todos los valores 3rd Qu. : 51.40 # Tercer cuartil (75% de los datos son ≤ este valor) Max. :209.80 # Valor máximo encontrado NA’s :1810 # Número de píxeles sin datos (missing values)
Min. : 0.910 # Valor mínimo 1st Qu. : 1.660 # Primer cuartil Median : 3.050 # Mediana Mean : 2.935 # Promedio 3rd Qu. : 3.860 # Tercer cuartil Max. :11.620 # Valor máximo NA’s :1810 # Píxeles sin datos
y cambiamos los nombres de los atributos:
(names(cund_soc) <- "soc")
## [1] "soc"
(names(cund_nitrogen) <- "nitrógeno")
## [1] "nitrógeno"
se introducen los valores se (SOC Y Nitrógeno) dentro de sus respectivas variables.
valores_soc <- values(cund_soc, na.rm = TRUE)
valores_nitrógeno <- values(cund_nitrogen, na.rm = TRUE)
Para los mapas, La función “colorNumeric()” crea una escala continua de colores basada en valores numéricos (en este caso, las concentraciones de carbono orgánico del suelo o el Nitógeno).El argumento “palette()” define los colores de transición; que junto a “domain = valores_soc”, especifica el rango de valores que cubrirá la paleta, por ultimo “na.color = transparent” evita que los valores nulos (NA) aparezcan como manchas de color, dejándolos invisibles en el mapa.
“leaflet(cund)” inicializa un mapa interactivo centrado en el objeto espacial cund, que contiene los límites municipales de Cundinamarca, que junto a “addTiles()” añade un fondo cartográfico, permitiendo contextualizar el mapa en el espacio geográfico. “setView(lng = -73.97, lat = 4.78, zoom = 9)” establece el punto de enfoque inicial y el nivel de zoom, siendo las coordenadas de longitud y latitud ubican el centro del mapa en Cundinamarca y un nivel de zoom 9 permite una vista departamental, mostrando claramente los límites municipales sin perder detalle regional.
“addPolygons()” dibuja los límites municipales contenidos en el objeto del departamento (cund) y junto a (color = “black”) define el borde de los municipios en negro para una mejor delimitación visual, y se le atribuyén ciertos parametros para su visualización; además la función (popup = paste(“Municipio:”, cund$MPIO_CNMBR)) genera una ventana emergente que muestra el nombre del municipio.”addRasterImage()” inserta la capa raster cund_soc o´cund_nitrogen, que contiene los valores de las variables de interes. (colors =) aplica la paleta de colores creada anteriormente, traduciéndolos en un gradiente visual de tonos.
opacity = 0.8 define el nivel de transparencia del raster, permitiendo visualizar simultáneamente la información de fondo (tiles y límites municipales). “addLegend()” añade una leyenda interpretativa al mapa.(pal =) asegura que la escala de colores de la leyenda coincida con la del raster. (values =) define el rango de valores a representar; y (title=) especifica el encabezado de la leyenda, aclarando a que corresponden los valores.
soc_sunset <- colorNumeric(
palette = c("#2D004B", "#7A1E6A", "#E35205", "#FDB863", "#FFF2B2"),
domain = valores_soc,
na.color = "transparent")
leaflet::leaflet(cund) %>%
addTiles() %>%
setView(lng = -73.97, lat = 4.78, zoom = 9) %>%
addPolygons(
color = "black",
weight = 1.0,
smoothFactor = 1.0,
opacity = 0.5,
fillOpacity = 0.10,
popup = paste("Municipio:", cund$MPIO_CNMBR)
) %>%
addRasterImage(
cund_soc,
colors = soc_sunset,
opacity = 0.8
) %>%
addLegend(
pal = soc_sunset,
values = valores_soc,
title = "Carbono Orgánico del Suelo (SOC) [g/kg]"
)