1 1. Análisis Espacial de Carbono Orgánico, Nitrógeno y su Relación C/N en Suelos de Cundinamarca mediante SoilGrids y R.

1.1 2. Introduución.

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.

1.2 3. Configuración.

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.

1.3 4. Definir las variables de interés.

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.

1.3.1 4.1. relación carbono/nitrógeno (C/N) y M.O.

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.

1.3.2 4.2. Descargar un archivo TIFF correspondiente a la región de interés (ROI).

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).

1.4 5. Descargar un archivo TIFF para una región de interés (ROI).

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.

1.4.1 5.1 Definición de la Bounding Box y Reorganización de Coordenadas.

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-.

1.4.2 5.2 Descarga y Recorte de Rasters de SoilGrids.

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

1.5 6. Exploración de datos.

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")

1.5.1 6.2. Distribución del Carbono Orgánico del Suelo (SOC).

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.

1.5.2 6.3. Distribución del Nitrógeno Total del Suelo (N).

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

1.5.3 6.4. ANÁLISIS CARBONO ORGANICO DEL SUELO -RASTER_:soc_jgh_15_30:

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)

1.5.4 6.5. ANÁLISIS NITROGENO RASTER:Anitrogen_jgh_15_30

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"

1.6 7. Mapas Representando las variables de interés.

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.

1.6.1 7.1 Mapa Carbono Orgánico del Suelo.

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

El mapa del SOC muestra la distribución espacial del contenido de carbono orgánico (g/kg) en los suelos del departamento de Cundinamarca, a profundidad de 15–30 cm. Los limites de la paleta de colores representa la variación del carbono Orgánico. Los Tonos morados y violetas representan valores bajos de carbono orgánico, asociados a suelos más degradados o con menor acumulación de materia orgánica.Por otra parte, los Tonos anaranjados, dorados y amarillos indican valores altos de SOC, característicos de suelos con buena cobertura vegetal, mayor humedad y menor alteración al medio del suelo, donde la materia orgánica se acumula más fácilmente.

Se observa que las zonas con mayores concentraciones de carbono (dorado) se localizan en sectores montañosos o de piemonte (como la Cordillera Oriental y áreas altas del occidente del departamento), mientras que los valores bajos (morado) predominan en las zonas centrales y planas, donde existe mayor actividad agrícola.El histograma del SOC complementa esta lectura, mostrando que la mayoría de los valores se concentran entre 20 y 60 g/kg, es decir, predominan los suelos con contenidos bajos a moderados, mientras que los valores más altos son minoritarios y corresponden a zonas puntuales con mejores condiciones edáficas o vegetales.

1.6.2 7.2 Mapa Nitrógeno.

n_pal <- colorNumeric(
  palette = c("#08306B", "#41B6C4", "#FFFF99"),  
  domain = valores_nitrógeno,
  na.color = "transparent")
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_nitrogen, colors = n_pal, opacity = 0.8) %>%
  addLegend(pal = n_pal, values = valores_nitrógeno, title = "Nitrogen [g/kg]")

El contenido total de nitrógeno (g/kg) en el suelo, también a 15–30 cm de profundidad. La paleta azul–verde utilizada indica los siguientes rangos; Los Tonos azul oscuro representan contenidos bajos de nitrógeno, comunes en áreas agrícolas intensivas o zonas con baja acumulación de residuos orgánicos. Mientras que, Tonos verdosos y turquesa reflejan contenidos medios a altos de nitrógeno, generalmente en suelos más fértiles o con vegetación más densa.

Al igual que en el caso del carbono, se aprecia una mayor concentración de nitrógeno en áreas montañosas o de piemonte y una disminución en las zonas planas agrícolas. El histograma del nitrógeno muestra una distribución más asimétrica hacia la derecha, con valores concentrados entre 2 y 5 g/kg, lo que indica que la mayoría de los suelos presentan niveles bajos a moderados de nitrógeno, y solo una fracción pequeña tiene concentraciones elevadas.

1.7 8. Visualización alternativa: velocidad potencial de descomposición de la materia orgánica (M.O.)

El contenido de carbono orgánico del suelo (SOC) y el nitrógeno total (N) son indicadores fundamentales de la calidad y fertilidad del suelo, estrechamente relacionados con la dinámica de la materia orgánica (M.O.). El SOC representa la reserva principal de energía y carbono para los microorganismos del suelo, mientras que el nitrógeno constituye un elemento esencial para su metabolismo y crecimiento.

La relación entre ambos (C/N) permite estimar la velocidad potencial de descomposición de la materia orgánica: Los valores bajos de la relación C/N (menores a 10–12) indican una alta disponibilidad de nitrógeno en relación con el carbono, lo que favorece una rápida descomposición de la materia orgánica y una liberación más eficiente de nutrientes hacia el suelo. En cambio, los valores altos de C/N (superiores a 20) reflejan una acumulación de carbono con escasez de nitrógeno, lo que ralentiza la actividad microbiana y retrasa los procesos de mineralización y descomposición. Por ello, el equilibrio del C/N es un indicador clave del estado biológico y funcional del suelo, determinando su capacidad para mantener la fertilidad y la estabilidad de la materia orgánica

1.7.1 8.1. Creación de un Raster con la Variable de Interes.

cund_cn <- cund_soc / cund_nitrogen
names(cund_cn) <- "C_N_ratio"
terra::writeRaster(cund_cn, "cund_CN_ratio.tif", overwrite = TRUE)

para ello, vamos a atribuir a un objeto la división del aster del carbono orgánico del suelo (SOC) entre el raster del nitrógeno total del suelo (N), celda por celda (cund_cn <- cund_soc / cund_nitrogen);Al dividirlos, se obtiene la relación carbono/nitrógeno (C/N) para cada píxel, que es un indicador de la calidad y del estado de descomposición de la materia orgánica del suelo. el argumento “names(cund_cn) <-”C_N_ratio”” cambia el nombre de la capa dentro del objeto raster a “C_N_ratio”, lo cual facilita su identificación y visualización cuando se maneja en R.(terra::writeRaster(cund_cn, “cund_CN_ratio.tif”, overwrite = TRUE)) exporta un raster resultante a un archivo llamado “cund_CN_ratio.tif” en el directorio de trabajo. Generando así un nuevo archivo GeoTIFF que puedes abrir en R, QGIS o cualquier software de análisis espacial.

1.7.2 8.2. Exploración de Datos.

summary(cund_cn)
##    C_N_ratio     
##  Min.   : 4.595  
##  1st Qu.:11.708  
##  Median :13.440  
##  Mean   :13.948  
##  3rd Qu.:15.531  
##  Max.   :51.211  
##  NA's   :1810

Min. : 4.60 # Representa zonas con una relación C/N muy baja, indicando materia orgánica altamente descompuesta, con abundante nitrógeno disponible para las plantas.

1st Qu. : 11.70 # Un 25% de los valores están por debajo de este número, lo que sugiere que una cuarta parte del territorio tiene una descomposición rápida o moderada.

Median : 13.40 # La mitad de los suelos tienen una relación C/N menor a este valor. Indica una condición intermedia entre descomposición moderada y lenta.

Mean : 13.90 # Promedio general del departamento. Es cercano a la mediana, lo que sugiere una distribución bastante equilibrada sin sesgos fuertes hacia valores extremos.

3rd Qu. : 15.50 # El 75% de los valores se encuentran por debajo de este nivel, lo que indica que solo una cuarta parte del territorio presenta materia orgánica más estable o de descomposición lenta. Max. : 51.20 # Representa puntos aislados con C/N muy alto, posiblemente suelos con gran acumulación de carbono y baja disponibilidad de nitrógeno (por ejemplo, zonas con restos vegetales poco degradados o suelos fríos/húmedos) NA’s :1810 # Número de píxeles sin datos (missing values)

El promedio y la mediana en torno a 13–14 indican que la mayoría de los suelos de Cundinamarca presentan relaciones C/N moderadas, lo cual sugiere una actividad microbiana activa y procesos equilibrados de descomposición de la materia orgánica.

1.7.3 8.3. Histograma Variable (M.O.).

ahora agregamos un histograma par a la nueva variable:

Con la función “terra::hist(cund_cn,)” se crea un histograma a partir de los valores de un objeto raster (SpatRaster), en este caso cund_cn.Cada píxel del raster tiene un valor de C/N, y el histograma muestra cómo se distribuyen esos valores. El argumento (main=) Este argumento define el título principal del gráfico. (xlab = “Relación C/N”) Define la etiqueta del eje X, indicando que los valores que se representan corresponden a la relación entre carbono y nitrógeno. (ylab = “Frecuencia de píxeles”) la etiqueta del eje Y, que muestra la frecuencia, es decir, cuántos píxeles del raster presentan un determinado rango de valores de C/N.

terra::hist(cund_cn,
            main = "Valor modal de la Relación C/N\nProfundidad 15–30 cm - Cundinamarca",
            xlab = "Relación C/N",
            ylab = "Frecuencia de píxeles",
            col = "sienna",
            border = "white")

El histograma de la relación C/N (Carbono/Nitrógeno) para los suelos de Cundinamarca, a una profundidad de 15–30 cm, muestra una distribución donde la mayoría de los píxeles presentan valores entre 10 y 16. Este rango concentra el valor modal (la clase más frecuente), lo que sugiere que la mayor parte de los suelos del departamento mantiene una relación C/N moderada, típica de suelos con un equilibrio relativamente estable entre la cantidad de carbono orgánico y nitrógeno disponible.

Los valores más bajos indican una mayor actividad microbiana y una descomposición más rápida de la materia orgánica, mientras que los valores altos, que aparecen con menor frecuencia, reflejan acumulación de carbono con limitaciones de nitrógeno, lo que puede reducir la velocidad de mineralización. En conjunto, el gráfico sugiere que los suelos de la región presentan una dinámica de descomposición activa, coherente con condiciones de fertilidad intermedia y un buen reciclaje de nutrientes.

1.7.4 8.4. Mapero de la Variable de Interes (M.O.).

Para ello:

atribuimos un objeto “cund_decomp” con la reagrupación de los valores del raster “cund_cn” en tres clases:

  • de o a 10 –> 1 –> Descomposición Rápida (alto N, baja C/N).

  • de 10 a 20 –> 2 –> Descomposición moderada.

  • 20 –> 3 –> Descomposición lenta (alto C, baja disponibilidad de N)

cund_decomp <- terra::classify(cund_cn,
                               rbind(c(0,10,1),
                                     c(10,20,2),
                                     c(20,999,3)))
niveles_decomp <- data.frame(
  ID = c(1, 2, 3),
  Descomposición = c("Lenta", "Moderada", "Rápida"))
  
levels(cund_decomp) <- list(niveles_decomp)

“niveles_decomp <- data.frame(ID = c(1, 2, 3), Descomposición = c(”Lenta”, “Moderada”, “Rápida”))” aqui se construye un data frame asocia el número del objeto (ID) con su etiqueta interpretativa, y se lo asocia a un objeto “niveles_decomp”.

Ahora, para el mapero:

se le atribuye a “colores_decomp” un vector de colores; a cada nombre corresponde a una categoría del raster y cada categoría recibe un color distintivo.

“ggplot() + geom_spatraster(data = cund_decomp)” Dibuja el raster clasificado las tres clases de descomposición, que junto a geom_spatvector(data = cund, …) Superpone los límites municipales de Cundinamarca, con sus excenticidades cómo contorno negro fino (sin relleno).

(scale_fill_manual()) aplica la paleta personalizada “colores_decomp.” y (name=) define el título de la leyenda. (na.value = “transparent”) evita que los valores nulos se muestren con color, que con “guide_legend(reverse = TRUE)” invierte el orden de los niveles en la leyenda (mostrando primero los valores altos o “rápidos”).

El argumento “labs()” Agrega el título, subtítulo y etiquetas de los ejes. Y Para su Estilo Visual, “theme_minimal(base_size = 12)”Aplica el tema minimalista de ggplot2 (sin fondos ni líneas innecesarias) y Establece el tamaño de fuente base para todo el gráfico (ejes, títulos, etc.) en 12 puntos. (Theme=) colora las excentricidades te titulos y leyendas que junto a panel.grid = element_line(color =) coloca las lineas de cuadricula.

colores_decomp <- c(
  "Lenta"    = "#8E5C42",  
  "Moderada" = "#F1C40F",  
  "Rápida"   = "#2ECC71"   
)


ggplot() +
  geom_spatraster(data = cund_decomp) +
  geom_spatvector(data = cund, fill = NA, color = "black", linewidth = 0.3) +
  scale_fill_manual(
    name = "Velocidad de Descomposición de la M.O.",
    values = colores_decomp,
    na.value = "transparent",
    guide = guide_legend(reverse = TRUE)
  ) +
  labs(
    title = "Velocidad Potencial de Descomposición de la Materia Orgánica",
    subtitle = "Relación C/N (15–30 cm) - Departamento de Cundinamarca",
    x = "Longitud", y = "Latitud"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    legend.position = "right",
    legend.title = element_text(face = "bold"),
  )
## <SpatRaster> resampled to 500634 cells.

El mapa de velocidad potencial de descomposición de la materia orgánica en el departamento de Cundinamarca revela tres categorías principales que reflejan la dinámica de los nutrientes en el suelo, directamente relacionadas con la relación Carbono/Nitrógeno (C/N) a profundidades de 15-30 cm. En condiciones de descomposición rápida se caracteriza por relaciones C/N bajas (generalmente inferiores a 12) e indica una intensa actividad microbiana que acelera la mineralización de la materia orgánica. En el territorio cundinamarqués, esta velocidad predominaría probablemente en los valles interandinos más cálidos, particularmente en zonas del Valle del Magdalena, donde las temperaturas más elevadas y la humedad constante favorecen la actividad biológica del suelo. Estas áreas representan suelos con alta fertilidad inmediata y rápida liberación de nutrientes para las plantas.

En condiciones de descomposición lenta, las relaciones C/N más elevadas (superiores a 16), indica limitaciones en la disponibilidad de nitrógeno que restringen la actividad microbiana. En Cundinamarca, estas condiciones se localizarían principalmente en las zonas altoandinas y páramos, donde las bajas temperaturas enlentecen los procesos biológicos, así como en áreas con suelos muy erosionados o con acumulación de materia orgánica recalcitrante. Estos suelos requieren mayor tiempo para liberar nutrientes asimilables por las plantas. Y para una descomposición moderada, Corresponde a relaciones C/N equilibradas (entre 12-16) y constituye el patrón predominante en el departamento. Esta condición refleja un balance óptimo entre el almacenamiento de carbono orgánico y la disponibilidad de nitrógeno, permitiendo una mineralización sostenida sin pérdidas excesivas de nutrientes. Se distribuye ampliamente en las zonas de ladera y mesetas de la cordillera Oriental, donde las condiciones climáticas templadas y los sistemas agrícolas tradicionales mantienen una dinámica estable de nutrientes.

1.8 9. Muestreando el mundo.

Supongamos que también nos interesa obtener una muestra del conjunto de datos SOC, Nitrógeno y MO.

Por ello, utilizamos (set.seed(123456)) que Garantiza que cualquier operación aleatoria posterior (como muestreos) produzca los mismos resultados cada vez que ejecutes el código, y junto a (geog = “+proj=longlat +datum=WGS84”) Define el sistema de referencia de coordenadas geográficas WGS84, “longlat” especifica coordenadas en grados (longitud, latitud) y “datum=WGS84” usa el datum WGS84, estándar global para GPS.

Para Cada variable, ve necesita hacer una reproyección “project()” a cada uno de los rasters de interés , ya que cada uno de ellos estan en proyección IGH, que es una proyección en metros y minimiza distorsiones globales. se pasan los archivos a coordenadas geográficas WGS84 para utilizar grados.

cada uno de los resulltados se les asigna un objeto geog.:

set.seed(123456)
geog ="+proj=longlat +datum=WGS84"
(geog.soc = project(cund_soc, geog))
## class       : SpatRaster 
## size        : 956, 863, 1  (nrow, ncol, nlyr)
## resolution  : 0.00220751, 0.00220751  (x, y)
## extent      : -74.90907, -73.00399, 3.72867, 5.839049  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        :        soc 
## min value   :   9.125371 
## max value   : 276.998291
(geog.nitrogen = project(cund_nitrogen, geog))
## class       : SpatRaster 
## size        : 956, 863, 1  (nrow, ncol, nlyr)
## resolution  : 0.00220751, 0.00220751  (x, y)
## extent      : -74.90907, -73.00399, 3.72867, 5.839049  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        :  nitrógeno 
## min value   :  0.8922093 
## max value   : 12.5459318
(geog.mo = project(cund_cn, geog))
## class       : SpatRaster 
## size        : 956, 863, 1  (nrow, ncol, nlyr)
## resolution  : 0.00220751, 0.00220751  (x, y)
## extent      : -74.90907, -73.00399, 3.72867, 5.839049  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : C_N_ratio 
## min value   :  4.823197 
## max value   : 52.883118

1.8.1 9.1. Procesamiento de Datos.

samples_soc <- terra::spatSample(geog.soc, 2000, method = "random", as.points = TRUE)
samples_nitrogen <- terra::spatSample(geog.nitrogen, 2000, method = "random", as.points = TRUE)
samples_mo <- terra::spatSample(geog.mo, 2000, method = "random", as.points = TRUE)

Lo que busca el codigo, es qué para cada Raster,Seleccionar 2000 ubicaciones aleatorias en el área de Cundinamarca,distribuidas uniformemente por todo el territorio, para extraer los valores de raster para cada ubicación.

(muestras_soc <- sf::st_as_sf(samples_soc) %>% st_as_sf())
## Simple feature collection with 2000 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.90576 ymin: 3.729774 xmax: -73.00509 ymax: 5.837946
## Geodetic CRS:  GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]
## First 10 features:
##         soc                   geometry
## 1  32.95347 POINT (-74.23909 5.707703)
## 2  27.11556 POINT (-74.42011 5.747438)
## 3  34.82338 POINT (-73.41128 5.562007)
## 4        NA POINT (-74.86602 5.316973)
## 5  45.38998 POINT (-74.36934 5.336841)
## 6  40.50652 POINT (-74.36271 5.676797)
## 7  42.67120 POINT (-73.98744 4.992469)
## 8  32.49107 POINT (-73.53269 5.067525)
## 9  49.98228 POINT (-74.43556 4.376574)
## 10 61.45358 POINT (-73.88589 4.164653)
(muestras_nitrogen <- sf::st_as_sf(samples_nitrogen) %>% st_as_sf())
## Simple feature collection with 2000 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.90797 ymin: 3.729774 xmax: -73.0073 ymax: 5.837946
## Geodetic CRS:  GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]
## First 10 features:
##    nitrógeno                   geometry
## 1   2.358443  POINT (-74.17949 5.70108)
## 2   2.664014 POINT (-74.80201 3.952732)
## 3   3.447291  POINT (-73.9168 4.374367)
## 4         NA POINT (-73.02055 3.983638)
## 5   5.088552   POINT (-74.3583 4.45163)
## 6   1.192616 POINT (-74.87706 4.707701)
## 7   1.819209 POINT (-74.16183 5.422934)
## 8   2.799183 POINT (-74.53711 5.241918)
## 9   4.087555 POINT (-73.93004 3.873262)
## 10  1.603402 POINT (-74.41349 4.626023)
(muestras_mo <- sf::st_as_sf(samples_mo)%>% st_as_sf())
## Simple feature collection with 2000 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.90797 ymin: 3.729774 xmax: -73.00509 ymax: 5.837946
## Geodetic CRS:  GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]
## First 10 features:
##    C_N_ratio                   geometry
## 1   11.35209   POINT (-73.153 4.767303)
## 2   12.35113 POINT (-74.13534 4.873264)
## 3   16.99043 POINT (-74.65852 4.065315)
## 4   10.72661 POINT (-73.33622 5.731985)
## 5   15.07742 POINT (-74.39803 5.703288)
## 6         NA POINT (-74.66956 5.445009)
## 7   12.06835 POINT (-73.34064 5.215428)
## 8         NA POINT (-74.85719 3.760679)
## 9   16.51412 POINT (-74.84836 5.197768)
## 10  14.34747 POINT (-74.84836 5.303728)

Y para cada dato, se eliminan los valores que no reistran, y se atribuyen a nuevos objetos.

nmuestras_soc <- na.omit(muestras_soc)
nmuestras_nitrogen <- na.omit(muestras_nitrogen)
nmuestras_mo <- na.omit(muestras_mo)

Ahora, visualicemos las muestras limpias:

longit_soc <- st_coordinates(nmuestras_soc)[,1]
lati_soc <- st_coordinates(nmuestras_soc)[,2]
soc <- nmuestras_soc$soc


longit_nitrogen <- st_coordinates(nmuestras_nitrogen)[,1]
lati_nitrogen <- st_coordinates(nmuestras_nitrogen)[,2]
N <- nmuestras_nitrogen$nitrógeno 


longit_mo <- st_coordinates(nmuestras_mo)[,1]
lati_mo <- st_coordinates(nmuestras_mo)[,2]
C_N_ratio <- nmuestras_mo$C_N_ratio 

Para cada raster “st_coordinates()” Extrae las coordenadas de objetos espaciales sf y devuelve una matriz con columnas de longitud (X) y latitud (Y) y asigna una columna para longitudes y latitudes; y nmuestras_x$x extrae la columna soc del data frame espacial y ccede a una columna específica por nombre.

como se puede apreciar:

summary(soc)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10.43   25.78   38.33   39.57   51.94   94.20
summary(N)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.064   1.664   3.068   2.940   3.903   7.598
summary(C_N_ratio)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   6.639  11.862  13.513  13.927  15.455  34.657

Estos códigos length() verifican la cantidad de observaciones/muestras que tienes en cada variable después de la limpieza.

length(soc)
## [1] 1894
length(N)
## [1] 1905
length(C_N_ratio)
## [1] 1892

Este código realiza un proceso completo de preparación y verificación de datos para asegurar la calidad y consistencia de tus muestras.

id <- seq(1, 1905, 1) 

cat("Longitud de los vectores originales:\n")
## Longitud de los vectores originales:
cat("longit_soc:", length(longit_soc), "\n")
## longit_soc: 1894
cat("longit_nitrogen:", length(longit_nitrogen), "\n")
## longit_nitrogen: 1905
cat("longit_mo:", length(longit_mo), "\n")
## longit_mo: 1892
min_length <- min(
  length(longit_soc), length(lati_soc), length(soc),
  length(longit_nitrogen), length(lati_nitrogen), length(N),
  length(longit_mo), length(lati_mo), length(C_N_ratio),
  1905
)


id <- seq(1, min_length, 1)

# Para SOC
sitios_soc <- data.frame(
  id = id,
  longit = longit_soc[1:min_length],
  latit = lati_soc[1:min_length],
  soc = soc[1:min_length]
)


sitios_N <- data.frame(
  id = id,
  longit = longit_nitrogen[1:min_length],
  latit = lati_nitrogen[1:min_length],
  N = N[1:min_length]
)


sitios_MO <- data.frame(
  id = id,
  longit = longit_mo[1:min_length],
  latit = lati_mo[1:min_length],
  C_N_ratio = C_N_ratio[1:min_length]
)


head(sitios_soc)
##   id    longit    latit      soc
## 1  1 -74.23909 5.707703 32.95347
## 2  2 -74.42011 5.747438 27.11556
## 3  3 -73.41128 5.562007 34.82338
## 4  4 -74.36934 5.336841 45.38998
## 5  5 -74.36271 5.676797 40.50652
## 6  6 -73.98744 4.992469 42.67120
head(sitios_N)
##   id    longit    latit        N
## 1  1 -74.17949 5.701080 2.358443
## 2  2 -74.80201 3.952732 2.664014
## 3  3 -73.91680 4.374367 3.447291
## 4  4 -74.35830 4.451630 5.088552
## 5  5 -74.87706 4.707701 1.192616
## 6  6 -74.16183 5.422934 1.819209
head(sitios_MO)
##   id    longit    latit C_N_ratio
## 1  1 -73.15300 4.767303  11.35209
## 2  2 -74.13534 4.873264  12.35113
## 3  3 -74.65852 4.065315  16.99043
## 4  4 -73.33622 5.731985  10.72661
## 5  5 -74.39803 5.703288  15.07742
## 6  6 -73.34064 5.215428  12.06835
str(sitios_soc)
## 'data.frame':    1892 obs. of  4 variables:
##  $ id    : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ longit: num  -74.2 -74.4 -73.4 -74.4 -74.4 ...
##  $ latit : num  5.71 5.75 5.56 5.34 5.68 ...
##  $ soc   : num  33 27.1 34.8 45.4 40.5 ...
summary(sitios_soc)
##        id             longit           latit            soc       
##  Min.   :   1.0   Min.   :-74.90   Min.   :3.730   Min.   :10.43  
##  1st Qu.: 473.8   1st Qu.:-74.40   1st Qu.:4.248   1st Qu.:25.84  
##  Median : 946.5   Median :-73.93   Median :4.790   Median :38.34  
##  Mean   : 946.5   Mean   :-73.94   Mean   :4.787   Mean   :39.59  
##  3rd Qu.:1419.2   3rd Qu.:-73.49   3rd Qu.:5.335   3rd Qu.:51.95  
##  Max.   :1892.0   Max.   :-73.01   Max.   :5.838   Max.   :94.20

(id <- seq(1, 1905, 1)) crea un vector de identificadores del 1 al 1905, Cada muestra tendrá un ID único para las observaciones.Y con cada comando “cat” mprime en la consola cuántos elementos tiene cada vectorn ente ellos:

*longit → longitudes (x)

*latit → latitudes (y)

*soc → contenido de carbono orgánico del suelo

*N → nitrógeno

*C_N_ratio → relación carbono/nitrógeno

“min_length <-” busca el menor tamaño entre todos los vectores y el número 1900. Así, si alguno tiene menos de 1905 datos, se ajusta todo a ese tamaño y se recorta cada vector para que todos tengan exactamente min_length elementos. Luego, vuelve a crear id con esa misma longitud y se crean 3 dataframes para cada variable.

sitios <- data.frame(
  longit = sitios_soc$longit,
  latit = sitios_soc$latit,
  soc = sitios_soc$soc,
  N = sitios_N$N,
  C_N_ratio = sitios_MO$C_N_ratio
)


m <- leaflet() %>%
  addTiles() %>%  # Corrección: addTiles() no addriles()
  addMarkers(
    lng = sitios$longit,
    lat = sitios$latit,  # Corrección: latit no lati
    popup = paste(
      "SOC: ", sitios$soc, "<br>",
      "N: ", sitios$N, "<br>",
      "C/N Ratio: ", sitios$C_N_ratio, "<br>"
    ),
    clusterOptions = markerClusterOptions()  # Corrección: clusterOptions con O mayúscula
  )
m

“sitios <- data.frame ()” une toda la información de los tres data frames anteriores (sitios_soc, sitios_N, sitios_MO) en uno solo. Cada fila representa un punto geográfico de muestreo con sus tres variables químicas. “leaflet()” crea un objeto de mapa vacío. “addTiles()” añade el fondo de mapa base.Y “addMarkers()” incluye marcadores para variables:

*lng → coordenadas de longitud (eje X, este-oeste).

*lat → coordenadas de latitud (eje Y, norte-sur).

*popup → texto que aparece al hacer clic en cada punto.

  • paste() para combinar etiquetas y valores.

*clusterOptions = markerClusterOptions() → agrupa los marcadores cercanos automáticamente (útil si hay muchos puntos).bEn vez de ver miles de puntos sobrepuestos, se verán burbujas agrupadas que se separan al hacer zoom.

1.9 10. Guardar Muestras.

sf::st_write(muestras_soc, "data/soc_samples.gpkg", driver = "GPKG", append=FALSE)
## Deleting layer `soc_samples' using driver `GPKG'
## Writing layer `soc_samples' to data source `data/soc_samples.gpkg' using driver `GPKG'
## Writing 2000 features with 1 fields and geometry type Point.
sf::st_write(muestras_nitrogen, "data/nitrogen_samples.gpkg", driver = "GPKG", append=FALSE)
## Deleting layer `nitrogen_samples' using driver `GPKG'
## Writing layer `nitrogen_samples' to data source 
##   `data/nitrogen_samples.gpkg' using driver `GPKG'
## Writing 2000 features with 1 fields and geometry type Point.
sf::st_write(muestras_mo, "data/cn_ratio_samples.gpkg", driver = "GPKG", append=FALSE)
## Deleting layer `cn_ratio_samples' using driver `GPKG'
## Writing layer `cn_ratio_samples' to data source 
##   `data/cn_ratio_samples.gpkg' using driver `GPKG'
## Writing 2000 features with 1 fields and geometry type Point.

sf::st_write() → Función del paquete sf utilizada para exportar datos espaciales (puntos, polígonos, etc.) a un archivo.

Primer argumento: el objeto que se desea escribir (p. ej., muestras_soc).

Segundo argumento: la ruta del archivo de salida (p. ej., data/soc_samples.gpkg).

driver = "GPKG" → especifica el formato de salida: GeoPackage, un formato moderno que permite almacenar capas vectoriales, atributos, sistemas de coordenadas, etc.

Cada línea crea (o reemplaza) un archivo GeoPackage dentro de la carpeta data/.

El parámetro append controla qué sucede si el archivo ya existe.

append = TRUE → agrega los nuevos datos al archivo existente (conserva las capas antiguas o agrega entidades a la misma capa).

append = FALSE → sobrescribe el archivo o la capa existente; comienza desde cero.

1.10 11.Conclusiones.

Los resultados revelan que los suelos del departamento presentan, en general, condiciones de fertilidad intermedia, con valores promedio de carbono orgánico de 39.65 g/kg y nitrógeno total de 2.935 g/kg en la capa de 15-30 cm de profundidad. La relación C/N promedio de 13.95 indica una dinámica de descomposición moderada en la mayoría del territorio, sugiriendo un equilibrio adecuado entre la acumulación de materia orgánica y la disponibilidad de nitrógeno para los procesos biológicos. Este equilibrio es fundamental para mantener la fertilidad del suelo y la sostenibilidad de los sistemas productivos.

La distribución espacial de estas variables muestra una clara diferenciación entre las zonas montañosas y las áreas de valle, reflejando la influencia de factores como la topografía, el clima y posiblemente las prácticas de uso del suelo. Las áreas con mayores elevaciones tienden a presentar mayores contenidos de carbono orgánico y relaciones C/N más elevadas, mientras que las zonas planas, sometidas a mayor intensidad de uso agrícola, muestran valores generalmente más bajos.

La metodología implementada, que combina el uso de datos globales de SoilGrids con herramientas de análisis espacial en R, ha demostrado ser efectiva para la caracterización de propiedades del suelo a escala departamental. Los mapas generados y el muestreo sistemático realizado constituyen una base valiosa para la toma de decisiones.

1.11 12. Referencias.

cuaderno origial por : Lizarazo, I. 2025. How to download soil data from SoilGrids. Available at: https://rpubs.com/ials2un/downloadSoilGrids

+Oklahoma State University Extension. (s.f.). Building Soil Organic Matter for a Sustainable Organic Crop Production. Recuperado de https://extension.okstate.edu/fact-sheets/building-soil-organic-matter-for-a-sustainable-organic-crop-production.html

+Lang, D. T. (2023). XML: Tools for parsing and generating XML within R and S-Plus [Software]. R package version 3.99-0.14. https://CRAN.R-project.org/package=XML +Pebesma, E. (2018). sf: Simple Features for R [Software]. R package version 1.0-14. https://CRAN.R-project.org/package=sf +Bivand, R. (2023). gdalUtilities: Wrappers for GDAL Utilities [Software]. R package version 1.0.1. https://CRAN.R-project.org/package=gdalUtilities +Hijmans, R. J. (2023). terra: Spatial Data Analysis [Software]. R package version 1.7-39. https://CRAN.R-project.org/package=terra +Wickham, H., François, R., Henry, L., Müller, K., & Vaughan, D. (2023). dplyr: A Grammar of Data Manipulation [Software]. R package version 1.1.3. https://CRAN.R-project.org/package=dplyr +Cheng, J., Karambelkar, B., & Xie, Y. (2023). leaflet: Create Interactive Web Maps with the JavaScript ‘Leaflet’ Library [Software]. R package version 2.1.2. https://CRAN.R-project.org/package=leaflet +Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, T. L., Miller, E., Bache, S. M., Müller, K., Ooms, J., Robinson, D., Seidel, D. P., Spinu, V., … Yutani, H. (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686. https://doi.org/10.21105/joss.01686 +Hernangomez, D. (2023). tidyterra: tidyverse Methods and ggplot2 Helpers for terra Objects [Software]. R package version 0.4.1. https://CRAN.R-project.org/package=tidyterra +R Core Team. (2023). R: A language and environment for statistical computing [Software]. R Foundation for Statistical Computing. https://www.R-project.org/

hecho por: +BOVVER+